Martin Dulac

Armelle Melet

Erwan Galardon

Reversible Detection and Quantification of Hydrogen Sulfide by Fluorescence using the Hemoglobin I from Lucina Pectinata

Keywords:

ABSTRACT: A new detection system for the endogenous gaseous transmitter and environmental pollutant hydrogen sulfide is presented. It is based on the modulation of the fluorescence spectrum of a coumarin dye by the absorption spectrum of the recombinant hemoglobin I from clam Lucina Pectinata upon coordination of the analyte. While we establish that the reported affinity of rHbI for H 2 S has been overestimated, the association of the protein with an appropriate fluorophore allows fast, easy and reversible detection and quantification of hydrogen sulfide in buffer as well as biological fluids such as human plasma, with a quantification limit around 200 nM at pH 7.4.

The demonstration by Snyder et al. [START_REF] Yang | H2S as a physiologic vasorelaxant: hypertension in mice with deletion of cystathionine gamma-lyase[END_REF] that hydrogen sulfide (H 2 S) is endogenously produced and acts as a vasodilator in mammals was a turning point in the history of this gas in biology. Indeed, this work confirmed the hypothesis that emerged in the late 90's that H 2 S could be the third gaseous transmitter, along with nitric oxide NO and carbon monoxide CO. [START_REF] Abe | The possible role of hydrogen sulfide as an endogenous neuromodulator[END_REF] It also boosted the development of techniques and molecular tools for studying the potential therapeutic role of H 2 S, and its signaling functions. Thus, numerous hydrogen sulfide donors mimicking the slow continuous biosynthesis of hydrogen sulfide were designed. [START_REF] Zhao | Cysteine-activated hydrogen sulfide (H2S) donors[END_REF][START_REF] Powell | A review of hydrogen sulfide (H2S) donors: Chemistry and potential therapeutic applications[END_REF] Simultaneously, various methods for measuring H 2 S in biological media were developed. At first, these methods were mostly derived from procedures used in environmental chemistry, since sulfur is a major physicochemical and geological indicator and hydrogen sulfide a natural and anthropogenic pollutant, but they were not appropriate for biological studies. For instance, these methods resulted in aberrant concentrations of hydrogen sulfide in biological fluids. [START_REF] Whiteman | Hydrogen sulfide and the vasculature: a novel vasculoprotective entity and regulator of nitric oxide bioavailability ?[END_REF] Later, more relevant techniques, using fluorescence, polarography or chromatography, were developed to detect or quantify H 2 S in biological samples. However, each offers its own advantages and disadvantages. [START_REF] Olson | A Practical Look at the Chemistry and Biology of Hydrogen Sulfide[END_REF][START_REF] Hartle | A practical guide to working with H2S at the interface of chemistry and biology[END_REF] For instance, chromatography displays a detection limit ~ 1 nM, but requires a derivatization step and a chromatographic separation, making it time consuming but well suited for analytical studies requiring high sensitivity. Polarographic sensors allow for real-time detection of hydrogen sulfide, [START_REF] Xu | Electrochemical hydrogen sulfide biosensors[END_REF] but they require costly apparatus. Last, fluorescent sensors are nowadays the best choice for bioimaging H 2 S. [START_REF] Guo | Fluorescence chemosensors for hydrogen sulfide detection in biological systems[END_REF] However, they are not suitable for real time detection of H 2 S, because of their reaction kinetics and lack of reversibility. Indeed, most specific fluorescent sensors are slow "reaction based" probes, which consume the analyte. Typical strategies include the reduction of azides [START_REF] Lippert | Reaction-Based Fluorescent Probes for Selective Imaging of Hydrogen Sulfide in Living Cells[END_REF] or nitro compounds [START_REF] Montoya | Selective turn-on fluorescent probes for imaging hydrogen sulfide in living cells[END_REF] , the use of the dualnucleophilicity of H 2 S, [START_REF] Liu | Capture and Visualization of Hydrogen Sulfide by a Fluorescent Probe[END_REF] or the formation of metal sulfides. [START_REF] Sasakura | Development of a Highly Selective Fluorescence Probe for Hydrogen Sulfide[END_REF] On another hand, several approaches using non-covalent interactions are currently explored to develop reversible detection systems. For instance, the reversible adsorption/desorption of H2S on molybdenum oxide was used to detect gaseous hydrogen sulfide. [START_REF] Woo | sensitive, and reversible detection of H2S using Mo-doped ZnO nanowire network sensors[END_REF] Pluth et al. recently reported a UV-absorbing supramolecular receptor that binds the hydrosulfide anion in organic solvents. Additionally, we and others are aiming at the development of reversible coordination-based fluorescent sensors for H2S in aqueous buffered solutions. However, these systems so far suffer from moderate affinities and / or detection limits. 15, 16, 17 Here, we present a new selective and reversible fluorescent system based on the "absorption into fluorescence", a strategy recently introduced by D'Auria et. al., [START_REF] Strianese | Absorption into fluorescence. A method to sense biologically relevant gas molecules[END_REF] to selectively detect hydrogen sulfide in real time with competitive detection limits of ~ 200 nM at pH 7.4. This strategy consists in using a sensor with a high molar attenuation coefficient to modulate the fluorescence properties of a fluorophore. To be efficient, the absorption band of the sensor should overlap with the excitation band of the fluorophore, and it should shift upon reaction with the analyte. As such, we selected the most avid biological system for hydrogen sulfide as a sensor, i.e. the hemoglobin HbI for the clam Lucina Pectinata. [START_REF] Kraus | Hemoglobins of the Lucina-Pectinata Bacteria Symbiosis .1. Molecular-Properties, Kinetics and Equilibria of Reactions with Ligands[END_REF] The unique affinity of HbI for H 2 S lies in the peculiar layout of its distal active site which contains a glutamine and two phenylalanines at the E7, B10 and E11 positions, respectively. [START_REF] Rizzi | Structure of the sulfide-reactive hemoglobin from the clam Lucina pectinata. Crystallographic analysis at 1.5 A resolution[END_REF] These unusual structural determinants allow for fast binding and slow dissociation of the protonated form H 2 S to the ferric heme center of HbI. [START_REF] Kraus | Hemoglobins of the Lucina-Pectinata Bacteria Symbiosis .1. Molecular-Properties, Kinetics and Equilibria of Reactions with Ligands[END_REF][START_REF] Pietri | Factors Controlling the Reactivity of Hydrogen Sulfide with Hemeproteins[END_REF] Despite the high affinity displayed by HbI for hydrogen sulfide (K' d † = 3.4 nM), [START_REF] Kraus | Hemoglobins of the Lucina-Pectinata Bacteria Symbiosis .1. Molecular-Properties, Kinetics and Equilibria of Reactions with Ligands[END_REF] its use has so far been limited to the determination of H 2 S solubility by spectrophotometry, [START_REF] Boffi | Determination of H2S solubility via the reaction with ferric hemoglobin I from the bivalve mollusc Lucina pectinata[END_REF] or to the tentative attempt to prepare a bioelectrode. [START_REF] Díaz-Ayala | Engineered (Lys)6-Tagged Recombinant Sulfide-Reactive Hemoglobin I for Covalent Immobilization at Multiwalled Carbon Nanotubes[END_REF] In addition to its high affinity for H 2 S, HbI is also perfectly suited for the biomolecular gating approach reported therein that necessitates a strong chromophore. Indeed, upon the coordination of H 2 S, the observed red-shift of the Soret band (ε ~ 178 000 L.mol -1 .cm -1. ) from 407 to 425 nm should results in a measurable modulation of the photophysical properties of a fluorophore with an appropriate excitation around 407 nm.

EXPERIMENTAL SECTION

Materials. Restriction enzymes and Quick-Ligation kits were purchased from New England Biolabs. Terrific broth medium (Sigma EZmixT9179) and hemin (Sigma 51280) were purchased from Sigma-Aldrich. Kanamycin and Isopropyl β-D-1thiogalactopyranoside (IPTG) were obtained from Euromedex. The BLi5 strain was a gift from Jacques Ballalou (Institut Pasteur, Paris, France). The Protein 200 Plus reagent kits and chips were provided by Biorad. Ultrafiltration devices Amicon Ultra-15 were obtained from Millipore. Anhydrous sodium hydrosulfide (NaSH) was purchased from Alfa Aesar, and kept in a glove box until use.

Solutions were prepared in buffer A previously degassed by purging with argon. The synthesis of N-[2(2-Hydroxyethoxy)ethyl]-6,8-difluoro-7hydroxycoumarin-3-carboxamide was carried out in DMF from pacific blue succidimidyl ester (ThermoFisher scientific), 2-(2aminoethoxy)ethanol and triethylamine (Aldrich) as previously described. [START_REF] Kerkovius | A Practical Synthesis of 6,8-Difluoro-7-hydroxycoumarin Derivatives for Fluorescence Applications[END_REF] All reactions were carried out in buffer A: potassium phosphate buffer (KPi) (50 mM, pH 7.4 supplemented with 1 mM of diethylenetriaminepentaacetic acid (dtpa) as a chelating agent), unless stated otherwise.

Cloning of HbI ORF into the pET28a expression vector. The cDNA sequence of HbI from Lucina Pectinata [START_REF] Antommattei-Perez | The cDNA-derived amino acid sequence of hemoglobin I from Lucina pectinata[END_REF] was obtained from Genbank database (GenBank: AF187049.1) and the rare codons changed to optimize expression in host E. coli. A NcoI restriction site was added in frame at the 5' end to facilitate cloning. The modified sequence was synthesized and cloned into a pUC57 vector by Proteogenix (Schiltigheim, France). The pUC57-HbI vector was then digested by NcoI / BamHI, and the resulting gel-purified inserted sequence NcoI/BamHI HbI cDNA quick-ligated into the pET28a expression vector previously digested by NcoI / BamHI, and dephosphorylated. The resulting pET28a-HbI vector was used to over-express the recombinant rHbI. Expression of recombinant rHbI in E. coli BLi5 strain. The culture protocol was derived from Leon et al. [START_REF] Leon | High-level production of recombinant sulfide-reactive hemoglobin I from Lucina pectinata in Escherichia coli -High yields of fully functional holoprotein synthesis in the BLi5 E-coli strain[END_REF] with modifications as follows. The vector pET28a-HbI was transformed into BLi5 E. coli strain for protein expression. An overnight culture was added at 1 / 50 volume to Terrific broth medium supplemented with Kanamycin 100 µg/mL and chloramphenicol 30 µg/mL. Bacteria were grown at 37°C until OD 600 = 2. After lowering the temperature to 30°C, hemin (30 µg/mL) and IPTG 1 mM were added to induce rHbI expression. After 16 h culture at 30°C and 150 rpm, bacteria were pelleted and stored at -80 °C until use.

Purification of recombinant rHbI. Bacterial pellets were thawed and re-suspended in 50 mM Hepes buffer (pH 7.5) containing Complete TM protease inhibitor cocktail (Roche) and DNaseI (60 µg /mL). After sonication of the solution with a Vibra-cell sonicator (Fisher Bioblock, Vibra-Cell 75115, 500 W) at 30 % amplitude for 6 min on ice (3 s on, 15 s off), the solution was centrifuged (18 000 rpm; 37 600 g for 45 min), and the clarified supernatant was then subjected to ammonium sulfate precipitation (0→55→90 %) and then centrifuged (15 000 rpm; 14 000 g) for 20 min at 4 °C. The 55→90% pellet was re-suspended and dialyzed overnight against 20 mM Tris.HCl buffer (pH 8.5) at 4°C, followed by two buffer changes. The dialysate was incubated in batch with Q-sepharose resin (1 mL for 15 mg proteins) for 1 h 30 at 4°C. The Q-sepharose resin was then washed with 3 volumes of 20 mM Tris-HCl buffer (pH 8.5) at 4°C and the flowthrough containing rHbI was concentrated with an Amicon ultra-15 mL filter (3 kDa cut-off). The solution was then loaded on a HiLoad 16/60 Superdex-75 column (GE Healthcare) equilibrated with 25 mM Hepes buffer (pH 7.5), 150 mM KCl, at a flow rate of 0.3mL/min. After analysis by microfuidic electrophoresis with Experion (Biorad), fractions containing pure rHbI were pooled, concentrated with an Amicon ultra-15 mL filter, and stored as aliquots at -80°C. Typical yields were 30-40 mg HbI per L of culture. According to Experion analysis (Biorad), the purity of rHbI was ≥95 %. The molecular mass of 15835 Da, assessed with the calibrated Superdex-75 was close to the theoretical value for the monomeric holoprotein (15429 Da). The molecular mass determined by ESI-MS on an LTQ Orbitrap instrument (14812.4 Da) confirmed to be in good agreement with the theoretical value for the apoprotein (14812.8 Da). Physical measurements. UV-vis spectra were recorded on a Shimadzu UV-2700 or Carry 300 at 25°C using 1 mL or 100 µL quartz cuvettes. Stopped-flow kinetic experiments were performed on a Biologic SFM-300 at 303 K. Fluorescence spectra were recorded on a Hitachi F-7000 spectrometer at 25°C using 100 µL quartz cuvettes. Excitation wavelength at 410 nm and emission wavelength at 455 nm were used. All kinetics and thermodynamic parameters were obtained by fitting the experimental data with the appropriate equation using SigmaPlot.

Kinetics studies. A) For stopped flow studies, two syringes (S1 & S2

) were used, with the following contents: S2: NaSH at various concentrations in buffer A, S1: 6 µM solution of rHbI in buffer A. Each run consisted of mixing 1/1 parts of S1/S2. An average of at least three runs was recorded for each substrate concentration. B) The determination of the association rate constant k' on for the fixation of the cyanide (0.025-0.375 mM) to rHbI was determined at 25°C by UV/Vis spectroscopy by monitoring the change in absorbance at 420 nm as a function of time. The data were fitted with an exponential equation to obtain the k' obs value at each concentration of cyanide. The k' obs values were then plotted against the concentrations of potassium cyanide to get the k' on and k off value by fitting the data with a linear equation. C) The determination of the k off of rHbI-H 2 S was carried out by UV/Vis spectroscopy. Typically: 5 µL of a 1.2 mM solution of rHbI were mixed in 935 µL of buffer A with 60 µL of a 1 mM solution of NaSH. A first UV/Vis spectrum was recorded to confirm the formation of rHbI-H 2 S, then 10 µL of a 2 M solution of potassium cyanide (20 mM final) were added to the reaction mixture and the displacement reaction was monitored by UV/Vis is spectroscopy. Fluorescence studies. A typical experiment is as follows: to a 100 µL cuvette were added 97 µL of buffer A, 1 µL of a 100 µM solution of PB in DMSO (1 µM final), 2 µL of a 1 mM solution of rHbI in buffer A (20 µM final). The average value of three recordings was used as the initial fluorescence value (F 0 ). For each concentration, the average of three measurements was used to determine the ratio F/F 0 . Each experiment was independently repeated 3 times to determine the S.D.. Plasma studies. Plasma from whole blood samples from healthy blood bank donors ("Etablissement Français du Sang", Convention # 07/CABANEL/106; Paris; France) was obtained by density centrifugation.

RESULTS

Characterization of rHbI and reactivity with H 2 S. The recombinant protein rHbI was produced in the BLi5 E. coli strain, after optimization of the protocol first described by López-Garriga et al. [START_REF] Leon | High-level production of recombinant sulfide-reactive hemoglobin I from Lucina pectinata in Escherichia coli -High yields of fully functional holoprotein synthesis in the BLi5 E-coli strain[END_REF] Derived from the BL21(DE3), this strain harbors a plasmid pDIA17 [START_REF] Munier | Isolation and Characterization of Catalytic and Calmodulin-Binding Domains of Bordetella-Pertussis Adenylate-Cyclase[END_REF] coding for the lacIq repressor, which maintains the recombinant protein expression strongly repressed until it is induced by IPTG. In addition, the BLi5 strain also efficiently incorporates hemin from the culture medium. [START_REF] Leon | High-level production of recombinant sulfide-reactive hemoglobin I from Lucina pectinata in Escherichia coli -High yields of fully functional holoprotein synthesis in the BLi5 E-coli strain[END_REF] Briefly, after optimization of the cDNA sequence [START_REF] Antommattei-Perez | The cDNA-derived amino acid sequence of hemoglobin I from Lucina pectinata[END_REF] (mainly the codon usage), the later was synthesized and cloned into a pU57 vector. Following a subcloning step into a pET28a expression vector, the BLi5 E. coli strain was transformed with the resulting pET28a-rHbI vector and rHBI protein expression was induced by addition of IPTG in the presence of hemin. The oxy-ferrous form of the protein, which displays the expected mass for the monomeric rHbI (see Experimental part), was then obtained after salting out and two chromatographic steps with typical yields of 30-40 mg rHbI per L of culture and a purity > 95%. The ferric form of rHbI was obtained upon oxidation with ferricyanide of the purified oxy-ferrous form of rHbI. It displays a UV/Vis spectrum identical to the one reported for the wild type HbI isolated from the clam or the recombinant HbI. [START_REF] Kraus | Hemoglobins of the Lucina-Pectinata Bacteria Symbiosis .1. Molecular-Properties, Kinetics and Equilibria of Reactions with Ligands[END_REF][START_REF] Leon | High-level production of recombinant sulfide-reactive hemoglobin I from Lucina pectinata in Escherichia coli -High yields of fully functional holoprotein synthesis in the BLi5 E-coli strain[END_REF] Upon addition of an excess of sodium hydrosulfide, the initial Soret and Q bands of the met-aquo form of rHbI (407, 502 and 629 nm) shift to new values (359, 425, 545 and 575 nm) corresponding to the met-hydrogensulfide form rHbI-H 2 S (Figure S1). The formation of the later was also corroborated by resonance Raman spectroscopy (data not shown), which indicates a spin change from 5/2 to 1/2 (the marker v 3 shifts from 1483 cm -1 to 1500 cm -1 ) and the formation of a Fe-S bond (v Fe-S = 374 cm -1 ). [START_REF] Cerda | Resonance Raman Studies of the Heme-Ligand Active Site of Hemoglobin I from Lucina pectinata[END_REF] To confirm that the binding properties of rHbI are identical to the ones displayed by wtHbI, we determined the relative kinetic constants k' on † and k off for its interaction with H2S. Stopped-flow absorbance experiments performed under pseudo-first order conditions allowed the determination of a k' on = (57.6 ± 2.0) x 10 3 M -1 .s -1 (Figure S2). Monitoring of the exponential decay (or rise) of the absorbance at 425 (or 407) nm after dilution of a concentrated solution of rHbI-H 2 S (Figure S3) gave k off values of (2.5 ± 0.1) x 10 -4 s -1 and (2.1 ± 0.1) x 10 -4 s -1 , respectively. These data correlate very well with those reported for the wild type enzyme (k' on ~ 65 x 10 3 M -1 .s -1 and k off = 2.2 x 10 -4 s -1 , respectively). [START_REF] Kraus | Hemoglobins of the Lucina-Pectinata Bacteria Symbiosis .1. Molecular-Properties, Kinetics and Equilibria of Reactions with Ligands[END_REF] Calculation of the dissociation constant from these values gives K' d values of 3.8 nM and 3.4 nM for rHbI and wtHbI, respectively. It has been shown that the determination of the dissociation constant under thermodynamic conditions by a titration experiment resulted in a K' d = 90 nM for wtHbI. [START_REF] Kraus | Hemoglobins of the Lucina-Pectinata Bacteria Symbiosis .1. Molecular-Properties, Kinetics and Equilibria of Reactions with Ligands[END_REF] We have performed a similar experiment with rHbI, during which we followed the increase in absorbance at 425 nm or the decrease at 407 nm observed upon addition of sodium hydrosulfide to rHbI (Figure 1). The corresponding data were fitted with the quadratic equation (1) (Figure S4), giving rise to K' d values of 332 ± 32 and 347 ± 59 nM at 425 and 407 nm, respectively. These values are within the same range as the 90 nM reported for wtHbI, although larger, most probably because of different experimental settings and of the use of sodium hydrosulfide instead of the hydrogen sulfide gas. Noteworthy, these values strongly differ from the low nM affinities obtained from the kinetic data.

Figure 2. Left: UV/Vis spectra of an 11 µM solution of rHbI in

Buffer A with 110 µM NaSH, recorded before (red) and immediately after (blue) dilution. For comparison purpose, the UV/Vis spectrum of the later curve has been multiplied by 10. Right: monitoring of the evolution of the UV/Vis spectra of a 6 µM solution of rHbI-H 2 S (from 6 µM rHbI and 60 µM NaSH) after addition of a 20 mM solution of potassium cyanide. Inset: monitoring of the absorbance at λ = 420 nm over time.

At first, we assumed that the origin of the aforementioned discrepancies originated from an overestimated assessment of the k off . Accordingly, we decided to focus on the dissociation step, and to begin with, we recorded the UV/Vis spectra of a solution of rHbI-H 2 S before and rapidly after its dilution. The results shown in Figure 2 (left) indicate that approximately 15% of rHbI-H 2 S has been converted to rHbI within the mixing time, which is clearly inconsistent with the half-life of 3150 s determined from the dissociation rate constant k off . Based on this experiment, we posited that a displacement experiment in the presence of cyanide as the competing ligand would be appropriate to obtain a more definite dissociation rate constant value. [START_REF] Dou | Mechanism of hydrogen cyanide binding to myoglobin[END_REF] Therefore, we next monitored the conversion of rHbI-H 2 S to rHbI-CN over time by UV/Vis spectroscopy at 420 nm (Figure 2 right). It follows pseudo-first order kinetics, and the dissociation rate constant for hydrogen sulfide was calculated with equation (2) (Figure S5), yielding a k off = (7.5 ± 2.5) x 10 -3 s -1 . This new value explains nicely the experiment described in Figure 2 (left). It also rationalized the previously high K' d values determined by titration experiments with rHbI and wtHbI (332 and 90 nM, respectively), which are now in the range of those calculated by the ratio k off / k' on (131 nM and 115 nM, respectively). Using HbI to detect H 2 S by fluorescence. With the idea to use the "absorption into fluorescence" strategy for developing a turnon fluorescent sensor to detect hydrogen sufide, we selected N-[2(2-Hydroxyethoxy)ethyl]-6,8-difluoro-7-hydroxycoumarin-3carboxamide, a derivative [START_REF] Kerkovius | A Practical Synthesis of 6,8-Difluoro-7-hydroxycoumarin Derivatives for Fluorescence Applications[END_REF] of the coumarin dye "Pacific Blue TM " (PB ‡ ), as fluorophore. We anticipated that the intense absorption of the Soret band of rHbI at 407 nm would efficiently filter the light at the maximum excitation wavelength of the coumarin (404 nm), and turn its fluorescence off in the absence of the analyte hydrogen sulfide. Because of the red-shift (18 nm) and the decrease in the attenuation coefficient associated with H 2 S binding to rHbI (178000 vs 102000 L.mol -1 .cm -1 , respectively), [START_REF] Boffi | Determination of H2S solubility via the reaction with ferric hemoglobin I from the bivalve mollusc Lucina pectinata[END_REF] the fluorescence should be partially restored in the presence of the sulfur-containing analyte, as schematized in Figure 3. To confirm our aforementioned hypothesis, we first added increasing amounts of rHbI to a buffered solution of the fluorophore (Figure S6). As expected, this leads to a decrease in the fluorophore's fluorescence intensity, both in the emission and excitation spectra (Figure S6). For instance, a 35-fold decrease in the excitation and emission spectra of a 1 µM solution of the fluorophore is observed in the presence of 20 µM of rHbI. Next, we added sodium hydrosulfide to the previous solution to decipher if the presence of the analyte partially restores the fluorescence of the fluorophore. As expected, the addition of NaSH immediately leads to a substantial increase in the fluorescence intensity (Figure 4). To substantiate the principle pictured in Figure 3, we first ruled out the direct reaction between PC and sodium hydrosulfide as a potential explanation for the gain of fluorescence intensity (Figure S7). We also confirmed that the red-shift of the Soret band of rHbI upon coordination of hydrogen sulfide occurs in the presence of PC (Figure S8). Finally, we observed that the substitution of rHbI by met-myoglobin (which possesses a Soret band at 408 nm but does not tightly binds hydrogen sulfide) results in a quenching of the fluorescence of PC, but not to its restoration in the presence of sodium hydrosulfide (Figure S9). All these data indicates that our system indeed operates by "absorption into fluorescence". To better define it, we performed experiments with a fixed concentration of the fluorophore PB (1 µM) and various rHbI and NaSH concentrations. For sodium hydrosulfide concentrations higher than 1 µM, we found that 1 µM of PB and 20 µM of rHbI allows a linear quantification of up to 25 µM sodium hydrosulfide (Figure 4 left). At saturation, the maximum recovery of fluorescence was a 17-fold increase in intensity. For sub-micromolar concentrations of sodium hydrosulfide (Figure 4 right), we observed that a lower concentration of the heme sensor (7 µM) was better suited. Under these conditions, we determined quantification limits (calculated as being 10 time the standard error / slope) of 226 and 294 nM for sodium hydrosulfide, using excitation and emission spectra, respectively. It must be noted that rHbI is specific for the protonated form of the analyte H 2 S, and that the aforementioned values are likely to be reduced at lower pH, as the affinity of rHbI for the analyte will increase (K d = [K' d /(1+(10 - pKa /10 -pH )], with pKa = 6.8 for H 2 S). Next we investigated the selectivity of our detection system. As shown in Figure 5 (left), our system is selective for hydrogen sulfide, since no significant variation of the fluorescence intensity is observed in the presence of physiological concentrations of glutathione (5 mM), cysteine (100 µM), or other possible interfer-ing species. Prominently, our sensor system is reversible, as depicted in Figure 5 (right). Indeed, after addition of a saturating concentration of sodium hydrosulfide to our gating system, the initial fluorescence ratio is recovered within minutes by carefully bubbling argon. In addition, the system can be cycled through several times. Last, to establish if our sensing system could permit H 2 S detection in more complex biological media, we tested it in human plasma. Although we were unable to detect a difference in the fluorescence intensity recorded with samples untreated or treated with lead acetate, to trap hydrogen sulfide as insoluble lead sulfide, [START_REF] Delwiche | A micromethod for the determination of hydrogen sulfide[END_REF] supplementing the plasma with increasing concentrations of sodium hydrosulfide results in an increase in the fluorescence intensity (Figure 6, left). Moreover, the sensor displays a good stability in the plasma (the fluorescence is unchanged after 10 minutes), does not exhibit any significant reactivity with high concentrations of glutathione (10 mM), and permits real time monitoring of H 2 S release by a glutathione-activated donor [START_REF] Artaud | A Persulfide Analogue of the Nitrosothiol SNAP: Formation, Characterization and Reactivity[END_REF] (Figure 6, right).

Figure 5. Left: selectivity of the fluorescent detection system towards potential interfering biomolecules at physiologically relevant concentrations (5 mM of glutathione, 1 mM for sodium thiosulfate and 100 µM for other substrates). Conditions: 1 µM of fluorophore, 20 µM HbI, in buffer A with 0.9 % DMSO, with (red) or without (black) 5 µM NaSH. Right: reversibility of the detection of hydrogen sulfide. Each cycle consists in i) the addition of 40 µM sodium hydrosulfide and 2) the bubbling of argon for 5 minutes. The reaction was carried out with 1 µM of fluorophore and 12 µM rHbI in buffer A with 0.9 % DMSO. The reaction mixture also contains catalase (500 U) that was found to stabilize the system during the bubbling phase.

