Bénédicte M Lefèvre

Marine Delvigne

Josué Vidal

Virginie Courtier-Orgogozo

Michael Lang
email: michael.lang@ijm.fr

A novel mechanism for left-right asymmetry establishment involving tissue remodeling and MyoID

How novel left-right asymmetries evolve in animals is unknown. Drosophila pachea males display a unique left-right asymmetry in the genital organs and mate in a right-sided position onto the female. In flies, male genitalia undergo a clockwise 360° rotation during development. To test if this tissue remodeling process is linked to asymmetry development in the male genitalia of Drosophila pachea, we developed CRIPSR and transgenesis in this singular species. As in D. melanogaster, mutants of the MyoID gene homolog undergo reverse rotation of the developing genitalia. In addition, they have mirror-image asymmetric genitalia. Although their genital asymmetry is reversed mutants still adopt a right-sided copulation posture. Our results reveal a novel mechanism for establishing left-right asymmetry in animals, involving complex tissue remodeling and MyoID.

One Sentence Summary: D. pachea genital left-right size asymmetry is established through MyoID-mediated tissue remodelling and does not determine the direction of mating posture.

While the body of most bilaterian animals is superficially left-right symmetric, some internal organs display left-right asymmetries in position, shape or size (1). The generation of asymmetric organs during development starts with an initial symmetry-breaking step, which can involve several mechanisms: motile cilia and directional fluid flow, orientation of the cell division plane and distribution of molecules within cells (1). How over the course of evolution initially symmetrical organs become asymmetric remains elusive. Males of the fruitfly Drosophila pachea have an asymmetric phallus with a right-sided opening for sperm release (gonopore) and a pair of asymmetric external genital lobes, with the left lobe being 1.5 times longer than the right lobe (Fig. 1 A,B) [START_REF] Pitnick | New Species of Cactus-Breeding Drosophila (Diptera: Drosophilidae) in the Nannoptera Species Group[END_REF][START_REF] Lang | Distinct copulation positions in Drosophila pachea males with symmetric or asymmetric: External genitalia[END_REF][START_REF] Acurio | Repeated evolution of asymmetric genitalia and right-sided mating behavior in the Drosophila nannoptera species group[END_REF]. In addition, D. pachea males mate in a right-sided posture with the male being shifted about 6-8° towards the right side of the female [START_REF] Lang | Distinct copulation positions in Drosophila pachea males with symmetric or asymmetric: External genitalia[END_REF][START_REF] Acurio | Repeated evolution of asymmetric genitalia and right-sided mating behavior in the Drosophila nannoptera species group[END_REF][START_REF] Rhebergen | Drosophila pachea asymmetric lobes are part of a grasping device and stabilize one-sided mating[END_REF]. The lobes are considered an evolutionary novelty since they are not present in the three closely related sister species [START_REF] Lang | Radiation of the Drosophila nannoptera species group in Mexico[END_REF]. The asymmetric genitalia of D. pachea likely evolved recently within the past 3-6 1 Ma from ancestral symmetry [START_REF] Lang | Radiation of the Drosophila nannoptera species group in Mexico[END_REF] and this species is therefore a suitable organism to search for the evolutionary origins of genital left-right asymmetry.

During development, the male genitalia of a wide range of shizophoran flies undergoes a 360° clockwise rotation [START_REF] Suzanne | Coupling of Apoptosis and L/R Patterning Controls Stepwise Organ Looping[END_REF][START_REF] Feuerborn | Das Hypopygium 'inversum' und 'circumversum' der Dipteren[END_REF]. In D. melanogaster, the direction of this rotation is driven by the myosin MyoID, encoded by the gene myo31DF [START_REF] Spéder | Type ID unconventional myosin controls left-right asymmetry in Drosophila[END_REF], for simplicity here referred to as myoID. Loss of function of myoID in D. melanogaster leads to a counter-clockwise 360° genitalia rotation direction, due to the function of another counter-acting myosin, MyoIC [START_REF] Suzanne | Coupling of Apoptosis and L/R Patterning Controls Stepwise Organ Looping[END_REF][START_REF] Spéder | Type ID unconventional myosin controls left-right asymmetry in Drosophila[END_REF][START_REF] Petzoldt | DE-Cadherin regulates unconventional Myosin ID and Myosin IC in Drosophila left-right asymmetry establishment[END_REF]. MyoID also controls the chirality of other organs, such as the gut and testes [START_REF] Spéder | Type ID unconventional myosin controls left-right asymmetry in Drosophila[END_REF][START_REF] Taniguchi | Chirality in Planar Cell Shape Contributes to Left-Right Asymmetric Epithelial Morphogenesis[END_REF][START_REF] Benitez | Both JNK and apoptosis pathways regulate growth and terminalia rotation during Drosophila genital disc development[END_REF][START_REF] González-Morales | The Atypical Cadherin Dachsous Controls Left-Right Asymmetry in Drosophila[END_REF][START_REF] Hozumi | An unconventional myosin in Drosophila reverses the default handedness in visceral organs[END_REF][START_REF] Taniguchi | D-JNK signaling in visceral muscle cells controls the laterality of the Drosophila gut[END_REF][START_REF] Hatori | Left-right asymmetry is formed in individual cells by intrinsic cell chirality[END_REF][START_REF] Okumura | Class I Myosins Have Overlapping and Specialized Functions in Left-Right Asymmetric Development in Drosophila[END_REF]. Furthermore, MyoID appears to be a conserved key factor among different mechanisms to establish left-right asymmetry in vertebrates and invertebrates [START_REF] Juan | Myosin1D is an evolutionarily conserved regulator of animal left-right asymmetry[END_REF][START_REF] Tingler | A Conserved Role of the Unconventional Myosin 1d in Laterality Determination[END_REF], through interaction with adherens junctions [START_REF] Petzoldt | DE-Cadherin regulates unconventional Myosin ID and Myosin IC in Drosophila left-right asymmetry establishment[END_REF][START_REF] Taniguchi | Chirality in Planar Cell Shape Contributes to Left-Right Asymmetric Epithelial Morphogenesis[END_REF], the actin cytoskeleton [START_REF] Petzoldt | DE-Cadherin regulates unconventional Myosin ID and Myosin IC in Drosophila left-right asymmetry establishment[END_REF][START_REF] Lebreton | Molecular to organismal chirality is induced by the conserved myosin 1D[END_REF][START_REF] Chougule | The Drosophila actin nucleator DAAM is essential for left-right asymmetry[END_REF], planar cell polarity pathway members [START_REF] González-Morales | The Atypical Cadherin Dachsous Controls Left-Right Asymmetry in Drosophila[END_REF][START_REF] Juan | Myosin1D is an evolutionarily conserved regulator of animal left-right asymmetry[END_REF] and Jun kinase pathway regulated apoptosis [START_REF] Benitez | Both JNK and apoptosis pathways regulate growth and terminalia rotation during Drosophila genital disc development[END_REF][START_REF] Taniguchi | D-JNK signaling in visceral muscle cells controls the laterality of the Drosophila gut[END_REF][START_REF] Macias | PVF1/PVR signaling and apoptosis promotes the rotation and dorsal closure of the Drosophila male terminalia[END_REF][START_REF] Rousset | The Drosophila serine protease homologue Scarface regulates JNK signalling in a negative-feedback loop during epithelial morphogenesis[END_REF]. Ectopically expressed, MyoID and MyoIC are sufficient to induce de novo directional twisting at the cell, organ, or body scale in D. melanogaster [START_REF] Lebreton | Molecular to organismal chirality is induced by the conserved myosin 1D[END_REF]. However, the evolutionary transitions to novel leftright size asymmetries remain to be uncovered at the molecular, cellular and tissue scale.

We wondered if MyoID is involved in the establishment of left-right asymmetry in D. pachea male genitalia. To test this, we generated a mutation in the coding sequence of the D. pachea myoID using CRISPR/Cas9 mediated gene editing (24, 25) (myoID mut allele, Fig. S1). The induced mutation caused a frameshift in the seventh coding exon, leading to a premature stop codon. In homozygous myoID mut/mut males, the length ratio of left and right genital lobes was reversed with a longer right lobe and a shorter left lobe in most (81/82) dissected individuals (Data S1, Figure 1 C). Furthermore, the phallus asymmetry was reversed in myoID mut/mut males with the gonopore located at the left tip of the phallus (29/29 dissected genitalia), compared to heterozygous myoID wt/mut and myoID wt/wt males, where the gonopore was located consistently on the right side (8/8 and 13/13 dissected genitalia, respectively; Fig. S2). The locus myoID is autosomal (Fig. S1) and the mutation appears to be recessive since myoID wt/mut males did not have visible alterations of male genitalia morphology (Fig. 1 B). Apart from phallus and lobe morphology, homozygous mutant males revealed a genitalia rotation phenotype, with varying orientations of the male genitalia in adults, with respect to the antero-posterior midline (Fig. 1 A,D). In wild-type D. pachea, male claspers and genital lobes are pointing towards the ventral side (Fig. 1 A), while myoID mut/mut males revealed variable orientations of genitalia, which varied among individuals between 0 and 270° (Fig. 1 D). This latter orientation phenotype was similar to previous observations in D. melanogaster myoID mutants [START_REF] Spéder | Type ID unconventional myosin controls left-right asymmetry in Drosophila[END_REF]. Our data shows that MyoID is essential to determine left-right identity in D. pachea male genitalia, including external genital lobes and the male phallus. This function might be mediated through its role in determining the direction of male genitalia rotation. In one myoID mut/mut males, genital lobe asymmetry was not reversed although genitalia were mis-oriented and in one individual lobes were strongly bent. However, it was unclear if their genitalia rotated clockwise or counter-clockwise during development. These mutants may have had remnant activity and expression of the myoID mut allele, leading to a partial rotation in clockwise direction. Alternatively, other unidentified environmental or genetic factors may also contribute to the determination of genital asymmetry in D. pachea and were variable in our experiment.