DISCUSSION

To design an efficient molecular gating system, the choice of both the sensor and the transducer are crucial. Thus, we decided to work with the recombinant hemoglobin rHbI from the clam Lucina Pectinata as a sensor. Indeed, rHbI has been reported to be the most avid biological system for hydrogen sulfide and it should accordingly enable low detection limits. In addition, HbI exhibits a strong absorbance band in the Soret region, which is modulated in wavelength and intensity by the coordination of H2S. While the high affinity of rHbI for H2S has been attributed to the combination of a very fast association and of a slow dissociation step, we observed inconsistencies between the dissociation constants determined either under kinetic or thermodynamic conditions. Accordingly, we reinvestigated the kinetics for hydrogen sulfide reactivity with the ferric center of rHbI. Our experiments validated the fast coordination of H2S to rHbI, but they revealed that the dissociation of H2S from rHbI takes place at a faster rate than the one previously reported. This discrepancy most probably lies in the slow recording interval used by Kraus et al. [START_REF] Kraus | Hemoglobins of the Lucina-Pectinata Bacteria Symbiosis .1. Molecular-Properties, Kinetics and Equilibria of Reactions with Ligands[END_REF] (600 s) for their kinetic study. The new half-life (92 s) for the unbinding of H2S allows us to reconcile the kinetic and thermodynamic data, and it is also fully consistent with the reversibility of our sensing system towards H 2 S (Figure 5B) which takes place within minutes. For the transducer, we selected a derivative of the coumarin Pacific Blue TM (PB) based on its photophysical properties. Its combination with rHbI results in significant improvements of our sensing system over those also based on the "absorption into fluorescence" approach to detect gaseous molecules like nitric oxide [START_REF] Strianese | Absorption into fluorescence. A method to sense biologically relevant gas molecules[END_REF][START_REF] Rapson | Fluorescent nitric oxide detection using cobalt substituted myoglobin[END_REF] or dioxygen. [START_REF] Strianese | Absorption into fluorescence. A method to sense biologically relevant gas molecules[END_REF] First, by choosing a fluorophore with a maximum excitation wavelength corresponding to that of the unliganted sensor rHbI, our system acts as a turn-on sensor, a strategy usually preferred over turn-off systems. Second, a large variation of the fluorescence signal intensity is observed between the off-form of the sensor rHbI and the turned-on form rHbI-H 2 S. Indeed, a 17fold variation in the fluorescence intensity F/F 0 is recorded upon saturation with H 2 S in a system containing 20 µM of rHbI and 1 µM PB. For comparison, variations of 1.3-and 7.1-fold were respectively obtained for the detection of nitric oxide with the couples cytochrome c peroxidase / fluoresceine and cobaltmyoglobin / cyan fluorescent protein (CFP). Third, our gated system enables the fast and easy quantification of hydrogen sulfide concentrations down to two hundred nanomolar. While this detection limit is middling compared to chromatographic methods that permit the detection of as low as 1 nM H 2 S concentration, [START_REF] Shen | Measurement of H2S In Vivo and In Vitro by the Monobromobimane Method[END_REF] our system presents several advantages compared to the later. Due to the extremely fast coordination step of the analyte to rHbI, H 2 S detection occurs in real time, a feature so far only proposed by polarographic sensors. [START_REF] Xu | Electrochemical hydrogen sulfide biosensors[END_REF] In addition, in contrast to classical fluorescent probes which are irreversibly consumed upon reaction with hydrogen sulfide, the coordination of H 2 S to rHbI is reversible under our experimental conditions. This allows the reversible detection of H 2 S, since no redox reaction is taking place between the analyte and the metal centre. Last, since the active site of rHbI only accommodates small uncharged ligands, [START_REF] Pietri | Factors Controlling the Reactivity of Hydrogen Sulfide with Hemeproteins[END_REF][START_REF] Boubeta | Access and Binding of H2S to Hemeproteins: The Case of HbI of Lucina pectinata[END_REF] the fluorescent detection of H 2 S is also selective, and physiological concentrations of other biological thiols, such as glutathione or cysteine, do not interfere with our gated system. Prominently, we applied our fluorescent sensing system to the quantification of hydrogen sulfide in human plasma. Since the regain of interest for the biological role of H 2 S in the late 1990's, plasmatic concentrations ranging from sub-micromolar to hundreds of micromolar have been allegated in the literature, [START_REF] Whiteman | Hydrogen sulfide and the vasculature: a novel vasculoprotective entity and regulator of nitric oxide bioavailability ?[END_REF] albeit mostly because of erroneous analytical procedures (such as the use of the methylene blue test [START_REF] Olson | A Practical Look at the Chemistry and Biology of Hydrogen Sulfide[END_REF] ). On our side, we were unable to unambiguously detect H 2 S in unspiked plasma samples with our sensor system. Thus, we did not observe any significant differences between untreated samples and samples treated with lead acetate to remove hydrogen sulfide. However, an immediate increase in fluorescence was recorded upon spiking plasma samples with concentrations as low as 1 µM hydrogen sulfide, indicating that our sensor operates in plasma and that H 2 S concentrations in our samples are below the quantification limit. This was confirmed by the monitoring of H 2 S release by a gluthathione activated H 2 S donor These results are in agreement with a very recent study that uses up to date LC-MS/MS techniques to detect H 2 S in biological matrices. [START_REF] Tan | New method for quantification of gasotransmitter hydrogen sulfide in biological matrices by LC-MS/MS[END_REF]

CONCLUSION

Overall, our results confirm the potential of the "absorption into fluorescence" strategy to detect analytes of biological interest. Thus, the combined use of the most avid biosystem for hydrogen sulfide and a rationally selected fluorophore allowed us to design a promising new fluorescence turn-on sensor for hydrogen sulfide. Our gated system detects down to two hundred nanomolar concentration of H 2 S, and is selective towards other biological thiols. Significantly, the detection occurs in real time, and reversibly. Although we focused on the application of the sensor in a biological media, we could also envision its use in more complex systems, like micro-fluidic devices to analyze environmental samples. To fully take advantage of our system, the immobilization of the protein on a transparent surface would be a great asset, and we are currently working on this in our group.

Figure 1 .

 1 Figure 1. UV-visible titration of a 1.5 µM solution of rHbI in buffer A (50 mM KPi, pH 7.4, 1 mM dtpa) by sodium hydrosulfide.

Figure 3 .

 3 Figure 3. Principle of the turn-on fluorescent system proposed in our work: the sensor rHbI plays the role of a filter at the excitation wavelength of fluorophore in absence of H 2 S, but not in its presence.

Figure 4 .

 4 Figure 4. Fluorescence ratio and corresponding standard deviations observed upon addition of NaSH to a 1 µM solution of fluorophore containing either 20 µM (right) or 7 µM rHbI (left) in buffer A (50 mM KPi, pH 7.4, 1 mM dtpa) containing 0.9 % DMSO). Excitation wavelength: 410 nm / Emission wavelength: 455 nm. Data are represented as the mean ± S.D. (n=3).