To monitor genitalia rotation we inserted a DE-Cadherin::YFP fusion construct (DE-Cad::YFP) (Fig. S3) into the myoID mut/mut stock using piggybac transposon mediated germline transformation (26) (Data S2, Fig. S4). We examined male genitalia rotation in pupae that were heterozygous for DE-Cad::YFP to reduce potentially deleterious effects of the piggybac insert 2 (Fig. S5). We observed that myoID wt/wt and myoID wt/mut males underwent clockwise genitalia rotation direction (n = 9/9, 25/25, respectively), as observed in most schizophoran flies [START_REF] Suzanne | Coupling of Apoptosis and L/R Patterning Controls Stepwise Organ Looping[END_REF][START_REF] Feuerborn | Das Hypopygium 'inversum' und 'circumversum' der Dipteren[END_REF] whereas in myoID mut/mut males, genitalia rotated counter-clockwise (n = 9/9) (Fig. 1E, Data S1). In 9 monitored males that had developed at least to the pharate adult stage, we were able to dissect male genitalia and to determine the orientation and the lobe length ratio of male genitalia (Fig. S2). All these myoID mut/mut males with counter-clockwise genitalia rotation direction had a reversed lobe length ratio with a longer right lobe and a shorter left lobe, although in one observation lobes were only rudimatarily developed and the measurements were not informative about the lobe length ratio (Fig. 1E). In contrast, myoID wt/wt and myoID wt/mut males had a wild-type lobe ratio (Fig. 1E). We conclude that in D. pachea males MyoID controls genitalia rotation direction, such as previously observed in D. melanogaster and that in addition, MyoID determines left-right identity of genitalia.

Apart from its establishment during development, the role of genital asymmetry with respect to reproduction is not well understood [START_REF] Huber | Mating positions and the evolution of asymmetric insect genitalia[END_REF]. Genital asymmetry has been hypothesized to co-evolve or to evolve in response to changes in mating position [START_REF] Huber | Mating positions and the evolution of asymmetric insect genitalia[END_REF]. Alternatively, asymmetric genitalia may cause lateralized mating behavior since the asymmetric parts mediate the coupling of the female and male abdomen and therefore potentially direct the relative orientation of the female and male bodies. We investigated if males with reversed genital asymmetry would adopt a reverted copulation posture by monitoring copulation behavior (Fig. S6, Data S3) of myoID wt/wt , myoID wt/mut , or myoID mut/mut males with wildtype myoID wt/wt females. We annotated mating behavior (Data S4) and copulation postures (Data S5, Fig. S7). A copulation was considered to occur when the male mounted the female and achieved to settle into an invariant copulation position for at least 30 seconds. Alternatively, males were sometimes observed to mount the female, but leg and abdomen movements did not settle and eventually the male separated from the female within a few minutes, which we considered to be a failed mounting attempt (106/140 total attempts, Fig. S6). Only 26% (10/38) of the couples with myoID mut/mut males that mounted the female also adopted a stable copulation position, compared to 11/11 and 13/13 for myoID wt/wt and myoID wt/mut males, respectively (Fig. S6). Genitalia orientation angles of myoID mut/mut males that managed to adopt a stable copulation posture were found to deviate at most by 16.7° from wild-type orientations, while it was higher in the majority of myoID mut/mut males (18/29) that failed to copulate (Fig. 1D, Fig. S6). Such defects in genitalia orientation were previously reported to be detrimental for copulation and efficient genitalia coupling for D. melanogaster myoID mutants [START_REF] Inatomi | Proper direction of male genitalia is prerequisite for copulation in Drosophila, implying cooperative evolution between genitalia rotation and mating behavior[END_REF]. Regardless of the genotype and direction of genital asymmetry, myoID wt/wt , myoID wt/mut and myoID mut/mut males consistently adopted a right-sided average copulation posture (Fig. 2). We conclude that the right-sided mating posture in D. pachea does not depend on the direction of male genital asymmetry.

Overall, our study reveals functional conservation of the unconventional myosin MyoID to control clockwise genitalia rotation in developing D. pachea male genitalia. In addition, the same protein got co-opted during evolution of D. pachea to contribute to determination of the direction of a novel male genital asymmetry. The latter function may be mediated through the genitalia rotation process itself since rotation direction and genital asymmetry direction are perfectly associated in our live-imaging trials (Fig. 1E). Mechanical forces between rotating genitalia and non-rotating outer tissue may act in opposite direction at the left and the right side and in function of the sense of rotation. Male genital asymmetry in D. pachea does not determine the posture itself since wildtype males and also mutant males with reversed genital asymmetry consistently mate in a right-sided copulation posture. Right-sided mating is therefore female-controlled or may be a hard-wired behavior that could have possibly favored the evolution of morphological asymmetry as an instructive cue to optimize genital contacts [START_REF] Huber | Mating positions and the evolution of asymmetric insect genitalia[END_REF]. This study provides the first experimental insights into de novo Evolution of left-right asymmetry. It involves recruitment of a conserved left-right regulator and is associated with a pre-existing tissue remodeling process. 8 Fig. 2. The right-sided D. pachea mating posture does not depend on the direction of male genital asymmetry. Mating angles (α) (Supplementary Materials) of D. pachea copulating with a male being either (A) myoID wt/wt with the wild-typle lobe length ratio, (B) myoID wt/mut heterozygous males with wildtype lobe length ratio, (C) myoID mut/mut males with reverted lobe asymmetry ratio. Females had a myoID wt/wt genotype in all conditions. The hypothesis angle = 0 was rejected for each male genotype (GLM fit angle ~ genotype, all p < 1 x 10 -10 ). Data acquisition and processing is illustrated in Figs. S6,S7. Positive angle values correspond to right-sided orientations of the male head relative to the female midline. Data points that correspond to the same trial are connected by grey lines. The number of trials are indicated by n, a drawing of the external male genitalia of each genotype is shown on the right.

Fig. 1 .

 1 Fig. 1. MyoID mutants display reversed genital asymmetry and genitalia rotation defects. (A) Wild-type D. pachea male in ventral view. The orientation of genitalia is measured as an angle between the midpoint of male genital claspers relative to the male midline, as depicted by the dashed polar diagram. The scale is 100 µm. (B) External adult male genitalia of (B) a wildtype male myoID wt/wt with asymmetric genital lobes at the ventral side (arrows). (C) External adult male genitalia of a myoID mut/mut mutant male with reversed genital asymmetry. (B-C) Lobes are indicated by arrows and lobe lengths measurement landmark are shown as red dots. The scale is 100 µm. (D-E) Male genitalia orientation and lobe length asymmetry ratio in wild-type myoID wt/wt males (blue), heterozygous myoID wt/mut mutant males (grey) and homozygous myoID mut/mut mutant males (orange). (D) The polar diagram illustrates genitalia orientation angles. (E) The direction of genital lobe asymmetry is associated with the direction of male genitalia rotation. Left and right lobe lengths are plotted for individual males monitored by time-lapse microscopy. Genitalia rotation direction is indicated as +: clockwise, and as -: counter-clockwise. The diagonal red dashed line indicates the 1:1 lobe length ratio. The star indicates an observation with rudimentarily developed lobes in form of buds on both sides.

290

Acknowledgments:

We thank Stephane Prigent for preparation of Drosophila genitalia. We thank the CNRS ImagoSeine Imaging platform at the Institut Jacques Monod for assistance in microscopy. We are grateful to Nicolas Gompel for initial help with Drosophila germline transfection. Funding: BL was supported by a pre-doctoral fellowship Sorbonne Paris Cité of the Université de Paris and by a fellowship from the Labex ''Who am I?'' ["Initiatives d'excellence", Idex ANR-18-IDEX-0001, ANR-11-LABX-0071]. This work was further

Supplementary Materials:

Materials and Methods

Figures S1-S7

Table S1 Movies S1-S2 External Databases S1-S5 7


François Jean-Noël Lorenzi
email: jean-noel.lorenzi@ijm.fr

François Graner

Virginie Courtier-Orgogozo

Guillaume Achaz

CNCA aligns small annotated genomes

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A naive nucleotide alignment of complete genomes usually result in many frameshifts and other oddities that do not exist in protein alignments.