Figure 6 .

 6 Figure 6. Left: Emission fluorescence ratio recorded after spiking samples of human plasma with sodium hydrosulfide (two samples from different donors were used, with n = 2 for each). Conditions: 1 µM of fluorophore and 12 µM of rHbI, in plasma (with 0.9 % DMSO final). Data are presented as the mean ± S.D. (n=4). Right: Continuous monitoring of the fluorescence intensity under various conditions in plasma (with 0.9 % DMSO final). 0-10 min: sensing system alone (recording every 2 minutes); 10-20 min: addition of 10 mM of GSH at t=10 min (recording every 2 minutes); 20 -50 min: addition of 25 µM of the donor at t= 20 min (recording every minute).

ABBREVIATIONS

  ABBREVIATIONS † K' d and k' refers to the pH-dependent equilibrium constant or rate constant, at pH 7.4, with the approximation that [NaSH] = [H 2 S]. H 2 S: hydrogen sulfide; rHbI: recombinant hemoglobin I from Lucina Pectinata; PB: N-[2(2-Hydroxyethoxy)ethyl]-6,8difluoro-7-hydroxy coumarin-3-carboxamide; dtpa: diethylenetriaminepentaacetic acid; NaSH: sodium hydrosulfide.

ACKNOWLEDGMENT

We thank Jacques Ballalou (Institut Pasteur, Paris, France) for the gift of the BLi5 strain. We thank Maryse Jaouen (UMR 8601) for practical support during the protein expression, and Nassima Bekkadour (UMR 8601) for the preparation of the plasma samples. We thank the Laboratoire d'Electrochimie Moléculaire (UMR CNRS -P7 7591, Paris, France) for the use of the Stopped-Flow spectrometer. We finally thank Pr Marie-Agnès Sari and Drs Dominique Padovani and Diana Over (UMR 8601) for useful discussions and for critically reading and editing of the manuscript.

ASSOCIATED CONTENT

Supporting Information. Additional Figures mentioned in the main text. This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION

Corresponding Author * erwan.galardon@parisdescartes.fr.

Author Contributions

The manuscript was written through contributions of all authors and all authors have given approval to the final version of the manuscript.

ACS Paragon Plus Environment


Piotr Sierant

Marco Schirò

Maciej Lewenstein

Xhek Turkeshi

Entanglement Growth and Minimal Membranes in (d + 1) Random Unitary Circuits

Understanding the nature of entanglement growth in many-body systems is one of the fundamental questions in quantum physics. Here, we study this problem by characterizing the entanglement fluctuations and distribution of (d+1) qubit lattice evolved under a random unitary circuit. Focusing on Clifford gates, we perform extensive numerical simulations of random circuits in 1 ≤ d ≤ 4 dimensions. Our findings demonstrate that properties of growth of bipartite entanglement entropy are characterized by the roughening exponents of a d-dimensional membrane in a (d + 1) elastic medium.

Characterizing entanglement evolution in many-body quantum systems is an essential problem in statistical mechanics, condensed matter theory, and high-energy physics [START_REF] Laflorencie | [END_REF][2][START_REF] Potter | Quantum Science and Technology[END_REF]. For instance, it has profound implications for quantum thermalization [START_REF] Deutsch | [END_REF][5][6][7][8][9][10][11] and its breakdown in disordered systems [12][13][14][15][16][17][18][19][20]. At the same time, it provides new perspectives on quantum chaos arising in semiclassical and long-range systems [21][22][23][24][25] or holographic models [26][27][28][29][30][31][32][33]. Yet, understanding entanglement dynamics is of outstanding difficulty due to the nonlocal nature of the quantum correlations. Insights in one-dimensional systems come from conformal field theory and integrability [34][35][36][37][38][39][40][START_REF] Alba | Proc. Natl. Acad. Sci[END_REF][START_REF] Alba | [END_REF][43][44], exactly solvable dualunitary dynamics [45][46][47][48][49][50][51], and extensive numerical investigations [52][53][54][55][56][57][58][59]. These results motivate the search for a phenomenological understanding of quantum information scrambling in generic many-body systems.

Random unitary circuits are ideal candidates for this purpose [60]. In suitable limits, Ref. [61][62][63] demonstrated that the evolution of entanglement entropy of (1+1)D circuits maps to the properties of a line defect in a two-dimensional elastic manifold. We note that those insights go beyond the leading order contribution to the entanglement entropy growth and are reflected in structural aspects of entanglement propagation. A non-trivial verification of this fact is the Kardar-Parisi-Zhang (KPZ) universal growth of entanglement fluctuations in (1+1)D Clifford and Haar random unitary circuits [61][62][63]. Inspired by these results, the authors of Ref. [61] conjectured that generic (d+1) dimensional systems correspond to a d-dimensional membrane in a (d+1) elastic medium. This defect is fixed by the minimal cut of entanglement bonds across space-time, resulting, for random unitary circuits, in the entanglement growth being governed by exponents of a membrane pinned by the disorder [64][65][66][67]. Without any numerical and analytical results for d > 1, it remains unclear whether the minimal membrane picture indeed provides insights into the structure of entanglement evolution in higher dimensional systems.

This work positively answers this question by investi- The dashed lines denote a power-law increase t β where β is an extrapolated exponent. To perform the extrapolation, we fit w(t) for a given L with an algebraic dependence t β L in time interval [t1, t2], chosen in such a way that the value of βL is maximized and t2 -t1 > 40. The values of βL are shown as functions of 1/L in (d) for d = 2, 3, 4, the extrapolation to the L → ∞ limit is performed with a first order polymomial in 1/L, yielding β = 0.361 (6), β = 0.401( 9), β = 0.500( 5) respectively for d = 2, 3, 4.

circuit starting from a weakly entangled state |Ψ 0 ⟩ on the lattice Λ. Each layer Û = ⟨i,j⟩∈I

(2) t Ûi,j i∈I (1) t
Ûi is composed out of random two-body gates Ûi,j acting on neighboring sites i, j and on-site random gates Ûi applied to random sites belonging to the sets of indices

I (1,2) t
. We quantify the entanglement via the von Neumann entropy, defined for a bipartition A ∪ B of the lattice Λ and the state

|Ψ t ⟩ = Ût |Ψ 0 ⟩ at time (circuit depth) t as S(A) ≡ -tr(ρ A log 2 ρ A ) with ρ A = tr B (|Ψ t ⟩⟨Ψ t |).
The minimal membrane picture conjectured in Ref. [61] states that, at a coarse-graining level and for the region A bounded by the (d -1) dimensional surface ∂A, the entanglement entropy at time τ is given by

S mc (A, τ ) = min Σ S(A ′ , 0) + Σ dtd d-1 xE(v) . (1)
The minimization is over the d-dimensional membrane Σ(x, t) with boundary conditions Σ(x, τ ) = ∂A and Σ(x, 0) = ∂A ′ , and involves the region A ′ describing the entanglement of the initial state and a local velocity v defined using the tangent vectors with the minimal angle with respect to the time axis. We note that the dependence of the energy density on the local velocity in E(v) is a strong assumption in light of the nonlocal nature of S(A).

Without fluctuations, by dimensional analysis, the minimal membrane hypersurface in Eq. ( 1) scales as |∂A|τ , compatibly with the leading order average entanglement in Ref. [11]. The randomness of the circuit is reflected in the disordered nature of Σ(x, t). Assuming uncorrelated and Gaussian noise, Σ(x, t) is fixed as the ground state of a classical Hamiltonian [64,65,70]

H[Σ(x, t)] = dtd d-1 x κ 2 (∇Σ(x, t)) 2 + V (x) , (2)
with V (x) a Gaussian noise with zero mean V (x) = 0 and variance V (x)V (x ′ ) = γδ(x -x ′ ), for some parameters κ and γ. This problem is fully characterized by the exponents θ and ζ, related by the hyperscaling relation 2ζ = θ + d -2 [67,68] and governing the energy and width fluctuations, respectively The minimal membrane picture conjecture, cf. Eq. ( 1), provides a link between the fluctuations of the membrane Σ(x, t) and of the entanglement entropy. In particular, entanglement entropy for a fixed realization is given by [61]

E(x, t) ≡ ⟨H[Σ(x, t)] 2 ⟩ -⟨H[Σ(x, t)]⟩ 2 ∝ L θ , W (x, t) ≡ ⟨Σ(x, t) 2 ⟩ -⟨Σ(x, t)⟩ 2 ∝ L ζ . (3)
S(t) = ṽt|∂A|t + bt θ+1-d χ(t), (4)
where χ(t) is a stochastic variable, ṽ and b are two constants depending on the microscopic details of the model, and we neglected sub-subleading contributions. For concreteness, from now on, we consider Λ to be a hyperrectanglular

L 1 × • • • × L d lattice, with L 1 = L and L 2 = . . . = L d = L/2
, assuming periodic boundary conditions in all directions. As the subsystem A we take a hyperrectangular sublattice of

Λ of size x × L 2 × . . . × L d .
Then, the membrane growth exponents imply that the spatiotemporal fluctuations obey

w(x, t) ≡ ⟨S(x, t) 2 ⟩ -⟨S(x, t)⟩ 2 ∝ L (d-1)/2 t β , G(r, t) ≡ ⟨(S(x, t) -S(x + r, t)) 2 ⟩ = r β/ζ g r t ζ , (5)
where β = θ +(1-d)/2 is the growth exponent, ⟨•⟩ is the disorder average and g is a scaling function. In the following, we numerically extract the growth exponent β (and hence θ) and ζ by studying the entanglement entropy fluctuations in random circuits and comparing them with the results for disordered elastic manifolds [67]. Numerical results.-Assuming the membrane picture conjecture, the above discussion is relevant for any choice of the unitary ensemble of the random gates. Nevertheless, numerical simulations of generic (Haar) random circuits requires computational resources that scale exponentially with the number of investigated qubits. Therefore, this work resorts to Clifford circuits whose dynamics for stabilizer states on N qubits are efficiently simulable via the Gottesmann-Knill theorem [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF][START_REF] Gottesman | [END_REF][77][78].

Denoting i |m µ i ) with elements in 0, 1 and whose rows define the stabilizing Pauli strings. Clifford gates are defined as unitary operators that map a Pauli string into a single Pauli string. Hence, the action of Clifford gates can expressed as operations in the Z 2 field on the matrix M . It follows that the unitary evolution under Clifford circuits is simulable in resources scaling polynomially with the number of qubits N . Furthermore, the entanglement entropy of stabilizer states is efficiently computable as S(A) = r A -N A , where r A = rk 2 M A is the rank of the submatrix M A = (n µ i |m µ i ) i∈A and N A is the number of qubits in A [61]. We note that earlier simulations in the literature were based on standard Jordan-Gauss elimination, with complexity O(N 3 ). Here, to calculate entanglement entropy we employ an asymptotically fast rank-revealing algorithm with complexity O(N 3 / log 2 N ) [81,82]. We combine an efficient implementation of this algorithm [83] with a state-of-theart STIM [84] library to perform large-scale computations for systems of up to N = 131072 qubits. This allows us to calculate entanglement entropy for the largest considered N around two-orders of magnitude faster than a naive implementation of Jordan-Gauss elimination in STIM.

In our implementation, I

t is a set of N/2 sites at which the single qubit unitaries act. The one body gates are selected randomly, with equal probability, to be either Hadamard, H, or phase, √ Z, gates. Instead, I

(2) t is a set of N/4 pairs of nearest-neighboring sites ⟨i, j⟩ where i is a site picked randomly with uniform distribution on Λ and j = i + e u where e u is a versor in a randomly chosen direction u = 1, . . . , d. The two-body gates acting on I

(2) t are chosen, with equal probability, between CNOT gates controlled on the left or right argument [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF]. For concreteness, we consider the fully polarized initial state |Ψ 0 ⟩ = ⊗ i∈Λ |0⟩. We have tested other low-entangled ini-tial states finding quantitatively consistent results.

First, we consider the fluctuations w(t) ≡ w(L/2, t) of the entanglement entropy, cf. Eq. ( 5). Our findings are summarized in Fig. 2, where we let the system evolve up to time t = 4L. To extrapolate the growth exponent β L for a fixed system size L, we fit the numerical results with an algebraic dependence t β L in the time interval [t 1 , t 2 ]. We consider t 2 -t 1 > 40, chosen in order to maximize β L . We extrapolate the thermodynamic limit L → ∞ using a first-order polynomial in 1/L, obtaining respectively β = 0.361 (6), β = 0.401 (9) and β = 0.500(5) for d = 2, 3, 4 in agreement with the exponents obtained for the pinned membrane, c.f. Table 1, where the hyperscaling relation is used to translate the obtained values of β to the exponents θ.

Furthermore, we study the spatiotemporal fluctuations G(r, t), cf. Eq. ( 5). We fix the system sizes to be equal to L = 360 and L = 64 for d = 2 and d = 3, respectively, and study the growth of the fluctuations varying r and t. (The d = 2 case is presented in Fig. 3(a), while the result for d = 3 is qualitatively similar and are not shown here.) The spatiotemporal dependence of G(r, t) can be analyzed with the hypothesis (5), where the growth exponents β and ζ obtained by rescaling of results for different circuit depths t with the aim of finding an optimal collapse the data onto a single curve. To perform the analysis in an unbiased fashion, we consider both β and ζ as unknown parameters, and optimize the cost function

C(β, ζ) = n-1 i=2 c(x i , y i , d i ) with x i = (r i /t ζ i ), y i = G(r i , t i )/r β/ζ i
and d i = σ(y i ) respectively the hyperparameters, data, and standard deviations. The indices i = 0, . . . , N data are ordered such that x i < x j when i < j [85,86], and the density c(x i , y i ,

d i ) = (y -y) 2 /∆ depends on y = [(x i+1 -x i )y i-1 -(x i-1 -x i )y i+1 ] /(x i+1 -x i-1 )
and on

∆ = d 2 i + x i+1 -x i x i+1 -x i-1 d i-1 2 + x i-1 -x i x i+1 -x i-1 d i+1 2 .
The optimal growth exponents are obtained via minimization β ⋆ , ζ ⋆ = arg min C(β, ζ), and lead to the data collapses in Fig. 3 Still, these values are compatible with the membrane predictions, cf. Table 1.

Lastly, we study the distribution of the variable S = (S -⟨S⟩)/σ(S), with σ(S) the standard deviation and ⟨S⟩ the average. In one dimension, this variable follows the centralized Tracy-Widom (TW) distribution [61], whereas, for d = d c = 4, the system is expected to reach a Gaussian distribution [67,70]. For d = 2, 3 the distribution is unknown but is expected to differ from the TW and Gaussian ones. Our results for d = 2, 4 are given in Fig. 3(d) at t = L for L = 320 and L = 32, respectively. We see that both distribution are closer to the Gaussian fit than to the TW. To quantify the closeness to Gaussianity, we investigate the skewness κ 3 of P ( S). The TW distribution relevant for the d = 1 case yields κ 3 ≈ 0.2241, while we find that the d = 2 distribution has κ 3 = 0.07 [START_REF] Laflorencie | [END_REF], showing that the distribution is nongaussian. For d = 3, we obtain that κ 3 = 0.004(3). Hence, while the distribution P ( S) approaches the Gaussian distribution with increasing system dimension, we still observe mild deviations from gaussianity at the level of the P ( S) for d = 3, consistently with the non-mean field values of the exponents θ and ζ. Instead, for d = 4, the numerically obtained P ( S) that is indistinguishable from Gaussian, and we find κ 3 = 0.000(4) confirming the membrane picture prediction for the upper critical dimension.

Discussion and conclusion.-This manuscript investigates the structural properties of entanglement propagation in (d + 1) dimensional random unitary circuits, as encoded in the fluctuations and distribution of the entanglement entropy. Focusing on Clifford circuit implementations, we compute the growth exponents θ and ζ for d ≤ 4 with high accuracy. We find they are compatible with the corresponding exponents governing the fluctuations of a membrane pinned by disorder in an elastic manifold, cf. Table 1. Furthermore, at d = 4, the entanglement fluctuations become Gaussian, in accordance with the upper critical dimension of membranes pinned by disorder in d + 1 manifolds. Overall, our findings demonstrate the effectiveness of the recently proposed minimal cut (membrane) conjecture in capturing the es-sential structure of entanglement dynamics [61].