Several methods have been developed to perform nucleotide alignments taking protein alignment into account. One approach is "backtranslation", where coding nucleotide sequences are translated into amino acids, that are aligned. Corresponding codons are then aligned in a final nucleotide alignment. The web-based tool web-prank (https://www.ebi.ac.uk/goldman-srv/webprank/;

(Löytynoja and Goldman, 2010)) is such an example. Other tools based on backtranslation proposes specific options like the choice of genetic codes (PAL2NAL [START_REF] Suyama | PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments[END_REF]), transAlign (Bininda-Emonds, 2005), RevTans [START_REF] Wernersson | RevTrans: multiple alignment of coding DNA from aligned amino acid sequences[END_REF]). Some are designed to consider cases in which frameshifts or stop codons can occur (MACSE [START_REF] Ranwez | MACSE: Multiple Alignment of Coding SEquences Accounting for Frameshifts and Stop Codons[END_REF][START_REF] Ranwez | MACSE v2: Toolkit for the Alignment of Coding Sequences Accounting for Frameshifts and Stop Codons[END_REF], PAL2NAL [START_REF] Suyama | PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments[END_REF]), transAlign (Bininda-Emonds, 2005)). TranslatorX [START_REF] Abascal | TranslatorX: multiple alignment of nucleotide sequences guided by amino acid translations[END_REF] checks the relevance of the amino acid alignment by finding regions of uncertainties in the amino acid alignment (masked by Gblocks [START_REF] Castresana | Selection of Conserved Blocks from Multiple Alignments for Their Use in Phylogenetic Analysis[END_REF]) and reports them in the nucleotide alignment. Others are optimized for virus gene sequences (NucAmino [START_REF] Tzou | NucAmino: a nucleotide to amino acid alignment optimized for virus gene sequences[END_REF], VIRULIGN [START_REF] Libin | VIRULIGN: fast codon-correct alignment and annotation of viral genomes[END_REF]). To the best of our knowledge, none of these methods processes genome alignment with both coding and noncoding regions. We have thus developed CNCA (Coding / Non-Coding Aligner), a genome-wide solution that returns a complete genome alignment compatible with the protein sequence alignment. The method was designed for small (up to 50 kb) homologous annotated syntenic genomes devoid of introns, such as virus genomes. It will ease the subsequent evolutionary analysis of these complete genomes.

CNCA takes as input two or more GenBank files of annotated genomes (Fig. 1). It first MAFFT-aligns [START_REF] Nakamura | Parallelization of MAFFT for large-scale multiple sequence alignments[END_REF] the nucleotide (nt) sequences of all genomes and produces a Multiple Sequence Aligmnent (MSAnt). It then generates MSAaa, the MAFFT-alignment of the concatenate of all protein sequences. As the concatenate takes protein sequence on the order of gene annotations, synteny must be conserved. The MSAnt is then updated using MSAaa for all coding regions where both alignments are not concordant. A final MSAcnca is returned that contains no contradiction with MSAaa and thus no frameshift (figure 1). We choose to implement a graphical web version of the pipeline to widen the potential users to non-experts. Results (logs and the three alignments MSAcnca, MSAnt, MSAaa in both nexus and fasta formats) are stored locally for a week. An email with a link to access the results is sent to the user at the end of the procedure.

As an illustration, the whole CNCA pipeline runs in 45 min for a dataset of 12 genomes closely related to SARS-CoV-2.

Figure 1 .

 1 Figure 1. Schematic view of the CNCA pipeline.

Funding

This work was supported by the Labex "Who Am I?", ANR-11-LABX-0071 and the Université Paris Cité, Idex ANR-18-IDEX-0001, funded by the French Government through its "Investments for the Future" program.

The interface is available at https://cnca.ijm.fr/ and the code is available at https://github.com/jnlorenzi/

Conflict of Interest: none declared.


Elsa Giraudat

Arnaud Burtin

Arthur Le Ber

Mathias Fink

Jean-Christophe Komorowski

Alexandre Aubry
email: alexandre.aubry@espci.fr.

Unveiling the deep plumbing system of a volcano by a reflection matrix analysis of seismic noise

In geophysics, volcanoes are particularly difficult to image because of the multiscale heterogeneities of fluids and rocks that compose them and their complex non-linear dynamics. By exploiting seismic noise recorded by a sparse array of geophones, we are able to reveal the magmatic and hydrothermal plumbing system of La Soufrière volcano in Guadeloupe. Spatio-temporal crosscorrelation of seismic noise actually provides the impulse responses between virtual geophones located inside the volcano. The resulting reflection matrix can be exploited to numerically perform an auto-focus of seismic waves on any reflector of the underground. An unprecedented view on the volcano's inner structure is obtained at a half-wavelength resolution. This innovative observable provides fundamental information for the conceptual modeling and highresolution monitoring of volcanoes.

In everyday life, a multitude of sensors surround us to monitor our environment. In wave physics, those sensors can be active and work together to control the wave-field at will whether

it be for focusing (1) or communication (2) purposes. For imaging, the problem is often illposed because of the medium complexity and/or the sensor array sparsity. This is particularly the case in seismology, where the topography of the site under investigation can be so irregular that it is illusory to deploy a large and dense network of geophones. This paper addresses the issue of seismic imaging in complex areas such as volcanoes or fault zones based on data recorded by a sparse array of seismometers. The goal is to provide high spatial resolution and in-depth imaging of such critical areas that are of paramount importance for Earth sciences. To that aim, we will build on a matrix imaging approach imported from other fields than geophysics, such as medical ultrasonics (3,[START_REF] Lambert | Proc. Nat. Acad. Sci. USA[END_REF] and optical microscopy [START_REF] Kang | [END_REF]6) that were designed for scales ranging from a few centimeters for ultrasonic waves to a few hundreds of nanometers for light. In contrast with concurrent seismic methods such as full waveform inversion (7), the strength of matrix imaging lies in the fact that: (i) it does not rely on a sophisticated wave velocity model whose knowledge is often limited and uncertain in geophysics; (ii) it is robust with respect to data quality which is a frequent issue in seismology.

Matrix imaging relies on the array response matrix that contains the set of impulse responses between each seismometer. Although a geophone is purely passive, cross-correlation of seismic noise received at two stations is known to converge toward the Green's function between receiving stations (8,9), as if one of them had been used as source, thus paving the way to passive matrix imaging (10)(11)(12). As surface waves dominate ambient noise, most past studies on the topic aimed at extracting surface wave properties from ambient noise correlations (13).

However, they also contain the contribution of body waves reflected by deep structures (14) and fluid reservoirs (10).

As a proof-of-concept, we here exploit seismic noise recorded by a sparse geophone network deployed at the surface of the La Soufrière volcano of Guadeloupe [START_REF]Institut De Physique Du Globe De Paris (IPGP[END_REF][START_REF] Burtin | Dense seismic monitoring of La Soufrière de Guadeloupe hydrothermal system[END_REF]. The covariance matrix of this seismic noise provides the reflection matrix that contains all the available information on the underground reflectivity. A numerical focusing process, often referred to as redatuming in seismology [START_REF] Schuster | Seismic Interferometry[END_REF], can then be applied to provide a confocal image of the subsoil reflectivity (10). This image is directly proportional to the axial fluctuations of the acoustic impedance associated with length scales typically of the order or smaller than the wavelength. It is therefore an extremely relevant observable for highlighting the presence of fluid-rock interfaces. However, the quality of the confocal image is drastically degraded by: (i) the mismatch between the wave velocity model and its real distribution that gives rise to a foggy image; (ii) the sparsity and finite size of the geophone network that limit its resolution. The former problem is solved by compensating wave distortions from the Earth surface, thereby revealing an helical conduit in the upper part of the volcano. The sparsity issue is addressed by an iterative phase reversal process driven from the k-space [START_REF] Najar | [END_REF][START_REF] Bureau | [END_REF] that resolves the deep reflectors with a transverse resolution of the order of a half-wavelength, thereby breaking the free space diffraction-limit usually limited by the array aperture. The inner structure of the volcano is revealed up to a depth of 10 km. It shows sub-horizontal bodies linked by thinner sub-vertical structures that match the current state-of-the-art conceptual and data-derived view of transcrustal magmatic systems. Such a mush-based model applies to numerous volcanic systems and has indeed been theorized for La Soufrière of Guadeloupe (20)(21)(22).

Canonical Reflection Matrix

Figure 1A shows the virtual network of 76 geophones whose distribution has been dictated by the topography of the volcano. It spans over a lateral extension d || = 1300 m and a vertical range d z = 500 m. The impulse response R(i, j, t) between each pair of stations (i, j) is estimated by cross-correlation of ambient seismic noise (9) (section S1). The set of impulse responses is stored in a time-dependent response matrix R(t).

This canonical reflection matrix is powerful since it enables a post-processing projection of seismic data into different mathematical bases. The reflection matrix can be investigated into the plane wave basis (or k-space) or any plane in the real space that sits between the Earth surface (u) and the expected focal plane (ρ) at a given time-of-flight t. To project the seismic data in these latter bases, a wave velocity model is nevertheless required.

Wave Velocity Model

As we consider only the vertical component of the impulse responses, collected echoes are assumed to be mainly associated with P-waves (14). Thus, we adopt in the following a homogeneous P-wave velocity model. More precisely, for each depth, we define a homogeneous velocity distribution whose value is calculated on the basis of the four-layer large scale velocity model (23) (Fig. 1C). This value ranges from c 0 = 2000 m.s -1 at shallow depth to c 0 = 4300 m.s -1 at depth z = 10 km below the surface. The detailed evolution of the wave velocity model c 0 (z) with respect to depth is given in Figure 1C. The assumed background velocity model is rough but is, as we will see, sufficient to image the volcano by leveraging the matrix approach.

Confocal Redatuming

In a first step, the velocity model is used to back-propagate in-depth the recorded echoes gathered in the canonical reflection matrix R in order to retrieve local reflectivity information at each depth of interest. Back-propagation is commonly achieved by applying appropriate time delays at emission and at reception to migrate echoes in post-processing. Such focusing operations are frequently used in imaging and are in particular known as redatuming in seismic exploration (24). The matrix formalism offers a convenient framework to easily perform such beamforming in post-processing, especially in the frequency domain where these operations are described using simple matrix products (10-12, 25) (section S2).

The result is a focused reflection matrix R ρρ (z) = [R(ρ out , ρ in , z)] at each depth z that contains the inter-element impulse responses between a set of virtual sources at r in = (ρ in , z)

and virtual receivers at r out = (ρ out , z) mapping the inner structure of the volcano (Fig. 1B). Its diagonal elements are associated with coincident input and output focusing points (ρ in = ρ out , see section S3). After compensation of wave attenuation with depth (section S4), a 3D confocal image of the volcano is obtained (Fig. 2A) with horizontal cross-sections shown for different depths in Fig. 2B : (i) z = 1.6 km i.e where the most abundant seismicity occurs at La Soufrière;

(ii) z = 6.9 km i.e at the level of the magma reservoir whose depth range is expected between 5.6 and 8.5 km (26); (iii) z = 9.1 km i.e beyond the magma reservoir.

Whether it be on the transverse or the vertical view (see Supplementary Movies 1, 2 and 3), some scattering structures seem to emerge at different locations in Figs. 2A,B but the overall structure appears to be fully blurred, suggesting a high level of aberrations. Such a raw confocal image is indeed very sensitive to aberrations and its interpretation should be extremely cautious.

Focusing Quality

The focusing quality can actually be assessed by considering the off-diagonal elements of R ρρ (z) (ρ in ≠ ρ out ) that provide an estimator of the point spread function in reflection (RPSF, see sections S5 and S6) as a function of the relative position ∆ρ = ρ out -ρ in (11,12). Figure 2C displays the evolution of RPSF for different depths and highlights a significant spreading of the back-scattered energy over off-diagonal elements of R ρρ (z). This is a direct manifestation of the gap between the wave velocity model and its real distribution in the volcano. In absence of aberration, all the back-scattered energy would be contained in a diffracted-limited confocal spot (white circle in Figure 2C) whose size is governed by the angle θ u = tan -1 (d || /2z) under which the geophone array is seen by the focusing point:

δρ u = λ/ (2 sin θ u ) . (1)
In Fig. 2C, the focused wave-field spans over a much larger area than this ideal focal spot and significant side lobes appear around the main central lobe, indicating that images suffer from a high level of aberration.

Overcoming Aberrations

To isolate and compensate for these aberration effects, we build upon a physical phenomenon referred to as the memory effect in wave physics (27). Waves produced by nearby points inside a complex medium can generate highly correlated, but tilted, reflected wave-fields (Fig. 3A).

To exploit this tilt-tilt memory effect, our strategy is thus the following (28) (section S7): (i)

project the reflection matrix between the focused basis (r) and the Earth surface basis (u) (Fig. 3A); (ii) highlight the angular correlations of the reflected wave-field by building a dualbasis matrix (the distortion matrix D) that connects any input focal point in the medium with the distortion exhibited at the Earth surface by the corresponding reflected wavefront (Fig. 3B);

(iii) take advantage of the angular correlations between those wave distortions to accurately estimate the aberration phase transmittance in the Earth surface basis through an iterative phase reversal algorithm (Fig. 3C); (iv) phase conjugate the resulting transmittance to tailor adaptive focusing laws that shall compensate for the volcano's heterogeneities (Fig. 3D).

Vertical and horizontal cross-sections of the resulting confocal image are displayed in Figs. 3E andF, respectively. The comparison with the initial image demonstrates the benefit of the correction process, especially at shallow depths (z < 4 km) where the twisted conduit of the volcano is revealed. The comparison of the original and the corrected RPSFs (Figs. 2C and3G) confirms that the focusing quality is significantly improved in this depth range: Whereas the original RPSF (top panel in Fig. 2B) spreads far beyond the diffraction-limited focal spot, the transverse extension of the corrected RPSF is drastically reduced. However, the gain in image and focusing quality is more modest at larger depths (Figs. 3E,F). The RPSFs still exhibit secondary lobes, a manifestation of residual aberrations (Fig. 3G). Moreover, the spatial extension of the central lobe is limited by the geophone network aperture (Eq. 1). As a consequence, the deep plumbing system of the volcano, in particular the deepest regions of the transcrustal magmatic system and its magma storage zones beyond 5 km depth, cannot be resolved.

Beating Diffraction

Strikingly, an analysis of wave distortions from the k-space will allow us to break this fundamental limit. In the plane wave basis, each distorted wave-field corresponds to the diffracted patterns of each laterally de-scanned output focal spot (4). In a far-field approximation, the contribution of each scatterer would emerge onto limited parts of the k-space because of the finite size of the geophone array (see Fig. 4A and section S8). However, the focal spots also exhibit a parabolic phase law scaling as exp(jk 0 |∆ρ| 2 /z) and resulting from the curvature of focused wave-fronts (see Fig. 3H and section S5). Projected in the k-space, the associated transfer function is thus a superposition of Fresnel rings associated with each reflector. The support of those Fresnel rings is not limited by the geophone network aperture (k 0 sin θ u ) but covers the whole diffraction disk of radius k 0 (Fig. 4B). Iterative phase reversal applied to the D-matrix expressed in the k-space leads to a focusing law (Fig. 4C) that realigns the phase of each spatial frequency component such that the focal spot size reduces to the diffraction limit δρ 0 ∼ λ/2 (sec- wave distortions (28) exhibiting a shift-shift memory effect (27). More importantly, the RPSF extension is now of the order of λ/2 over the whole considered depth range, thereby beating the usual aperture-limited resolution (Eq. 1) displayed by conventional imaging methods.

One necessary condition for this striking performance is the sparsity of the volcano reflectivity with only a few reflectors emerging at each depth (Fig. 4E). As the signature of each reflector is independent, we are able to focus simultaneously on each scatterer provided that they are not too numerous. More precisely, the contrast of the confocal image will typically scale as the ratio between the number of independent geophones and the number of reflectors lying at each depth. Unveiling the plumbing system of La Soufrière Transcrustal magmatic systems consist of vertically-arranged piles of lenses of magmatic mushes more or less ductile (intricate network of crystals and intersticial melt fraction), eruptible melt, and magmatic fluids that extend laterally. This model of a magmatic plumbing system has been described at many other volcanoes (32,33). The internal image of the volcano revealed by Fig. 5 strikingly matches the complex structure described by recent studies on La Soufrière of Guadeloupe (21,22,34,35).

Last but not least, the seismic confocal image of La Soufrière shows that the main magmatic plumbing systems extends from about 5 km below the surface to a depth of at about 8.5 km, values in agreement with this those determined by independent petrological studies (26) who showed that, for the last magmatic eruption of La Soufrière in 1530 CE, the top of the magma storage zone was located between 5.6 and 7 km and the base could not exceed 8.5 km.

Discussion

Inspired by pioneering works in optical microscopy (6,[START_REF] Najar | [END_REF], ultrasound [START_REF] Lambert | Proc. Nat. Acad. Sci. USA[END_REF][START_REF] Bureau | [END_REF] and passive seismology (10)(11)(12), a novel matrix approach to volcano imaging is proposed in this paper. Exploiting Green's functions retrieval from seismic noise (8,9), it takes advantage of bulk seismic waves reflected by heterogeneities in depth to map the internal structure of La Soufrière of Guadeloupe. This volcano is currently in a state of significant long-lasting unrest that could result from dynamic changes in the magmatic plumbing system. Such changes can propagate to the shallow hydrothermal system and could herald eruptive unrest in the future.

The seismic matrix image of the magmatic plumbing system reveals, for the first time, high-resolution features of the magma storage zone, its geometry and dimensions, its complex layered structure, its relative connectedness with other regions of the multi-layer transcrustal magmatic system, and the size and geometry of the upper final eruptive conduit.

The impedance contrast in this complex image also offers the potential, upon further analysis, to distinguish zones of mush from those of eruptible melt, their relative volume, their position in the system. Hence, it can lead to the estimation of parameters such as pressure, temperature, volatile saturation, density contrast, and the connectivity to the surface in evolving magmatic systems, parameters that drive volcanic eruptions.

The strength of this new imaging method lies in its robustness with respect to sparsity of the geophone array and inaccuracy of the wave velocity model. In the future, it will be combined with time-lapse ability resulting from reiteration surveys at active unresting volcanoes and can be coupled with multiparameter data analysis from other classic monitoring networks. Matrix imaging can therefore become a revolutionary game changer in the way scientists understand and model volcanic systems and how volcano observatories monitor their evolving dynamics to forecast their potential for hazardous eruptive activity that threatens the lives of 800 million people living within 100 km from a dangerous volcano (36). responses forms the canonical reflection matrix R gg (t) = [R(g out , g in , t)] that is used to image the inner structure of La Soufrière volcano in the accompanying paper.

S2 Broadband focused reflection matrix

To that aim, a temporal Fourier transform is first applied to R gg (t) to get the set of monochromatic canonical reflection matrices R gg (f ) over the desired frequency bandwidth [10 -20] Hz. The monochromatic matrices are then propagated at emission and at reception towards a focal plane at depth z using the corresponding free-space Green propagator G 0 (z, f ):

G 0 (ρ, g, z, f ) = e -i2πf √ ∥ρ-g || ∥ 2 +∥z-gz∥ 2 /c 0 (z) 4π √ ∥ρ -g || ∥ 2 + ∥z -g z ∥ 2 . ( S1
)
G 0 (z, f ) describes the causal 3-D propagation of waves between any geophone g = (g x , g y , g z ) = (g || , g z ) and any focusing point ρ = (x, y) in the focused basis at depth z in a supposed homogeneous medium with a wave velocity c 0 (z). The evolution of the wave velocity c 0 (z) with respect to depth is provided in Fig. 1C of the accompanying paper. Within the framework of matrix imaging, the projection of R gg (f ) at each depth z is described by the following matrix product:

R ρρ (z, f ) = G * 0 (z, f ) × R gg (f ) × G † 0 (z, f ) (S2)
or in terms of matrix coefficients :

R(ρ out , ρ in , z, f ) = ∑ gout G * 0 (ρ out , g out , z, f ) ∑ g in R(g out , g in , f ) G * 0 (ρ in , g in , z, f ) (S3)
where the symbols * , † and × stand for phase conjugate, transpose conjugate and matrix product respectively. It leads to the set of monochromatic focused reflection matrices

R ρρ (z, f ) = [R(ρ out , ρ in , z, f )].
Physically, each coefficient of R ρρ (z, f ) contains the inter-element impulse response between a virtual source located at r in = (ρ in , z) and a virtual detector at r out = (ρ out , z) (see Fig. 1B of the accompanying paper). In the single scattering regime an in absence of reverberations, the axial dimension of this focal spot, δz 0 , is only limited by diffraction, such that

δz 0 ∼ 2λ z / sin 2 θ u (S4)
with λ z = c(z)/f , the wavelength at depth z.

In order to enhance this axial resolution, a broadband focused reflection matrix R ρρ can be derived at each depth by coherently summing the monochromatic matrices over the frequency bandwidth:

R ρρ (z) = ∫ f + f - df R ρρ (z, f ) (S5)
with f ± = f 0 ± ∆f /2, f 0 = 15 Hz, and ∆f = 10 Hz. The operation amounts to a ballistic time gating of singly-scattered echoes at times t ∼ 2z/c 0 (z). Thanks to this operation, the axial dimension of virtual geophones is greatly reduced and only limited by the frequency bandwidth ∆f :

δz 0 ∼ c 0 (z)/∆f. (S6)
In the single scattering regime, the coefficients of R ρρ (z) can be theoretically expressed as follows (25):

R(ρ out , ρ in , z, f ) = ∫ dρH(ρ, ρ out , z)γ(ρ, z)H(ρ, ρ in , z) (S7)
where γ(ρ, z) is the medium reflectivity at depth z. H(ρ, ρ in/out , z) corresponds to the pointspread-function (PSF), that is to say the spatial amplitude distribution of the focal spot around the focusing point r in/out . Its support defines the characteristic size of each virtual source at r in and detector at r out .

S3 Confocal Imaging

A confocal image of the medium can be easily retrieved from the focused reflection matrix at each depth z by considering the diagonal elements which verify ρ c = ρ in = ρ out :

I(ρ c , z) = R(ρ c , ρ c , z) (S8)
Injecting Eq. S7 into the last equation leads to the following expression for the confocal image:

I(ρ c , z) = ∫ dρH 2 (ρ, ρ c , z)γ(ρ, z) (S9)
Each line of the confocal image results from the sample reflectivity γ and the confocal PSF H 2 . The raw confocal image is displayed in Fig. 2 but note that a time gain compensation is priorly applied to compensate for geometrical spreading as well as scattering and absorption losses that drastically decreases the backscattered energy as a function of depth.

S4 Depth gain compensation of the confocal image

The raw confocal image exhibits a strong amplitude drop with depth (Fig. S1). This attenuation is due to the decay of energy experienced by seismic waves while they propagate. Without compensation, this attenuation strongly degrades the contrast of the confocal image at large depths.

The depth attenuation of the confocal signal can be caused by several factors such as geometrical spreading, scattering and absorption (intrinsic or anelastic attenuation) [START_REF] Shapiro | [END_REF]39). In the present case, the geometrical spreading of waves is compensated, at least partially, by the focusing process performed both at input and output of the reflection matrix. The attenution of the confocal image is thus mainly due to scattering and absorption. In a statistically homogeneous disordered medium, the mean intensity, ⟨|I(ρ, z)| 2 ⟩, shall scale as exp(-2z/ℓ ext ). ℓ ext is the extinction length that combines the scattering and absorption losses as follows:

ℓ -1 ext = ℓ -1 s + ℓ -1
a , with ℓ s , the scattering mean free path and ℓ a , the absorption length.

To retrieve such an exponential decay, the random-like fluctuations of the confocal image due to lateral reflectivity variations should be priorly smoothed out by averaging. The resulting mean confocal intensity, ⟨|I(ρ, z)| 2 ⟩ ρ , is displayed in log-scale as a function of effective depth z in Fig. S1. It highlights four depth ranges with distinct decay rates. For each depth range, the decrease of the mean confocal intenity is fitted by an exponential curve whose decay provides an estimation of ℓ ext reported in Tab. S1.

Depth range Extinction length 0 -0.5 km 1665 m 0.5 -1.9 km 310 m 1.9 -3.8 km 765 m 3.8 -10.5 km 3070 m Table S1: Extinction length. Estimation of the extinction length ℓ ext from the depth decay of the confocal intensity displayed in Fig. S1.

Beyond the near-field (0 -0.5 km) that is probably be polluted by surface waves, the confocal intensity decays rapidly with an exponential fit (red dashed line in Fig. S1) yielding a relatively small value of ℓ ext ∼ 310 m , a manifestation of the strong heterogeneities encountered by the seismic waves in the upper part of the volcano. From depth z = 1.9 km to 3.8 km, this decay slows down. As a consequence, the fitting curve (green dashed line in Fig. S1) is associated with a larger value for the extinction length: ℓ ext ∼ 760 m. Beyond z = 3.8 km (purple dashed line in Fig. S1), the mean confocal intensity decays even more slowly and the value of ℓ ext increases drastically: ℓ ext ∼ 3075 m. In this region (magma storage zone), there is more back-reflection at interfaces between rocks and fluids but less scattering in average (see Fig. 4D and Fig. 5). This point will be discussed in more details in Sec. S10

The overall fitting curve, exp[-β(z)], displayed in Fig. S1, can be used to normalize at each depth the confocal images shown in the manuscript, such that:

I N (ρ, z) = exp[β(z)/2]I(ρ, z), ( S10
)
with I N (ρ, z), the normalized confocal image displayed in the accompanying manuscript. Despite our effort to improve the image contrast by compensating attenuation, the raw confocal image displayed in Fig. 2A appears to be fully blurred, suggesting a high level of aberrations. To understand the effects of diffraction and aberrations on this confocal image, the theoretical expression of the transmit PSF is now investigated.

S5 Transmit point-spread function

The PSF H(ρ, ρ in/out , z) can be expressed using the real Green's function G(ρ, g i , z) between the geophones and the focused basis (25):

H(ρ, ρ in/out , z) = ∑ g i G(ρ, g i , z, f 0 )G * 0 (ρ in/out , g i , z, f 0 ) (S11)
This discrete equation can be rewritten under a continuous form as a function of a coordinate u describing the Earth surface at a depth origin z = 0 defined by the average elevation of the seismic stations g i :

H(ρ, ρ in/out , z) = ∫ duO(u)G(ρ, u, z, f 0 )G * 0 (ρ in/out , u, z, f 0 ) (S12) with O(u) = N -1 g ∑ i δ(u-g ||,i
), the distribution of geophones, N g , their number and δ the Dirac distribution.

In the absence of aberrations, i.e if the wave velocity model is valid (G ≡ G 0 ), the expression of the PSF becomes under the Fresnel approximation:

H 0 (ρ, ρ in/out , z) = exp [i k 0 2z (|ρ| 2 -|ρ in/out | 2 )] O ( ρ -ρ in/out λz ) (S13)
In absence of aberration, the reference PSF H 0 is thus the product between: (i) a geometrical spreading term; (ii) a parabolic phase law that accounts for the curvature of the focused wavefront; (iii) the Fourier transform O of the geophone network aperture O, such that

O ( ρ-ρ in/out λz ) = ∫ duO(u) exp [i k 0 z u.(ρ -ρ in/out )].
The transverse dimension of the focal spot, δρ u , is then only limited by diffraction:

δρ u (z) ∼ λ z /[2 sin[θ u (z)] (S14)
with θ u , the mean angle under which the geophone network is seen by the focusing point r in/out . In the presence of aberrations, i.e., if the velocity model is inaccurate, there is a mismatch between the true Green's matrix G and its model G 0 . If aberrations are moderate (δρ < √ λz), they can be accounted for, at each depth z, by a phase screen A(u, z) at the Earth surface, such that

G(ρ, u, z, f 0 ) = A(u, z)G 0 (ρ, u, z, f 0 ) (S15)
Equation S12 then simplifies into:

H(ρ, ρ in/out , z) = ∫ duF (u, z)G 0 (ρ, u, z, f 0 )G * 0 (ρ in/out , u, z, f 0 ) (S16)
with F (u, z) = O(u)A(u, z), the overall transmittance that combines the array aperture O and the aberration phase screen A. Under the Fresnel approximation,the previous equation becomes:

H(ρ, ρ in/out , z) = exp [i k 0 2z (|ρ| 2 -|ρ in/out | 2 )] ∫ duF (u, z) exp [i k 0 z u.(ρ -ρ in/out )] F ( ρ-ρ in/out λz ) (S17)
The PSF H is thus the product of a parabolic phase law that results from the curvature of focused wave-fronts and a focusing function F , that results from the convolution between the network PSF O that accounts for diffraction and the aberration PSF A defined as the Fourier transform of the aberration transmittance A:

F ( ρ -ρ in/out λz ) = O ⊗ A ( ρ -ρ in/out λz , z) (S18)
where the symbol ⊗ stands for the convolution product.

S6 Reflection point spread function

Interestingly, the focused reflection matrix can provide a local assessment of the focusing quality. Lambert et al. (25) and Touma et al. (11) showed that the amplitude distribution along each antidiagonal of R ρρ (z) provides a key quantity that we will refer to as the reflection pointspread function (RPSF):

RP SF (∆ρ, ρ c , z) = R(ρ c -∆ρ, ρ c + ∆ρ) (S19)
Along an antidiagonal of R ρρ (z), all couple of points on a given antidiagonal share the same midpoint ρ c = (ρ out + ρ in )/2 but with a varying relative position ∆ρ = (ρ out -ρ in )/2. This RSPF is a direct indicator of the local focusing quality. For a sparse scattering medium like a volcano, the RPSF at a scatterer position ρ s scales as:

RSP F (∆ρ, ρ s , z) = exp (i k 0 z |∆ρ| 2 ) F ( ∆ρ λz ) F (- ∆ρ λz ) . (S20)
Therefore, the energy spreading in the vicinity of each scatterer position shall enable one to probe the spatial extension of the PSF. As the scatterer position is a priori unknown, the RPSF is, in practice, probed by considering the antidiagonal whose common mid-point exhibits the maximum confocal signal.

S7 Iterative phase reversal from the Earth surface basis

As highlighted in the previous section, the focused basis is the proper framework for imaging and quantification of focusing quality. However, a dual basis is a better framework to analyse and compensate for aberrations. In the accompanying paper, aberrations are unscrambled by projecting the reflection matrix between the Earth surface (u) and focused basis (r).

The broadband focused reflection matrix can be projected in a geophone basis (here at ouput) by the following matrix product:

R uρ (z) = G ⊺ 0 (z, f 0 ) × R ρρ (z) (S21)
where the symbol ⊺ stands for matrix transpose. An angular de-scan of the input focusing points as sketched in Fig. 3B can be performed by a Hadamard product between R uρ and its ideal counterpart G 0 (z):

D uρ (z) = G † 0 (z, f 0 ) ○ R u,ρ (z) (S22)
Each column of the resulting distortion matrix D uρ (z) maps the phase-distortions withe respect to the ideal wave-front that would be obtained for a point-like source at r in . Injecting Eqs. S7 and Eq. S1 leads to the following expression of the D-matrix coefficients under the Fresnel approximation:

D(u out , ρ in , z) = F (u out , z) transmittance ∫ dρ γ(ρ, z) exp (i k 0 z |ρ| 2 )
apparent reflectivity

F ( ρ -ρ in λz ) transmit PSF exp (i k 0 z u out .(ρ -ρ in ))
angular de-scan (S23) In previous works, the aberration transmittance was extracted through a singular value decomposition (SVD) of D uρ or, equivalently, an eigenvalue decomposition (EVD) of the correlation matrix C uu = D uρ D † uρ (11,28). This result can be understood if we assume a point-like transmit PSF in the previous equation. In that case, we have

D(u out , ρ in , z) = F (u out )γ(ρ in , z) exp (ik 0 |ρ in | 2 /z) .
The D-matrix is then of rank 1 and its first singular vector U 1 directly provides the aberration transmittance at depth z.

Physically, the first eigenvector U 1 is the result of a virtual iterative time reversal (ITR) experiment on a guide star whose reflectivity corresponds to the transmit PSF intensity (28).

Mathematically, the time reversal invariant can be found by solving the following iterative relation (41):

λ (n+1) 1 U (n+1) 1 = C uu × U (n) 1 (S24) with U 1 = lim n→∞ U (n) 1 and λ 1 = lim n→∞ λ (n)
1 , the first eigenvalue of C uu . ITR converges towards the wave-front that maximizes the energy backscattered by the virtual guide star. If this guide star is point like, U (n) 1 thus converges towards the aberration transmittance. However, in reality, the transmit PSF is of course not point-like and its blurring biases the estimation of the aberration transmittance with a time reversal invariant that concentrates on the central part of the geophone array and vanishes on its edge (28).

To circumvent this problem, a related approach consists in an iterative phase reversal (IPR) process [START_REF] Najar | [END_REF][START_REF] Bureau | [END_REF] that forces a uniform amplitude for the phase reversal invariant W, such that

W (n+1) = exp (iarg {C × W (n) }) (S25)
with W (0) = [1⋯1] T chosen arbitrarily as a unit wave-front. The resulting wave-front, W = lim n→∞ W (n) , provides a satisfying estimation of the aberration transmittance over the whole array aperture, as shown in Fig. 3D of the accompanying paper.

Mathematically, the success of IPR can be explained by expressing the correlation matrix C uu . In the speckle regime (42) or in a sparse medium made of a few point-like reflectors, its coefficients can be expressed as follows:

C(u out , u ′ out , z) = F (u out , z)F * (u ′ out , z) [F uout ⊗ F ] (u out -u ′ out , z) (S26)
The correlation term of the right hand side results from the Fourier transform of the input focal spot intensity distribution |F (ρ -ρ in /λz) | 2 . This formulation is reminiscent of the Van Cittert Zernike theorem for an aberrating layer, which links the spatial correlation of a wavefield to the Fourier transform of the intensity distribution from the virtual guide stars (here the input focal spots). In other words, the support of the coherence function scales as the inverse of the input focal spot size. Injecting Eq. S26 into Eq. S25 leads to the following equation:

W (u out , z) = exp ⎛ ⎝ jarg ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ F (u out , z) ∑ u ′ out W (u ′ out , z)F * (u ′ out , z) [F uout ⊗ F ] (u out -u ′ out , z) ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ ⎞ ⎠ .
(S27) For a real autocorrelation function F ⊗F , or equivalently, a symmetric input PSF |F | 2 , the solution of the previous equation is

W (u out , z) = F (u out , z). (S28)
If the previous condition is not fulfilled, the estimation of the aberration transmittance suffers from a bias that can be reduced by iterating the aberration correction process, thereby gradually reducing the size of the virtual guide star and flattening the autocorrelation function F ⊗F .

To do so, the phase conjugate of the estimator W is first used as a focusing law to compensate (partially) for wave distortions. An updated focused reflection matrix is obtained through the following relation:

R ρρ (z) = G † 0 (z, f 0 ) × [G 0 (z, f 0 ) ○ W * (z) ○ D uρ (z)] (S29)
The whole process is then iterated to improve the estimation of the aberration transmittance by alternating aberration correction at input and output (4). In practice, two iterations of the aberration correction process were enough to converge in the present case.

At the end of the process, a novel confocal image is obtained by considering the diagonal elements of the updated focused reflection matrix [Figs. 3 E andF]. The fine compensation of wave distortions is highlighted by the RSPF (Eq. S20) deduced from the updated focused R-matrix (Fig. 3G). As expected theoretically, compensation of aberrations in the geophone basis enables the recovery of a resolution only limited by the geophone aperture. As explained in the acompanying paper, this is nevertheless not sufficient to have a contrasted image of volcano in depth.

S8 Iterative phase reversal driven from the k-space

To go beyond, we will now tackle the parabolic phase law exhibited by the transmit PSF (Eq. S17) that was not addressed in the geophone basis and highlighted by the phase of the RPSF in Fig. 3H. To that aim, the R-matrix can be investigated from the k-space.

As previously done for the geophone basis, the first step consists in a projection of the focused reflection matrix in the plane wave basis:

R kρ (z) = T 0 × R ρρ (z), (S30)
with T 0 , the Fourier transform operator

T 0 (k || , ρ) = exp (-ik || .ρ) . (S31)
Then, the D-matrix is built by isolating the difference between each reflected wave-field and the reference wave-field that would be obtained for a point-like guide star at r in :

D kρ (z) = R kρ (z) ○ T * 0 (z) (S32)
To derive an expression for the D-matrix coefficients in the k-space, one can inject Eqs. S7, S13, S30 and S31 into the last equation. Assuming that aberrations have been fully compensated in the geophone basis (F ≡ O and F ≡ O), this expression writes as follows:

D(k out , ρ in , z) = ∫ dρ transmittance K(k out , ρ, z) reflectivity γ(ρ, z) exp (-i k 0 2z |ρ| 2 ) O ( ρ -ρ in λz ) exp (i k 0 2z |ρ in | 2 ) input PSF exp [ik out .(ρ -ρ in )]
lateral de-scan

, (S33) with K(k out , ρ, z) = ⊮ |kout|<k 0 {O(ρ + λzk out ) kout ⊗ exp (iπλz|k out | 2 )} , (S34)
the aperture transmittance projected in the k-space. The symbol ⊮ |kout|<k 0 accounts for the low pass filter operated by diffraction in the spatial frequency domain: Only spatial frequencies k out whose magnitudes are smaller than |k 0 | can propagate into the Earth; higher spatial frequency components are evancescent and cannot probe in depth the Earth beyond a wavelength. If we compare Eq. S33 with the D-matrix in the geophone basis (Eq. S23), the main difference lies in the fact that the aperture transmittance is no longer isoplanatic from the k-space. Indeed, the angular component of the wave-field recorded by the geophone network and induced by one scatterer in the field-of-view strongly depends on its position, as shown by the term O(ρ + λzk out ) in Eq. S34 and highlighted by Fig. 4A. While this would be an issue in the speckle regime (random reflectivity), this property can become an asset in a sparse scattering medium made of a few reflectors: γ(ρ, z) = ∑ s γ s δ(ρ-ρ s (z)). Under this assumption, Eq. S33 becomes:

D(k out , ρ in , z) = ∑ s K(k out , ρ s , z)γ s exp (-i k 0 2z |ρ s (z)| 2 )O ( ρ s (z) -ρ in λz ) exp (i k 0 2z |ρ in | 2 ) exp [ik out .(ρ s (z) -ρ in )] . (S35)
The coefficients of the associated correlation matrix, C kk (z) = D k,ρ (z)D † kρ (z), are given by:

C(k out , k ′ out , z) = ∑ s |γ s | 2 K(k out , ρ s , z)K * (k ′ out , ρ s , z) [O kout ⊗ O] (k out -k ′ out ) (S36)
For analytical tractability, we will consider, in first approximation, the correlation term

O kout ⊗ O as constant: C(k out , k ′ out , z) ≃ ∑ s |γ s | 2 K(k out , ρ s )K * (k ′ out , ρ s ) (S37)
Under this assumption and provided that the scatterers belong to different resolution cells (K † K ∝ I, with I the identity matrix), Eq. S36 has the form of an eigenvalue decomposition of with σ s , the singular values of D k,ρ , and U s = [U s (k out )], the eigenvectors of C kk , or equivalently the output singular vectors of D k,ρ (z). The identification of Eqs. S37 and S38 shows a one-to-one association between each eigenstate of C kk and each scatterer. The singular values σ s are proportional to each scatterer reflectivity γ s , while each eigenvector U s provides the far-field transmittance K(k out , ρ s ) of each scatterer. Figure S2 confirms this conjecture by showing the result of the ITR processing applied to D kρ (z) at depth z = 6.9 km. As seen before, this process is mathematically equivalent to the SVD of D kρ (z). Figure S2A shows the singular value spectrum of D kρ (z) dominated by three singular values eigenvalues. Each corresponding eigenvector U s covers a distinct part of the k-space (Fig. S2D), as predicted by the term O(ρ s + λzk out ) in Eq. S34. The phase of these eigenvectors (Fig. S2B) shows the Fresnel rings corresponding to the parabolic phase term exp (iπλz|k out | 2 ) in Eq. S34.

C kk : C(k out , k ′ out ) = ∑ s σ 2 s U s (k out )U * s (k ′ out ) (S38)
Note that this observation also enables to revisit the results of Ref. (11) that showed a similar feature in the fault area of San Jacinto. As shown in that previous paper, the phase conjugate of each eigenstate can provide the focusing law to image each scatterer (Figure S2F). A compound image can be built by combining the result provided by each eigenstate. However, this approach only provides an image of the main structures of the volcano. As we will see further, it actually fails in highlighting smaller reflectors.

The IPR process that we previously introduced above in the geophone basis can provide a much more complete view of the inner volcano structure. By forcing a transmittance estimation with the same weight over the whole k-space, the resulting focusing law W (Fig. S2G) can address simultaneously all scatterers in the field-of-view. The angular spectrum addressed by this focusing law can be estimated by considering the modulus of the vector C kk × W (Fig. S2H). The comparison with the angular spectrum covered by the three first eigenvectors of C k,k (z) (Fig. S2D) shows the benefit of IPR for tailoring a focusing law operating over the whole kspace. Also applied at input of the R-matrix, the focusing law derived by IPR leads to the confocal image displayed in Fig. S2I. Compared with the ITR process that only focuses on the three main reflectors at the considered depth (Fig. S2E), the IPR algorithm provides a full-field image of the subsoil highlighting six main structures at the same depth (Fig. S2I).

One could argue that this difference comes from the fact that we did not consider enough eigenstates in the ITR process. However, the higher-order eigenstates cannot be used for imaging purposes (11). Figure S3C illustrates this fact by showing the fourth eigenstate U 4 of D kρ . Its support emerges in the same angular range as U 2 . The link between the second and fourth eigenstates is confirmed by Fig. S3D that displays the phase difference between U 4 and U 2 . U 4 is an higher-order eigenstate associated with the same reflector as U 2 . Each reflector gives actually rise to a set of eigenmodes induced by the autocorrelation term in Eq. S36. Only the fundamental modes corresponding to the highest singular values can be considered for imaging. Higher-order eigenmodes as U 4 in Fig. S3 correspond to smaller singular values and give rise to distorted PSFs (28). They cannot be used for imaging but they pollute the singular value spectrum of D kρ . Hence they can prevent from imaging scatterers of smaller reflectivity within the framework of an ITR process. The proposed IPR process allows to circumvent this limit by finding a phase reversal invariant over the whole angular spectrum (Fig. S2G).

S9 Overcoming diffraction to attain a half-wavelength resolution

As highlighted by the final confocal image (Fig. 4E andF) and corresponding RPSF (Fig. 4G and F), the IPR algorithm driven from the k-space leads to a resolution of the order of λ/2 much thinner than the usual diffraction limit dictated by the geophone aperture: δρ u ∼ λz/d 0 . Mathematically, this can be understood by the convolution product between the geophone aperture and Fresnel rings exhibited by the transmittance K(k out , ρ s ) in the k-space (Eq. S34). These Fresnel rings originate from the parabolic phase law exhibited by the focal spots in real space (Eq. S20). Encoded in the secondary lobes of the PSF, this Fresnel phase law exhibits spatial

S10 Depth evolution of the maximum confocal signal

Figure S4 shows the depth evolution of the maximum confocal signal at each depth without normalization. This curve is plotted at the end of the matrix imaging process. The compensation of aberrations and diffraction operated by matrix imaging allows a finer analysis than the preliminary study provided in Sec. S4. It shows different behaviors in each main part of the volcano. While a strong attenuation is observed in the superior part of the volcano (z =0-3.5 km, see also Sec. S4), the deeper part of the volcano (Fig. S4B) exhibits fluctuations around a relatively constant reflectivity. The magma storage zone (z=5-8.5 km) shows a weaker reflectivity probably due the presence of extended magma volumes. This region is surrounded by two areas of larger reflectivity around 4 and 9 km. This larger confocal signal is probably a manifestation of the important impedance mismatch existing at the boundaries of the magma storage zone (fluid-rock interface). Above the outer carapace of the magma storage system, there are lenses of supercritical acid fluids/brines, and then closer to the surface these supercritical fluids become zones with gases and/or liquid hydrothermal fluids that are present in the pores of the host-rocks, along special zones of porosity-permeability (20). On the one hand, this porous region may account for the enhancement of the confocal signal observed between z = 3.5 and 5 km in Fig. S4. On the other hand, the increase of reflectivity at the bottom of the magma storage system (z = 8.5 -10 km) is a priori due to the back-reflection echo induced by the interface between eruptive melt of the magma storage system and the deeper host-rock.

Figure 1 :

 1 Figure 1: Passive imaging of La Soufrière volcano. (A) Map of the 76 geophones installed above La Soufrière (Guadeloupe, France). Both permanent stations (red) and temporary nodal array (green) are used. (B) Covariance matrix of seismic noise acquired during 2 months is post-processed to obtain the impulse responses between a set of virtual geophones identified by their position r in/out = (ρ in/out , z) and mapping the inside of the volcano. (C) 1D-velocity model (23) used for the seismic data redatuming process.

Figure 2 :

 2 Figure 2: Confocal redatuming (A) Vertical slice of the 3D confocal image along the South-North direction. This image is shown after depth compensation of seismic wave attenuation (section S4). (B) Horizontal slices and (C) associated RPSFs at depths z = 1.6 km, 6.9 km and 9.1 km below the summit. The spatial extension δρ u (Eq. 1) of the theoretical diffraction-limited focal spot is denoted as a white circle.

Figure 3 :

 3 Figure 3: Exploiting memory effect for overcoming aberrations. (A) Schematic view of back-scattered wave-fronts in the Earth surface basis (u) generated by virtual sources (r in ) at a given depth z. (B) The extraction of wave-distortions amounts to an angular de-scan of each input focal spot. (C) Iterative phase reversal applied to the D-matrix at each depth provides an aberration transmittance whose phase is here shown at depths z = 1.6 km, 6.9 km and 9.1 km and whose modulus is encoded as a transparency mask. Black dots indicate the lateral position of geophones. (D) Phase conjugation and tilt of such aberration phase laws enable an adaptive focusing process on each point of the subsoil. (E) Vertical slice of the 3-D confocal image obtained by means of these optimized focusing laws. This image is shown after a depth compensation of seismic wave attenuation (section S4). (F) Corresponding horizontal slices at depths z = 1.6 km, 6.9 km and 9.1 km. (G) Modulus and (H) phase of the resulting RPSFs at the same depths. The spatial extension δρ u (Eq. 1) of the theoretical diffraction-limited focal spot is denoted as a white circle.

Figure 4 :

 4 Figure 4: Overcoming diffraction by operating the focusing process from the k-space. (A) The distorted wave-field in the k-space amounts to a lateral de-scan of each reflected echoes. Under the Franhofer approximation, the far-field projection of each focal spot is limited by the de-scanned geophone network aperture. (B) In the present case, Fresnel diffraction gives rise to a modulation of each focal spot by a parabolic phase law. The diffraction pattern of each reflector then corresponds to Fresnel rings that cover the whole diffraction disk of radius k 0 . For sake of clarity, only the contribution of a central reflector is displayed. (C) Iterative phase reversal applied to D extracts those diffraction patterns here shown at depths z = 1.6 km, 6.9 km and 9.1 km. (D) Vertical slice of the resulting image (same view as in Figs. 2A and 3D). This image is shown after depth compensation of seismic wave attenuation (section S4). (E) Corresponding horizontal slices and (F) RPSFs at the same depths as panel C. The spatial extension δρ u (Eq. 1) of the aperture-limited focal spot is denoted as a white circle.

  tion S9). It leads to a new confocal image whose several cross-sections are displayed in Figs. 4D and E. In particular, a complex multi-lens melt reservoir is revealed by Fig 4D beyond a depth of 5 km. The comparison with the previous image (Figs. 3E and F) highlights the spectacular gain in terms of contrast and resolution provided by a k-space analysis of the D-matrix. This observation is confirmed by the new RPSFs displayed in Fig. 4F. Compared to their previous version (Fig. 3G), the diffuse background has been suppressed by a compensation of residual

Figure 5 :

 5 Figure 5: Three-dimensional view of the first 10 km of the hydrothermal and magmatic system of La Soufrière. (A,B) Isosurface plots of the three-dimensional image of the volcano viewed from East and North, respectively. The isosurface is fixed to be -15 dB. (C,D) Corresponding zooms on the first 3 km depth. The isosurface is fixed to be -10 dB. This image is shown after a depth compensation of seismic wave attenuation (section S4).

Figure 5

 5 Figure5shows two perpendicular views of La Soufrière down to a depth of 10.5 km below the summit (see also Supplementary Movies 1, 2 and 3). Based on an analysis of the P-wave reflected wave-field, it displays the iso-surfaces of the confocal image obtained at the end of the matrix imaging process. As outlined above, the superior part of the volcano, from a depth of 5 km up to the surface, exhibits the clear signature of a tortuous conduit that finds its way through the hostrock forming the upper part of the volcano. On the contrary, its deep structure, between ca. 5 and 8.5 km depth, induces a more diffuse scattering that is compatible with the existence of a vertical succession of several sub-horizontal and irregular globular coalescing structures. Those elements are superimposed over a distance of a few kilometers and linked together by narrow sub-vertical diffuse structures. The sub-horizontal structures extend

  figures. A.B. and J.-C. K. provided the geophysical and volcanological interpretation for the obtained image. E.G., J.-C. K. and A.A. prepared the manuscript. E.G., A.B., A.L.B., J.-C. K., M.F., and A.A. discussed the results and contributed to finalizing the manuscript.

Figure S1 :

 S1 Figure S1: Scattering intensity vs. depth. The logarithm of the mean confocal intensity, ⟨|I ′ (ρ, z) 2 ⟩ ρ , is plotted as a function of depth. The depth decay of intensity is fitted with an exponential curve over four depth ranges summarized in Tab. S1.

Figure S2 :

 S2 Figure S2: Comparison between iterative time reversal (left) and iterative phase reversal (right) processes for aberration phase law extraction -Illustration at depth z = 6.9 km. (A) Singular value histogram of the distortion matrix D kr . (B) Phase and (C) modulus of the first three singular vectors U s = [U s (k out )]. (D) Sum of the modulus of the first three eigenvectors. (E) Final image resulting from the combination of each image displayed in panel F. (F) Image obtained using the phase conjugate of U s as the focusing law. (G) Phase ϕ of the IPR invariant W and (H) associated angular spectrum obtained by considering the modulus of C kk × W. (I) Confocal image using the phase conjugate of W as the focusing law.

Figure S3 :

 S3 Figure S3: Limit of iterative time reversal processing -Illustration at depth z = 6.9 km. (A) Singular value histogram of the distortion matrix D kr . (B) Modulus and (C) phase of the second and fourth singular vectors, U 2 and U 4 . (D) Phase of the Hadamard product between U 2 and U * 4 .

Figure S4 :

 S4 Figure S4: Depth evolution of the maximum confocal intensity at each depth: (A) between z = 1 and 10.5 km; (B) between z = 3 and 10.5 km.

Acknowledgments

The authors wish to thank R. Touma, M. Campillo and A. Derode for initial discussions on the project; the colleagues at the OVSG-IPGP for field assistance in installing and running the node network, SISMOB network of RESIF (now EPOS-FRANCE); the Parc National de Guadeloupe for permission to install the node in the field.

Funding

The authors are grateful for the funding provided by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 819261, REMINISCENCE project). This work has also been supported by the AO-IPGP 2017 project "Dense seismic monitoring of the hydrothermal system of La Soufrière de Guadeloupe", the AO-TelluS-INSU 2017 action ALEAS (coord.: A. Burtin), the project "Vers la Plateforme Régionale de Surveillance Tellurique du futur" -(PREST) co-funded by INTERREG Caraïbes V for the European Regional Development Fund, and the European Union's Horizon 2020 research and innovation programme (grant agreement, no. 731070, EUROVOLC project). The authors also thank IPGP for general funding to the Observatoires Volcanologiques et Sismologiques (OVS), the INSU-CNRS for funding provided by the Service National d'Observation en Volcanologie (SNOV), and the Ministère pour la Transition Ecologique et Solidaire (MTES) as well as the IdEx project "Université Paris Cité" (ANR-18-IDEX-0001) for financial support.

Data and Materials Availability

Seismic data used in this manuscript has been deposited at the Data collection of the seismological and volcanological observatory of Guadeloupe [START_REF]Institut De Physique Du Globe De Paris (IPGP[END_REF][START_REF] Burtin | Dense seismic monitoring of La Soufrière de Guadeloupe hydrothermal system[END_REF]. The seismic noise correlation data generated in this study and codes used to post-process it are available at Zenodo [START_REF] Giraudat | Unveiling the deep plumbing system of a volcano by a reflection matrix analysis of seismic noise[END_REF] (https://zenodo.org/record/10066910).

Author Contributions

Competing Interests

The authors declare no competing interests.

Supplementary materials

Supplementary Text

Figs. S1-S2-S3-S4

Reference [START_REF] Shapiro | [END_REF](39)(40)(41)(42)

Supplementary Information

This document provides further information on: (i) the seismic data and noise correlation processing; (ii) the building of the focused reflection matrix; (iii) the formation of the confocal image ; (iv) its depth gain compensation for display purpose; (v) the transmit point spread function; (vi) the reflection point spread function; (vii) the iterative phase reversal algorithm for compensation of wave distortions from the Earth surface basis; (viii) iterative phase reversal driven from k-space and its comparison with a singular value decomposition approach; (ix) the spatial resolution of the final image; (x) the depth evolution of the maximum confocal signal.

S1 Seismic data and noise correlation processing

The seismic data used in this study consists of a temporary nodal array of 65 geophones ( 16) and 6 permanent stations [START_REF]Institut De Physique Du Globe De Paris (IPGP[END_REF] operated by the OVSG-IPGP (Volcanologic and Seismologic Observatory of Guadeloupe). The geophone sensors were Zland 3C Gen2 (Fairfieldnodal) with a natural frequency of 5 Hz, recording at 500 samples per second and along 3 orthogonal directions (Vertical, North and East). The 6 OVSG seismic stations are 3 components broadband sensors, all having a flat response in the [1-50] Hz frequency band. The seismic records are sampled at 100 Hz. For this study only vertical components are used. The temporary nodal array was deployed from mid-November 2017 to mid-January 2018 during 2 sessions in order to download seismic data and recharge batteries. Since we moved the location of 5 geophones between both acquisition sessions, we ended with a virtual network of 76 sites, for which we could applied the computation of seismic noise correlation functions (NCFs).

The procedure to compute the seismic NCF mainly follows the stages detailed in Ref. (40). Here, we summarize each step that we apply on seismic recordings whether it was a temporary geophone or a permanent seismic sensor. (1) We detrend each hourly vertical seismic record and removed the mean. (2) We remove the instrument response to homogenize the seismic signals and we applied a band-pass filter between 1 Hz and 20 Hz. (3) We resample the seismic record to a unique sample frequency of 100 Hz. (4) We apply a spectral and temporal normalization by proceeding to a spectral whitening followed by a 1-bit normalization to only keep the sign of the seismic signal. [START_REF] Kang | [END_REF] We end with the computation of the NCF by cross-correlating hourly seismic records at each stations pair for time delays ranging from -30 to +30 seconds. To increase the signal to noise ratio of a NCF, we apply some quality checks and a waveform summation by first averaging the 24 hourly NCFs in a daily one, for which we discard hourly segments that were not coherent with the raw daily average (correlation coefficient threshold of 0.5). The average over each daily NCF estimated during the 2 months of nodal array deployment finally provides an estimation of the impulse response R(g j , g i , t) between each couple of geophones i and j whose positions are identified by vector g. The set of the estimated 2850 vertical impulse