Beyond numerical considerations, suitable limits (e.g., large local Hilbert space dimension) may be analytically treatable for Haar and Clifford circuits even in higher dimensions [62,63,87]. At the same time, it would be interesting to investigate the properties of the volumelaw phase in (d + 1) monitored systems [85,86, and long-range systems [109][110][111][112][113]. While recent works for (1+1)D circuits revealed a rich structure related to polymer defects [114][115][116][117], the case of higher dimensions has been left essentially unexplored. A first step in (d + 1) dimension was recently presented by Ref. [118] for interacting Majorana circuits, where nontrivial subleading corrections due to the topological nature of the system were found. An interesting idea is to relate the L ln L subleading behavior to Fermi, surface corrections in the spirit of the Widom conjecture [119], or to the statistical mechanics of loop models [120][121][122][123]. We leave these questions for future investigation.

Figure 1 .

 1 Figure1. Cartoon: the entanglement growth of (d+1) dimensional random unitary circuits is conjectured to behave as a ddimensional disordered membrane in a (d+1) elastic manifold. Table:summary of the manuscript results (Entanglement) and comparison with the best numerical estimates[68,69] for random bond Ising model (Membrane), and with the renormalization group analysis[70][71][72][73][74] (RG). Our results on entanglement propagation are compatible with predictions for membranes pinned by disorder in (d + 1) elastic manifolds.

Figure 2 .

 2 Figure 2. Dynamics of entanglement entropy fluctuations in (d+1) dimensional random Clifford circuits. The standard deviation of entanglement entropy w(t) rescaled by L (d-1)/2 for system size L is shown, as a function of time t, in panel (a), (b), (c) respectively for d = 2, 3, 4.The dashed lines denote a power-law increase t β where β is an extrapolated exponent. To perform the extrapolation, we fit w(t) for a given L with an algebraic dependence t β L in time interval [t1, t2], chosen in such a way that the value of βL is maximized and t2 -t1 > 40. The values of βL are shown as functions of 1/L in (d) for d = 2, 3, 4, the extrapolation to the L → ∞ limit is performed with a first order polymomial in 1/L, yielding β = 0.361(6), β = 0.401(9), β = 0.500(5) respectively for d = 2, 3, 4.

  Renormalization group and numerical simulations predict ζ ≃ 0.21(4 -d) and θ = 0.84(3) (θ = 1.45(4)) in d = 2 (d = 3), while the fluctuations become Gaussian at the upper critical dimension d c = 4 [70].

  the Pauli strings as P = e iϕ i∈Λ Xn i i Ẑm i i , any pure stabilizer state |Ψ⟩ is determined by N = |Λ| independent Pauli strings Pµ as |Ψ⟩⟨Ψ| = N µ=1 ( 1 1 + Pµ )/2. Hence, the stabilizer state is encoded in a N × (2N + 1) matrix M = (ϕ µ |n µ

  (b) and (c) for d = 2 and d = 3 respectively. In the insets, we present the cost function landscape C(β, ζ). In d = 2, our estimates β = 0.351(12) and ζ = 0.435(33) correspond to a relatively narrow dip in the cost function leading to exponents that are compatible with the pinned membrane ones. In d = 3, the cost function, cf. Fig. 3(c,inset), possesses a wide minimum which leads to larger error bars on the estimates β = 0.39(3) and ζ = 0.3(1).

Figure 3 .

 3 Figure 3. Spatiotemporal fluctuations of entanglement entropy in random Clifford circuits. The entanglement entropy correlation G(r, t) is plotted in (a) as a function of r for various times t for the (2+1) dimensional circuit of size L = 360. The data for different times t are collapsed according to (5) in panel (b), the optimal obtained values of the exponents are β = 0.351(12), ζ = 0.435(33), the cost function w is shown in the inset. The same is shown for d = 3 in panel (c), the optimal exponents are β = 0.39(3), ζ = 0.3(1). Distribution of S -the rescaled entanglement entropy, at time t = L is shown in (d) for (2 + 1) and (4 + 1) dimensional circuits and compared with Gaussian and standardized Tracy-Widom distributions (dashed lines).

Acknowledgments.-We thank A. Nahum and T. Zhou for discussions. We acknowledge the workshop "Dynamical Foundation of Many-Body Quantum Chaos" at Institute Pascal (Orsay, France) for hosting P.S. and X.T. during the finalizing stage of this manuscript's writing. X.T. and M. S. acknowledge support from the ANR grant "NonEQuMat" (ANR-19-CE47-0001) and computational resources on the Collǵe de France IPH cluster. P.S. and M. L. acknowledge support from: ERC AdG NOQIA; Ministerio de Ciencia y Innovation Agencia Estatal de Investigaciones (PGC2018-097027-B-I00/10.13039/501100011033, CEX2019-000910-S/10.13039/501100011033, Plan National FIDEUA PID2019-106901GB-I00, FPI, QUANTERA MAQS PCI2019-111828-2, QUAN-TERA DYNAMITE PCI2022-132919, Proyectos de


Patrick C Chaumet

Guillaume Maire

Anne Sentenac

Accelerating the discrete dipole approximation by initializing with a scalar solution and using a circulant preconditioning

Keywords: Scattering, Discrete dipole approximation, Iterative method, preconditioner

The discrete dipole approximation (DDA) is a method of choice for simulating the electromagnetic scattering by objects of arbitrary shape and permittivity. To recover the field inside the object, it requires the iterative solving of a dense linear system which can be time consuming. To ease this task, we propose to start the inversion with the solution of the recently introduced scalar approximation [Chaumet et al. J. Opt. Soc . Am. A, 39, 1462 (2022)]. This initial guess allows a reduction of the time required for the solving of the linear system up to 50%. In addition, we study the interest of preconditioning the system to accelerate convergence. We show that the gain can be up to a factor of 5, especially for homogeneous objects on a plane substrate.

Introduction

Simulating the interaction of an electromagnetic field with a three dimensional object of arbitrary shape and relative permittivity is mandatory in many Preprint submitted to Journal of L A T E X Templates October 2, 2023 areas of photonics and surface science. There exist various techniques for solving Maxwell's equations in a complex environment in the time-harmonic regime [START_REF] Kahnert | Numerical methods in electromagnetic scattering theory[END_REF][START_REF] Gallinet | Numerical methods for nanophotonics: standard problems and future challenges[END_REF], including the finite element method (FEM), the multiple multipole method (MMP), the fast-multipole method (FMM), the method of moments (MoM) and the discrete dipole approximation (DDA).

In this article, we focus on the DDA, which is particularly adapted to the study of the light scattered by clusters of possibly anisotropic particles of arbitrary shapes [START_REF] Draine | The discrete-dipole approximation and its application to interstellar graphite grains[END_REF][START_REF] Chaumet | Coupled dipole method for scatterers with large permittivity[END_REF][START_REF] Moreno | The scattering matrix for size distributions of irregular particles: An application to an olivine sample[END_REF][START_REF] Yurkin | Accuracy of the discrete dipole approximation for simulation of optical properties of gold nanoparticles[END_REF]. Reviews on the potentiality of DDA can be found in [START_REF] Yurkin | The discrete dipole approximation: An overview and recent developments[END_REF][START_REF] Yurkin | The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength[END_REF], and many DDA codes are freely available on the net [START_REF] Draine | Discrete-dipole approximation for periodic targets: theory and tests[END_REF][START_REF] Yurkin | The discrete-dipole-approximation code adda: Capabilities and known lim itations[END_REF][START_REF] Chaumet | Ifdda, an easy-to-use code for simulating the field scattered by 3d inhomogeneous objects in a stratified medium: tutorial[END_REF].

The DDA can be described as a two step process. The first step, which is the main bottleneck of the approach, consists in estimating the field existing inside the object by solving a self-consistent equation discretized into a dense linear system. The second step consists in computing the scattered field at the required observation points. The second step has been solved recently with a Fourier based method which permits a significant reduction of the computation time [START_REF] Chaumet | Fast far-field calculation in the discrete dipole approximation[END_REF]. For the first step, combining the use of Fourier transforms (FFT) for the matrix vector products (MVP) [START_REF] Flatau | Light scattering by rectangular solids in the discrete-dipole approximation: a new algorithm exploiting the Block-Toeplitz structure[END_REF] and iterative solvers [START_REF] Da Cunha | The Parallel Iterative Methods (PIM) package for the solution of systems of linear equations on parallel computers[END_REF][START_REF] Flatau | Improvements in the discrete-dipole approximation method of computing scattering and absorption[END_REF][START_REF] Yurkin | The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength[END_REF][START_REF] Chaumet | Efficient iterative solution of the discrete dipole approximation for magneto-dielectric scatterers[END_REF], has permitted to accelerate this calculation. Yet, for large objects, the computation time of the MVP becomes important and any means for reducing the number of iterations of the solvers is welcome. To this aim, different implementations of the DDA were investigated [START_REF] Skorupski | Using the dda (discrete dipole approximation) method in determining the extinction cross sectionof black carbon[END_REF][START_REF] Yurkin | Performance of iterative solvers in the discrete dipole approximation[END_REF]. Thus, using the filtered coupled-dipole technique [START_REF] Piller | Increasing the performance of the coupled-dipole approximation: a spectral approach[END_REF][START_REF] Yurkin | Application of the discrete dipole approximation to very large refractive indices: Filtered coupled dipoles revived[END_REF] or the integration for the Green tensor [START_REF] Chaumet | Coupled dipole method for scatterers with large permittivity[END_REF][START_REF] Smunev | Rectangular dipoles in the discrete dipole approximation[END_REF] instead of the default point-dipole formulation was shown to reduce the number of iteration.

In this work, we will investigate two additional techniques for accelerating DDA simulations: starting the iterative process with an initial guess as close as possible to the actual solution and preconditioning the linear system.

Principle of the DDA

The discrete dipole approximation (DDA) was introduced by Purcell and Pennypacker in 1973 [START_REF] Purcell | Scattering and absorption of light by nonspherical dielectric grains[END_REF]. Since the DDA is a well known method [START_REF] Draine | The discrete-dipole approximation and its application to interstellar graphite grains[END_REF][START_REF] Flatau | Light scattering by rectangular solids in the discrete-dipole approximation: a new algorithm exploiting the Block-Toeplitz structure[END_REF][START_REF] Draine | Beyond Clausius-Mossotti: Wave Propagation on a Polarizable Point Lattice and the Discrete Dipole Approximation[END_REF][START_REF] Chaumet | Coupled dipole method for scatterers with large permittivity[END_REF][START_REF] Chaumet | Efficient iterative solution of the discrete dipole approximation for magneto-dielectric scatterers[END_REF], we shall outline only its main features. The object under study is represented by a cubic array of N polarizable subunits. The field at each subunit can be written:

E(r i ) = E ref (r i ) + N j=1 G(r i , r j )α(r j )E(r j ), (1)
where E ref (r i ) is the field at the position r i in the absence of the scattering object, G is the dyadic tensor associated to the free space [START_REF] Jackson | Classical Electrodynamics[END_REF] or a multilayer system [START_REF] Paulus | Accurate and efficient computation of the Green's tensor for stratified media[END_REF]. α(r j ) is the polarizability of the j th subunit. We use the polarizability of the Clausius-Mossotti relation with the radiative reaction term given by Draine [START_REF] Draine | The discrete-dipole approximation and its application to interstellar graphite grains[END_REF]:

α(r j ) = α 0 (r j ) 1 -(2/3)ik 3 0 α 0 (r j ) , ( 2
)
where α 0 holds the usual Clausius-Mossotti relation α 0 (r j ) = a 3 (ε(r j )-1)/(ε(r j )+

2). Equation ( 1) can be written symbolically as:

E = E ref + AD α E, (3)
where E, E ref are 3N vectors representing the unknown field and the reference field, respectively. A is a 3N × 3N matrix which contains the Green tensor, D α is a diagonal matrix of size 3N × 3N which contains the polarizability. The linear system is solved iteratively, which means that we define a residual as

r = ||E e -E ref -AD α E e || ||E ref || , ( 4
)
where E e is the field estimated iteratively and the iterative process is terminated once r < η where η is a prescribed tolerance. This tolerance is the parameter that quantifies the accuracy of the final field estimation. The iterative technique that is used in this article to estimate E e is the generalized product bi-conjugate gradient solver [START_REF] Zhang | Gpbi-cg: Generalized product-type methods based on bi-cg for solving nonsymmetric linear systems[END_REF] which is an alternative version of the bi-conjugate gradient stabilized solver, see Ref. [START_REF] Chaumet | Efficient iterative solution of the discrete dipole approximation for magneto-dielectric scatterers[END_REF][START_REF] Skorupski | Using the dda (discrete dipole approximation) method in determining the extinction cross sectionof black carbon[END_REF][START_REF] Yurkin | Performance of iterative solvers in the discrete dipole approximation[END_REF] for the efficiency of different iterative methods. In the iterative process numerous MVP are required. The matrix A being Toeplitz, placing the dipoles on a cubic lattice: i = (i x , i y , i z ), where

i x = 3 0, • • • , n x -1, i y = 0, • • • , n y -1 and i z = 0, • • • , n z -1 with r i = (i x d, i y d, i z d)
and N = n x n y n z permits to calculate the MVP using Fast Fourier transform (FFT). Notice that objects with arbitrary shape can still be represented by a cubic lattice, by simply setting to zero the polarizations at the lattice sites lying outside the object boundaries [START_REF] Goodman | Application of fast-fourier-transform techniques to the discrete-dipole approximation[END_REF]. The main advantage of DDA is that only a small portion of space containing the object is meshed and there is no need for absorbing boundary conditions, thus it is very memory efficient. Its main problem lies in the convergence of the iterative solver which can be very slow, or worse, fail, when the object is large compared to the illumination wavelength and/or is highly contrasted.

In the following, we study the influence of the initial estimate and the preconditioner on the time required for the iterative solver to converge. All the simulations were performed with the freely accessible code IFDDA which can be found in Ref. [START_REF] Chaumet | Ifdda, an easy-to-use code for simulating the field scattered by 3d inhomogeneous objects in a stratified medium: tutorial[END_REF]. Other codes exist for DDA (see Ref. [START_REF] Chaumet | The discrete dipole approximation: A review[END_REF] appendix B for a list of Open-Source DDA), but they are all based on the same principle, so the results presented in this article are general. The Fast Fourier Transforms (FFT) are computed with the Fastest Fourier Transform in the West [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. The calculations are parallelized (OpenMP) on a computer with 12 Intel(R) Xeon(R) CPU E5-2687W v2 at 3.40GHz.

Influence of the initial guess of the iterative solver on its computation time

The iterative estimation of the field inside the sample usually starts with a null field, i.e. E ini = 0. Yet, taking an initial guess that is closer to the final solution, such as the field estimated by a fast approximate method, is a better choice for diminishing the number of iterations. Hence, if one does not depart too much from the validity domain of Born approximation , i.e. for objects with small contrast and small size compared to the wavelength, the initial guess can be the field that would exist without the sample,

E ini ≈ E ref . ( 5
)
For larger objects that mainly diffract in the forward direction, one can use the field given by the beam propagation method (BPM) [START_REF] Kamilov | Learning approach to optical tomography[END_REF][START_REF] Kamilov | Optical tomographic image reconstruction based on beam propagation and sparse regularization[END_REF][START_REF] Chaumet | Reflection and transmission by large inhomogeneous media. validity of b orn, rytov and beam propagation methods[END_REF] or the Rytov approximation [START_REF] Carminati | Phase properties of the optical near field[END_REF]. When the sample polarizability varies slowly compared to the wavelength, one can use the Classical Scalar Approximation (CSA) [START_REF] Born | Principles of Optics[END_REF].

This technique consists in approximating the Green's tensor appearing in the linear system Eq. ( 1) by the diagonal tensor

G = k 2 0 e ik 0 r
r I where I is the identity tensor [START_REF] Carney | Inverse scattering for near-field microscopy[END_REF][START_REF] Carney | Theory of total-internal-reflection tomography[END_REF][START_REF] Carney | Three-dimensional total-internal reflection microscopy[END_REF]. This approach yields a scalar linear system for each field component. If the incident field polarization is constant throughout the object (directed along û, with û * .û = 1), and the field inside the object is only weakly depolarized, another scalar approximation, named uGu, can be used. By assuming that the field inside the object remains polarized along û, E = eû, the uGu method transforms the vectorial linear system given by Eqs. ( 1) and ( 3) into a scalar system where the Green's tensor is replaced by the function û * .Gû and the unknown vector E by the unknown scalar component e [START_REF] Chaumet | Scalar approximation of maxwell equations: derivation and accuracy[END_REF]:

e(r i ) = E ref (r i ).û * + N j=1 û * .G(r i , r j )ûα(r j )e(r j ). (6)
It was shown in Ref. [START_REF] Chaumet | Scalar approximation of maxwell equations: derivation and accuracy[END_REF] that the uGu method was significantly more accurate than CSA for a similar computational cost.

In the following, we consider objects that are illuminated by a plane wave so that the incident field inside the sample is directed along û. We study the computation time for solving the linear system Eq. ( 3) when the initial guess is the field given by Born, Rytov, BPM, CSA or uGu method. Note that while the field given by the Born, Rytov and BPM techniques is obtained almost instantly, that given by the CSA and uGu requires the solving of a linear system of size N × N . The latter is performed with the same conjugate gradient technique as that used for solving the vectorial system, Eq. ( 3), except that the stopping criterion is different, η s for the scalar problem, η for the vectorial one. Obviously, the time required for calculating the initial guess is taken into account when estimating the global computation time.

The first object under study consists of a cube of relative permittivity ε = 1.6

and increasing side a, the mesh-size being set to d = λ/10. The iterative solving of the vectorial system Eq. ( 3) was conducted with a stopping criterion η = 10 -4 while that of the scalar CSA or uGu systems was η s = 0.01. We plot in Fig. 1(b) the computation time of the different methods divided by the time obtained when uGu provides the initial guess, as a function of the cube size, from a = 10λ to a = 20λ. It is seen that the uGu field is the best initial guess for reducing the computation time, especially for large objects. Note that decreasing the computation time for small objects, a ≤ λ, is of minor interest, as the calculation is very quick. For a = 20λ, using the uGu field as initial guess yields a time gain of at least 25% compared to all the other approaches.

To complete this analysis, we plot in Fig. 1(a) the number of MVPs required by the iterative solver to converge depending on the initial guess. As expected, starting with the uGu field permits a significant decrease in the number of MVP.

This result shows the interest of providing an accurate initial guess for starting the iterative inversion scheme, even though the latter also requires the solving of a linear system. At this point, it is worth studying the influence of the stopping criterion η s on the time gain. Obviously, large η s diminishes the calculation time of the scalar field but increases its inaccuracy while small η s improves the scalar field estimation but increases its computation time. Comparing the number of MVPs required for the solving of Eq. ( 3) when the uGu field is estimated with η s equal to 0.1, 0.01, 0.001 shows that the best compromise is η s = 0.01, see Fig.

In a second example, we consider the same object deposited on a glass substrate of refractive index 1.5 and illuminated under normal incidence from the glass side or under total internal reflection (classic microscopy configuration [START_REF] Chaumet | Superresolution of threedimensional optical imaging by use of evanescent waves[END_REF]). Few approximate methods are able to simulate the field inside an object above a substrate and we could only resort to the null, Born and uGu fields for providing the initial guess. Figures 3 and4 show the normalized computation time together with the number of MVPs required for solving the linear system with these different initial guesses. They confirm that, despite the computation cost of its calculation, the uGu field is the best option for initializing the iterative solver as it yields a time gain up to 20% under normal incidence and 50% under total internal reflection.

In an attempt to ameliorate further the initial guess, we implemented the Hagedüs trick which consists in optimizing a scalar so that the initial field times this scalar minimizes the error on the linear system [START_REF] Strakoš | On numerical stability in large scale linear algebraic computations[END_REF]. Unfortunately, except when taking the Rytov field as initial guess, this technique did not yield any improvement (and even though faster, the Rytov solution remained the slowest).

Figure 5 studies the influence of the relative permittivity ε for a cuboid of size (a×a×a/2) with a/λ = 20 on the computation time. Taking the uGu solution as initial guess permits a significant reduction of the computation time, up to 40%, as long as the cuboid permittivity remains smaller than 2. Indeed, cuboid with permittivities higher than 2 tend to depolarize significantly the incident beam and the field obtained with the uGu approximation is too far from the actual one [START_REF] Chaumet | Scalar approximation of maxwell equations: derivation and accuracy[END_REF] to be useful. Note that very recently, Inzhevatkin and Yurkin show that the approximation of WKB (Wentzel Kramers Brillouin) as an initial estimate could also reduce the number of iterations for objects with low permittivity (ε < 1.3) [START_REF] Inzhevatkin | Uniform-over-size approximation of the internal fields for scatterers with low refractive-index contrast[END_REF]. This confirms the fact that DDA is mainly suitable for optically soft particles, when the particles become larger than wavelength of illumination and that other techniques should be preferred when the permittivity exceeds 2.5 [START_REF] Liu | Comparison between the pseudo-spectral time domain method and the discrete dipole approximation for light scattering simulations[END_REF]. In this case, a possible solution for speeding up the DDA consists in preconditioning the matrix of interaction to change its spectrum [START_REF] Lv | A modified t. chan's preconditioner for toeplitz systems[END_REF]. This is what we will study in the next section.

Influence of the preconditioner on the computation time

In this section, we investigate the interest of preconditioning the linear system for ameliorating the convergence of the iterative solver. This study is conducted for both the vectorial and scalar linear systems.

Construction of the preconditioner

We first rewrite the linear system to be solved Eq. ( 3) in the form,

AE = E ref , (7)
where A = I -AD α . Instead of solving the original linear system above, one may consider the right preconditioned system,

AP -1 (PE) = E ref , ( 8
)
and solve AP -1 y = E ref with E = P -1 y. Alternatively, one may solve the left preconditioned system

P -1 AE = P -1 E ref . ( 9
)
Both systems give the same solution as the original system as long as the preconditioner matrix P is non singular. Recently, Groth et al. [START_REF] Groth | Accelerating the discrete dipole approximation via circulant preconditioning[END_REF] proposed a left preconditioning strategy based on the two level circulant preconditioner of Chan and Olkin [START_REF] Chan | Circulant preconditioners for toeplitz-block matrices[END_REF] for the DDA in the homogeneous configuration. The idea consists in finding a two level circulant matrix, the closest possible to A, which can be easily invertible. To this aim, they consider a configuration in which the object polarizability, α, is assumed to be constant over the cubic lattice.

For inhomogeneous samples, α is the average of the polarizability over the cubic lattice which corresponds to the Maxwell Garnett 'homogenized' medium.

When α is a constant, the matrix A can be written, in the general case of an object placed in a multilayered environment, as a block Toeplitz matrix (along i y ) of Toeplitz sub-matrices (along i x ). To obtain an easily invertible matrix,

A is approximated by a circulant matrix along these two dimensions. More precisely, the preconditioning matrix reads

P ix,iy,iz1,iz2 = δ(i -0) -α (n y -i y ) n y (n x -i x ) n x g ix,iy,iz1,iz2 -α (n y -i y ) n y i x n x g -(nx-ix),iy,iz1,iz2 - i y n y (n x -i x ) n x g ix,-(ny-iy),iz1,iz2 - i y n y i x n x g -(nx-ix),-(ny-iy),iz1,iz2 , ( 10
)
where g is one of the nine components of the dyadic Green's tensor and δ is the Dirac's delta function.

To obtain P -1 , each 2-level circulant-block matrix P(i z1 , i z2 ) of size n x n y is diagonalized using 2D-FFTs [START_REF] Chan | Fft-based preconditioners for toeplitz-block least squares problems[END_REF]. Then, the 9 components of the Green tensor are gathered to create a diagonal block matrix with n x n y blocks of size 3n z × 3n z . The inversion of each block is done thanks to a LU method. For the LU method we use the zgetrs routine from LAPACK [START_REF] Anderson | LAPACK Users' Guide, 3rd Edition[END_REF]. More details are given in Ref. [START_REF] Groth | Accelerating the discrete dipole approximation via circulant preconditioning[END_REF]. Once P and P -1 are calculated, one solves iteratively the left or right conditioned linear system.

Testing the preconditioner for an object in homogeneous space

To check the interest of the preconditioning, we first consider a homogeneous cuboid of side a = b = 10λ, relative permittivity ε = 2, with an increasing height l z from λ/5 to 10λ. The tolerance of the iterative method is fixed to η = 10 -4 and a null field is taken as initial guess. decomposing n z are small, the 3D FFTs that are used for solving the linear system in absence of preconditioner (we recall that, in homogeneous space,

A is block Toeplitz along the three directions (i x , i y , i z ) while P is circulant only along i x and i y ) are very fast and the preconditioning is useless. On the contrary, when the prime factors of n z are large, the 3D FFTs are less efficient and preconditioning the system is interesting. It allows a faster convergence of the iterative solver, as shown in Fig. 7(a) where the residue of the iterative solver is plotted versus the number of MVPs for NP, LP and RP for a cuboid with a = b = 10λ, l z = a/5 and relative permittivity 2.

In Fig. 7(b) we investigate the interest of the preconditioning when the permittivity of the cuboid of Figure 7(a) is increased. It is shown that preconditioning is interesting when the permittivity contrast is moderate (ε between 1.2 to 2.3). For weak permittivities, NP requires few MVPs and preconditioning is useless, for strong permittivities (above 2.5) the convergence of the preconditioned system appears less stable than the non-reconditioned one. We now consider an inhomogeneous sample corresponding to a cuboid of side a = b = 10λ and l z = a/5 and random relative permittivity with Gaussian probability density of mean ε bg and variance σ 2 , with a Gaussian correlation function [START_REF] Zhang | Electromagnetic wave imaging of targets buried in a cluttered medium u sing a hybrid inversion-dort method[END_REF], defined by ε(r), ε(r

′ ) = ε 2 bg + σ 2 exp -r-r ′ 2 l 2 c
. Figure 8 plots the number of MVP and normalized time versus σ for l c = λ and ε bg = 2.

We observe that the preconditioning reduces the computation time by a factor of two when the inhomogeneity is moderate but is clearly inappropriate when the inhomogeneity is strong. We tried to improve its performance by extending to the two-dimensional case the optimal circulant preconditioner presented in [START_REF] Remis | Circulant preconditioners for domain integral equations in electromagnetics, Electromagnetics in advanced applications (ICEAA)[END_REF] for one-dimensional problems. Unfortunately, this modification did not yield any amelioration.

Testing the preconditioner for the scalar approximation

Recently we have introduced a scalar approximation that yields accelerate estimates of the field calculations with good accuracy for permittivities typically less than two [START_REF] Chaumet | Scalar approximation of maxwell equations: derivation and accuracy[END_REF]. Under this approximation, the 2-level Chan preconditioner is easy to implement as the Green's function has only one component. It becomes a diagonal block matrix made of n x n y blocks of size n z × n z the inversion of which is fast. We study the time gain brought by the preconditioning on the same cuboid of increasing height as before. Figure 9 shows that whatever the height, the number of MVPs and the computation time is drastically reduced by the preconditioner. Another study conducted with the inhomogeneous cuboid revealed the same behavior as with the full vectorial problem : the preconditioner is efficient for moderate inhomogeneity.

Testing the preconditioner for an object in presence of a substrate

For an object in the presence of a multilayer, the A matrix in Eq. ( 3) is only 2D Toeplitz in the x and y direction. This means that, except when the Figure 9: Same as Fig. 6 but the vectorial linear system was replaced by the scalar linear system using the uGu approximation.

object is in the substrate or superstrate [START_REF] Yurkin | Rigorous and fast discrete dipole approximation for particles near a plane interface[END_REF], the matrix vectors product in the z-direction can not be done by FFT [START_REF] Schmehl | Discrete-dipole approximation for scattering by features on surfaces by means of a two-dimensional fast fourier transform technique[END_REF]. This configuration seems thus better adapted to the use of the two-level circulating preconditioner. To test its performance, we study the same configurations as previously but the cuboids are now deposited on a substrate of permittivity 2.25 and illuminated under normal incidence. We observe in Figs. 10(a) and 10(b) that, in this case, the preconditioner reduces the number of MVPs and allows a significant time gain whatever the height of the cuboid.

Similarly, when the permittivity of the cuboid is moderately increased, Fig. 11, the time gain brought by the preconditioner is important and can reach a factor of four. Yet, similarly to the homogeneous case, the preconditioner renders the iterative solver unstable for permittivity above ε = 2.5.

Conclusion

The new formulation of the scalar approximation that we have introduced in [START_REF] Chaumet | Scalar approximation of maxwell equations: derivation and accuracy[END_REF] is an excellent initial estimate for the iterative method that solves the DDA linear equation system. It reduces the computation time by 50% for moderately contrasted objects (permittivity below 2) such as those encountered in biological applications. We observed that this initial estimate was particularly efficient when the illumination was inhomogeneous, as shown for evanescent illumination.

In addition, we have studied the interest of the 2-level circulating left preconditioner introduced by Groth et al. [START_REF] Groth | Accelerating the discrete dipole approximation via circulant preconditioning[END_REF] for further diminishing the computation time. Groth et al. showed that this preconditioner was efficient for homogeneous flat object. We extended it to inhomogeneous objects and showed its general interest for speeding up the solving of Maxwell equation under the scalar approximation. In the vectorial case, its interest was more limited and concerned weakly inhomogeneous objects with moderate permittivity deposited on a substrate.

Figure 1 :Figure 2 :

 12 Figure 1: Cuboid of side (a × a × a/2) and relative permittivity ε = 1.6 with a meshsize d = λ/10. (a) Number of MVP versus a for the different initial guess. (b) Time of computation normalized to the time used by the scalar approximation uGu.

Figure 3 :

 3 Figure 3: Cuboid of side (a × a × a/2) and relative permittivity ε = 1.6 with a meshsize d = λ/10 deposited on a glass substrate. (a) Number of MVP versus a for the three different initial guess: Born, null and uGu. (b) Computation time using the Born or null field as initial guess divided by the time obtained when the scalar approximation uGu is the initial guess.

Figure 4 :

 4 Figure4: Same legend as in Fig.3, but the illumination is done in total internal reflection.

Figure 5 :

 5 Figure 5: Cuboid of side (a×a×a/2) with a/λ = 20 with a meshsize d = λ/10 in homogeneous space. (a) Number of MVP versus ε for the three different initial guess: Born, null and uGu. (b) Computation time using the Born or null field as initial guess divided by the time obtained when the scalar approximation uGu is the initial guess.

Figure 6 (Figure 6 :

 66 Figure 6: Solving the diffraction by a cuboid of width a = b = 10λ and increasing height, lz/a (a) Number of MVPs required by the iterative solver in log scale with no preconditioner (NP), left preconditioner (LP) and right preconditioner (RP). (b) Computation time required by the solving of the left or right preconditioned linear system divided by the time without preconditioner (NP).

Figure 7 :

 7 Figure 7: Cuboid of side a = b = 10λ and lz = a/5. (a) Evolution of the residue in log scale as a function of the number of MVPs of the iterative method for no preconditioning (NP), with a right (RP) and left preconditioner (LP) for ε = 2. (b) Number of MVPs in log scale versus the permittivity for η = 10 -4 .

Figure 8 :

 8 Figure 8: Inhomogeneous cuboid of size (10λ×10λ×2λ) with ε bg = 2 and lc = λ. (a) Number of MVPs versus σ for no preconditioning (NP) with a right (RP) and left preconditioner (LP). (b) Time of computation with left (LP) and right (RP) preconditioner divided by the computation time obtained with no preconditioner (NP).

Figure 10 :

 10 Figure10: Same as Fig.6but the cuboid is placed on a glass substrate located at z = 0 and the illumination is done along the z axis from the glass side.

Figure 11 :

 11 Figure 11: Cuboid of size (10λ × 10λ × 2λ) placed on a glass substrate located at z = 0 and the illumination is done along the z axis from the glass side. (a) Number of MVPs in log scale as a function of ε with no preconditioner (NP), left preconditioner (LP) and right preconditioner (RP). (b) Time of computation normalized with respect to the computation time without preconditioner (NP) for the right and left preconditioning.