Xiong Wang

Junfeng Zhang
email: jfzhang@cug.edu.cn

Andrea Tommasi
email: andrea.tommasi@umontpellier.fr

Zhicheng Jing

Maoshan Yuan

Andréa Tommasi

Microstructure and seismic properties of amphibole-rich rocks from the deep crust in southern Tibet

Keywords:

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Seismic anisotropy is a powerful tool to decipher the deformation patterns in the deep Earth. Geophysical observations show that seismic anisotropy is ubiquitous in the lithosphere, asthenosphere, and lowermost mantle (D ′′ ) [START_REF] Ferreira | Ubiquitous lower-mantle anisotropy beneath subduction zones[END_REF][START_REF] Wookey | Mid-mantle deformation inferred from seismic anisotropy[END_REF][START_REF] Wu | Complex seismic anisotropy beneath western Tibet and its geodynamic implications[END_REF]. Seismic anisotropy is commonly ascribed to crystallographic preferred orientations (CPOs) of mica [START_REF] Lloyd | Mica, deformation fabrics and the seismic properties of the continental crust[END_REF] and amphibole [START_REF] Tatham | Amphibole and lower crustal seismic properties[END_REF] in the deep crust, olivine [START_REF] Nicolas | Formation of anisotropy in upper mantle peridotitesa review[END_REF] in the upper mantle, wadsleyite [START_REF] Kawazoe | Seismic anisotropy in the mantle transition zone induced by shear deformation of wadsleyite[END_REF][START_REF] Tommasi | Strain-induced seismic anisotropy of wadsleyite polycrystals and flow patterns in the mantle transition zone[END_REF] in the mantle transition zone, post-perovskite in the D ′′ layer [START_REF] Merkel | Deformation of (Mg, Fe) SiO3 post-perovskite and D ′′ anisotropy[END_REF][START_REF] Tommasi | Deformation, crystal preferred orientations, and seismic anisotropy in the Earth's D ′′ layer[END_REF][START_REF] Wu | Seismic anisotropy of the D ′′ layer induced by (001) deformation of post-perovskite[END_REF]. The basic principle behind the interpretation of deformation patterns from seismic anisotropy data is to be able to link the anisotropies in seismic velocities to CPOs of major rock-forming minerals and the CPOs to a deformation geometry. Thus, we can infer the deformation patterns from seismic anisotropy data if the process producing the CPOs is well known, as in dislocation creep. However, a large body of evidence indicates that strong CPOs can also be produced by other mechanisms, for instance, cataclastic flow [START_REF] Ko | Crystal preferred orientation of an amphibole experimentally deformed by simple shear[END_REF], anisotropic growth during diffusion creep [START_REF] Getsinger | Amphibole fabric formation during diffusion creep and the rheology of shear zones[END_REF][START_REF] Miyazaki | Olivine crystals align during diffusion creep of Earth's upper mantle[END_REF], oriented crystallization in a rock deforming in the presence of fluids [START_REF] Imon | Development of shape-and lattice-preferred orientations of amphibole grains during initial cataclastic deformation and subsequent deformation by dissolution-precipitation creep in amphibolites from the Ryoke metamorphic belt, SW Japan[END_REF][START_REF] Stokes | Deformation of amphibolites via dissolution-precipitation creep in the middle and lower crust[END_REF][START_REF] Giuntoli | Replacement reactions and deformation by dissolution and precipitation processes in amphibolites[END_REF] or melt topotactic reactions [START_REF] Boudier | Serpentine mineral replacements of natural olivine and their seismic implications: oceanic lizardite versus subductionrelated antigorite[END_REF][START_REF] Higgie | Feedbacks between deformation and melt distribution in the crust-mantle transition zone of the Oman ophiolite[END_REF][START_REF] Padrón-Navarta | Oriented growth of garnet by topotactic reactions and epitaxy in high-pressure, mafic garnet granulite formed by dehydration melting of metastable hornblende-gabbronorite (Jijal Complex, Kohistan Complex, north Pakistan)[END_REF][START_REF] Tommasi | Hydrous melts weaken the mantle, crystallization of pargasite and phlogopite does not: Insights from a petrostructural study of the Finero peridotites, southern Alps[END_REF], and crystal rotations and interactions during magmatic flow [START_REF] Berger | Deformation mechanisms and reaction of hornblende: examples from the Bergell tonalite (Central Alps)[END_REF][START_REF] Ildefonse | Effect of mechanical interactions on the development of shape preferred orientations: a two-dimensional experimental approach[END_REF][START_REF] Nicolas | Kinematics in Magmatic Rocks with special Reference to Gabbros[END_REF]. However, in these cases, the rules linking the CPO-induced seismic anisotropies to the deformation patterns may differ from those for CPOs produced by dislocation creep. The determination of processes responsible for the development of CPOs is, therefore, an essential step for the interpretation of seismic anisotropy data, especially in tectonic settings with abundant magmatic activities, such as orogenic belts.

The P to S converted phase from the Moho (PmS) permits the characterization of the azimuthal anisotropy in the crust [START_REF] Liu | Estimating crustal seismic anisotropy with a joint analysis of radial and transverse receiver function data[END_REF]. Strong shear wave anisotropy, quantified by the polarization orientation of the fast shear wave (ϕ) and the delay time (δt) between the fast and slow shear waves, is observed in the crust of southern Tibet (Fig. 1a). The fast shear wave polarization directions are mainly subparallel to the India-Asia suture boundary. The average shear wave delay times are 0.50 ± 0.15 s in central Tibet [START_REF] Wu | Deformation of crust and upper mantle in Central Tibet caused by the northward subduction and slab tearing of the Indian lithosphere: New evidence based on shear wave splitting measurements[END_REF], 0.50 ± 0.40 s in east Tibet [START_REF] Sun | Crustal structure and deformation under the Longmenshan and its surroundings revealed by receiver function data[END_REF], and 0.48 ± 0.13 s in southeastern Tibet [START_REF] Cai | Crustal anisotropy and deformation of the southeastern margin of the Tibetan Plateau revealed by Pms splitting[END_REF][START_REF] Kong | Crustal anisotropy and ductile flow beneath the eastern Tibetan Plateau and adjacent areas[END_REF][START_REF] Sun | Crustal structure and deformation of the SE Tibetan plateau revealed by receiver function data[END_REF][START_REF] Zheng | Crustal azimuthal anisotropy beneath the southeastern Tibetan Plateau and its geodynamic implications[END_REF]. In southwestern Tibet, fast polarizations are at a high angle to the suture and delay times more variable, ranging between 0.45 and 1.3 s [START_REF] Wu | Complex seismic anisotropy beneath western Tibet and its geodynamic implications[END_REF]. Strong radial anisotropy is also observed using surface waves in the deep crust (peak of sensitivity at 35 km) of southeastern Tibet (Fig. 1b) [START_REF] Hu | Direct surface wave radial anisotropy tomography in the crust of the Eastern Himalayan Syntaxis[END_REF]. These studies propose that the observed crustal anisotropy mainly originates from CPOs formed during ductile deformation of the middle-lower crust [START_REF] Agius | Complex, multilayered azimuthal anisotropy beneath Tibet: evidence for co-existing channel flow and pure-shear crustal thickening[END_REF][START_REF] Chen | Low velocity crustal flow and crust-mantle coupling mechanism in Yunnan, SE Tibet, revealed by 3D S-wave velocity and azimuthal anisotropy[END_REF][START_REF] Shapiro | Thinning and flow of Tibetan crust constrained by seismic anisotropy[END_REF]. Amphibole, because of its abundance in the deep crust and high intrinsic seismic anisotropies, has often been proposed to explain strong crustal anisotropies [START_REF] Ko | Crystal preferred orientation of an amphibole experimentally deformed by simple shear[END_REF][START_REF] Li | Petrofabrics and seismic properties of Himalayan amphibolites: implications for a thick anisotropic deep crust beneath southern Tibet[END_REF][START_REF] Tatham | Amphibole and lower crustal seismic properties[END_REF]. However, observations of naturally and experimentally deformed amphibolites show that most of the processes cited above may produce strong amphibole CPOs [START_REF] Elyaszadeh | Different slip systems controlling crystallographic preferred orientation and intracrystalline deformation of amphibole in mylonites from the Neyriz mantle diapir, Iran[END_REF][START_REF] Getsinger | Amphibole fabric formation during diffusion creep and the rheology of shear zones[END_REF][START_REF] Imon | Development of shape-and lattice-preferred orientations of amphibole grains during initial cataclastic deformation and subsequent deformation by dissolution-precipitation creep in amphibolites from the Ryoke metamorphic belt, SW Japan[END_REF][START_REF] Ko | Crystal preferred orientation of an amphibole experimentally deformed by simple shear[END_REF][START_REF] Skrotzki | Defect structure and deformation mechanisms in naturally deformed hornblende[END_REF][START_REF] Stokes | Deformation of amphibolites via dissolution-precipitation creep in the middle and lower crust[END_REF][START_REF] Tommasi | Hydrous melts weaken the mantle, crystallization of pargasite and phlogopite does not: Insights from a petrostructural study of the Finero peridotites, southern Alps[END_REF]. Determination of the processes producing amphibole CPOs in the deep crust of Tibet is therefore important for the interpretation of the seismic anisotropy data.

In southern Tibet, the deep crust underwent a series of deformation, metamorphic and magmatic events during subduction and collision [START_REF] Zhu | Magmatic record of India-Asia collision[END_REF]. Deep crustal amphibolite xenoliths carried by mantlederived magmas usually display strong amphibole CPOs [START_REF] Ji | Magnitude and symmetry of seismic anisotropy in mica-and amphibole-bearing metamorphic rocks and implications for tectonic interpretation of seismic data from the southeast Tibetan plateau †[END_REF], but the processes responsible for these CPOs have not been constrained. Due to fluid release during subduction, amphibole-rich rocks may form by hydration reactions in the deep crust. In such a setting, amphibole CPOs may form by either oriented or topotactic growth. Fractional crystallization of hydrous basaltic magma may also produce amphibolerich cumulates [START_REF] Davidson | Amphibole "sponge" in arc crust?[END_REF][START_REF] Zhou | Amphibole-rich cumulate xenoliths in the Zhazhalong intrusive suite, Gangdese arc: Implications for the role of amphibole fractionation during magma evolution[END_REF]. Due to the strongly anisometric crystal shape of amphiboles, CPOs of amphibole may form during the crystallization of these magmas. These CPOs will show a large-scale (km to tens of km) organization if the magmatic flow is controlled by a large-scale (tectonic) stress field, as observed in synkinematic plutons (e.g., [START_REF] Nicolas | Kinematics in Magmatic Rocks with special Reference to Gabbros[END_REF][START_REF] Vauchez | Slow cooling and crystallization of the roots of the Neoproterozoic Araçuaí hot orogen (SE Brazil): Implications for rheology, strain distribution, and deformation analysis[END_REF].

To further constrain the processes responsible for the observed seismic anisotropies in the middle-lower crust of southern Tibet, we have characterized the compositions, microstructures, CPOs, and seismic properties of amphibole-rich deep-crust xenoliths and exhumed rocks from a deep-crust terrane in southern Tibet.

Geological setting and samples

The Tibetan-Himalayan orogenic system formed by collision between the Indian and Asian continents during the Cenozoic [START_REF] Yin | Geologic evolution of the Himalayan-Tibetan orogen[END_REF]. The Tibetan Plateau is composed of, from south to north, the Himalaya, Lhasa, Qiangtang, and Songpan-Ganzi terranes (Fig. 1a). The Lhasa terrane lies between the Bangong-Nujiang suture zone in the north and the Indus-Tsango suture zone (ITS) in the south. It can be divided into northern, central, and southern subterranes based on different sedimentary cover rocks [START_REF] Zhu | Magmatic record of India-Asia collision[END_REF]. The southern Lhasa subterrane hosts the widespread Cretaceous-early Tertiary Gangdese batholith and the Linzizong volcanic succession [START_REF] Zhu | Magmatic record of India-Asia collision[END_REF]. Three amphibole-rich rocks (#521, #431-1, and #431-2) were collected from enclaves in granitic gneisses of the Gangdese batholith (Fig. S1c-d and see locations in Fig. 1b). Isotopic data indicated these xenoliths may represent cumulates or residuals of basaltic magma ponding at the base of crust [START_REF] Xu | Constructing the early Mesozoic Gangdese crust in southern Tibet by hornblende-dominated magmatic differentiation[END_REF].

The Eastern Himalayan Syntaxis (EHS), located at the eastern part of Himalayan orogen (black box in Fig. 1b), is a strongly deformed and fast exhumed terrane subjected to migmatitization [START_REF] Booth | Constraints on the metamorphic evolution of the eastern Himalayan syntaxis from geochronologic and petrologic studies of Namche Barwa[END_REF]. The EHS includes three major tectonic units: the Lhasa terrane, representing the southern segment of the Asian continent; the ITS, forming the residual of the Neo-Tethyan Ocean between the Asian and Indian plates; and the Himalayan sequences, representing the northern margin of the Indian continent [START_REF] Yin | Geologic evolution of the Himalayan-Tibetan orogen[END_REF]. The Higher Himalayan Crystalline (HHC) unit, locally referred to as the Namche Barwa Complex (NBC), mainly consists of high-grade gneiss, amphibolite, schist, marble, granulite, and migmatite, and is interpreted as an exhumed middle-lower crust, which has been subjected to extensive partial melting [START_REF] Kang | Prolonged partial melting of garnet amphibolite from the Eastern Himalayan Syntaxis: Implications for the tectonic evolution of large hot orogens[END_REF]. Amphibolite occurs interlayered with schist and gneiss and usually shows clear deformation structures (foliation and lineation). A recent study showed that the protoliths of these amphibolites were mafic magmatic rocks of Late Paleoproterozoic age [START_REF] Kang | Prolonged partial melting of garnet amphibolite from the Eastern Himalayan Syntaxis: Implications for the tectonic evolution of large hot orogens[END_REF] and that the amphibolites recorded a clockwise P-T path with peak metamorphism of high-pressure granulite-facies, leading to extensive partial melting associated with amphibole dehydration under conditions of 1.6-1.8 GPa and <800 • C at 40-20 Myr, followed by retrogression to 0.5-0.6 ± 0.1 GPa and 695-740 ± 35 • C between 20 and 10 Myr [START_REF] Kang | Prolonged partial melting of garnet amphibolite from the Eastern Himalayan Syntaxis: Implications for the tectonic evolution of large hot orogens[END_REF]. Three amphibolite samples (16lzwx1-3) were collected from a thick layer of amphibolite in the core of the EHS (see locations in Fig. 1b). [START_REF] Wu | Complex seismic anisotropy beneath western Tibet and its geodynamic implications[END_REF], central south Tibet [START_REF] Wu | Deformation of crust and upper mantle in Central Tibet caused by the northward subduction and slab tearing of the Indian lithosphere: New evidence based on shear wave splitting measurements[END_REF], and east Tibet [START_REF] Kong | Azimuthal anisotropy and mantle flow underneath the southeastern Tibetan Plateau and northern Indochina Peninsula revealed by shear wave splitting analyses[END_REF][START_REF] Sun | Crustal structure and deformation of the SE Tibetan plateau revealed by receiver function data[END_REF][START_REF] Sun | Crustal structure and deformation under the Longmenshan and its surroundings revealed by receiver function data[END_REF][START_REF] Zheng | Crustal azimuthal anisotropy beneath the southeastern Tibetan Plateau and its geodynamic implications[END_REF]. The color bar at the bottom is elevation in kilometers. The white box denotes the location of Figure1b. b) Radial anisotropy in southeastern Tibet at 35 km depth (modified from [START_REF] Hu | Direct surface wave radial anisotropy tomography in the crust of the Eastern Himalayan Syntaxis[END_REF]. The color bar at the bottom shows the range of radial anisotropy. Red and blue colors refer to positive and negative radial anisotropies, respectively. The black box marks approximately the location of the Eastern Himalayan Syntaxis (EHS). Stars indicate locations of studied samples: 16lzwx1-3 are the samples collected in the NBC complex while #521, #431-1, and #431-2 are collections from the xenoliths. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

X. Wang et al.

Methods

The foliation and lineation of exhumed samples from EHS were determined based on the shape preferred orientation (SPO) of amphibole grains and field observations; thin sections were prepared normal to foliation and parallel to the lineation. Thin sections for xenoliths #431-1 and #431-2 were cut in a random direction normal to the compositional layering, as the lineation could not be determined. Xenolith #521 was cut randomly because it is composed of nearly pure amphibole and no SPO can be determined. All thin sections were polished until the sample surfaces were even and clean for microstructural and microprobe analyses. For EBSD mapping, vibration polishing with 0.05 μm colloidal silica was applied for at least 4 h before measurements.

Analysis of mineral compositions

Quantitative compositional analyses of amphibole, plagioclase, and garnet were performed using a JEOL JXA-8230 electron probe microanalyzer at the Center for Global Tectonics, School of Earth Sciences, China University of Geosciences (Wuhan) with an accelerating voltage of 15 kV, a beam current of 20 nA, and a beam size of 1-3 μm. The following standards were used: Sanidine (K), Pyrope Garnet (Fe, Al), Diopside (Ca, Mg), Jadeite (Na), Rhodonite (Mn), Olivine (Si), Rutile (Ti) and Apatite (P).

CPOs data acquisition and processing

Crystal orientation data were obtained using a JEOL JSM 5600 scanning electron microscope equipped with EBSD system at Geosciences Montpellier, Université de Montpellier, France (Samples: 16lzwx1-3) and a Quanta 450 Field Emission Gun (FEG)-SEM-EBSD system at the School of Earth Sciences, China university of Geosciences, Wuhan (Samples: #521, #431-1, and #431-2). Large area mapping of samples was performed to acquire a representative orientation dataset. Conditions during EBSD mapping were an accelerating voltage of 17 kV, a spot size of 5.5-6, and a working distance of 24 mm. The automatic data were collected with a step size of 15-28 μm depending on grain size. Raw indexing rates range between 78 and 85%. The raw data have been noise-reduced using a "wild spike" correction with a threshold mean angular deviation (MAD) values of 1.3 • . Non-indexed pixels with at least six neighboring coherent orientations have been filled with the average orientation of the neighbors using the CHANNEL 5 software.

The crystallographic data analyses were performed using the MTEX toolbox in MATLAB [START_REF] Bachmann | Texture analysis with MTEXfree and open source software toolbox[END_REF][START_REF] Hielscher | A novel pole figure inversion method: specification of theMTEXalgorithm[END_REF]. CPOs data are displayed as one point per grain in pole figures presented as lower hemisphere stereographic projections. The CPOs data of sample #521 was rotated by aligning the clusters of the [001] and [100] axes of amphibole parallel to lineation and normal to foliation of the other sample, respectively. Densities of pole figures were normalized to a uniform distribution. J-index (one point per grain) was used to quantify the strength of CPOs. The automatic grain detection routine in the MTEX toolbox [START_REF] Bachmann | Grain detection from 2d and 3d EBSD geology[END_REF] has been used for quantitative microstructure analyses. Grain boundaries are defined by misorientations between neighboring pixels >10 • . Grains composed by less than 10 pixels were not considered in the microstructural analysis (grain sizes and shapes). The misorientation of each pixel relative to the mean orientation of the grain (M2M) and the grain orientation spread (GOS) was calculated to quantify the intracrystalline orientation gradients, which is a proxy for dislocation density.

Calculation of seismic properties

The seismic properties of amphibole aggregates and amphibolite with varying amphibole/plagioclase ratios were computed using the MTEX tool (https://mtex-toolbox.github.io) with the measured CPOs data of constituent minerals, single-crystal elastic constants, the Christoffel equation, and the Voigt-Reuss-Hill averaging scheme [START_REF] Mainprice | Calculating anisotropic physical properties from texture data using the MTEX open-source package[END_REF]. The elastic constants (C ij ) and densities (ρ) of amphibole [START_REF] Brown | Elasticity of calcium and calcium-sodium amphiboles[END_REF] and plagioclase (Brown et al., 2016) used in the calculation of seismic properties were shown in Table S1. The parameters (C ij and ρ) of Ab 75 An 25 and Ab 52 An 48 were used for plagioclase in samples 16lzwx2 and #431-1 respecting the different chemical compositions of these samples.

Results

Microstructures

Photomicrographs and EBSD phase maps illustrating the microstructures of the six studied amphibolites are displayed in Figs. 2 and3. The volume fractions of constituent minerals and microstructure parameters from EBSD mapping are presented in Table 1. The amphibole content is ≥65 vol% in all samples. Other constituent minerals are plagioclase, garnet, quartz, and minor oxides.

The amphibolites (16lzwx1-3) from EHS are fresh samples from the cores of the lenticular amphibolites, which were only slightly affected by subsequent retrogression during exhumation. The amphibolite 16lzwx1 has a bimodal grain size distribution. It contains coarse-grained garnet and amphibole porphyroblasts enclosed in a medium-grained to finegrained matrix of amphibole and minor plagioclase and quartz (Fig. 3a). Coarse amphibole grains have lenticular shapes (Fig. 2a). The foliation, which is marked by the SPO of the amphibole, molds the garnet porphyroblasts (Fig. 2a). Minor fine-grained plagioclase (10%) and quartz (1%) are mainly located along grain boundaries of amphibole, locally forming patches flattened in the foliation plane (Figs. 2a,3a). The SPO of the coarse amphibole grains also marks a lineation.

The amphibolite 16lzwx2 is garnet-free and contains more plagioclase (20%) and quartz (9%) than sample 16lzwx1. Amphibole and plagioclase have similar mean grain sizes (0.3-0.5 mm) and, despite their irregular shapes, both show a strong SPO (maximum aspect ratio is presented in Table 1), which delineates the foliation and lineation (Figs. 2b,3b). This sample also shows an incipient banding marked by alternating enrichment in amphibole or plagioclase (Fig. 3b). Quartz occurs mainly in flattened aggregates along amphibole grain boundaries.

The amphibolite 16lzwx3 is composed of coarse garnet porphyroblasts (5-8 mm) enclosed in a medium-grained (~0.4 mm) amphibole matrix (Fig. 3c). The garnet grains have nearly rounded shapes and are molded by the foliation formed by the amphibole SPO (Fig. 3c). Amphibole has subhedral to euhedral shapes and shows no evidence of intracrystalline deformation (Fig. 2c). Plagioclase and quartz form pressure-shadow structures around the garnet porphyroblasts, indicating retrogression during exhumation (Fig. 3c). Isolated small plagioclase and quartz grains can also be observed in the matrix.

The three amphibole-rich xenoliths have typical magmatic textures without evidence for either solid-state deformation or metamorphism (Fig. 2d and3d-f). Sample #521 is composed of nearly 100% coarsegrained amphibole (~1.7 mm). It shows a very weak SPO (the trace of the foliation is roughly vertical in Fig. 3d). Sample #431-1 has an equigranular texture. It is composed of amphibole (65%), plagioclase (30%), and minor quartz (3%). Amphibole and plagioclase are coarsegrained (0.4-0.7 mm; Fig. 2d). Amphibole is subhedral to lens-shaped (Fig. 2d). Its SPO, together with the shape of the plagioclase aggregates, defines the foliation (Fig. 3e). Plagioclase grains have irregular shapes, but a weak SPO. They display well-developed polysynthetic growth twins (Fig. 2d). Variations in the amphibole and plagioclase proportions mark rough layering parallel to the foliation. Minor quartz occurs as isolated interstitial grains (Figs. 2d,3e). It has a weak undulatory extinction (Fig. 2d). Sample #431-2 has a higher amphibole content (76%) and a lower plagioclase content (13%) than sample #431-1. Amphibole has a weak SPO, marking a foliation (Fig. 3f). Both plagioclase and quartz occur as isolated grains in the amphibole matrix (Fig. 3f). Coarser plagioclase grains have a weak SPO parallel to the amphibole one.

Mineral composition

The major element compositions of amphibole, plagioclase, and garnet are reported in Table 2. The amphiboles are classified as pargasite, edenite, or ferropargasite. The plagioclases are oligoclase in 16lzwx2 and 16lzwx3 and andesine in 16lzwx1, #431-1, and #431-2, respectively. The garnets are aluminum-and iron-rich. The BSE images and composition profiles from core to rim within single grain show that amphibole and plagioclase grains have approximately homogeneous compositions, without compositional zoning (Figs. S2 andS3). The compositions of contacting amphibole and plagioclase grains are used for the estimation of equilibrium pressure and temperature conditions. The pressure is estimated with the Al content using an amphibole geothermometer [START_REF] Johnson | Experimental calibration of the aluminum-inhornblende geobarometer with application to Long Valley caldera (California) volcanic rocks[END_REF]). The temperature is estimated with an amphibole-plagioclase geobarometer [START_REF] Holland | Non-ideal interactions in calcic amphiboles and their bearing on amphibole-plagioclase thermometry[END_REF]. The estimated equilibrium pressure (~1 ± 0.05 GPa) and temperatures (696 ± 20 • C, 630 ± 18 • C, and 567 ± 11 • C) of the three EHS samples are similar to those proposed for the retrogression stage (M3) of the EHS garnet amphibolites [START_REF] Kang | Prolonged partial melting of garnet amphibolite from the Eastern Himalayan Syntaxis: Implications for the tectonic evolution of large hot orogens[END_REF]. The temperature of 16lzwx3 may be underestimated because the plagioclase was formed by retrogression of garnet. For the amphibole-rich xenoliths, the estimated equilibrium pressures are 1.05 ± 0.05 GPa in sample #521, 0.7 ± 0.05 GPa in sample #431-1, and 0.54 ± 0.05 GPa in sample #431-2, respectively. The equilibrium temperatures of #431-1 and #431-2 are 612 ± 9 • C and 587 ± 10 • C, respectively. The equilibrium temperature of sample #521 was not calculated due to the lack of plagioclase.

CPOs

The CPOs of the major constituent minerals were plotted as pole figures presenting orientations of the [100] axes, poles to the (010) Modal compositions were estimated as the area fraction occupied by each mineral on the EBSD maps, recalculated for a 100% indexation. Grain sizes (2D) were weighted (by the grain area) averages over the entire EBSD map. J-ppg was calculated using the mean orientation of each grain.

a Plagioclase grain sizes of 16lzwx1 and 16lzwx3 are not reliable due to large step sizes (15 μm) relative to fine grain size of plagioclase.

b Plagioclase J-ppg of 16lzwx3 is probably overestimated due to the minor content of plagioclase in this sample.

plane, and the [001] axes for amphibole and the [100] axes, poles to the (010) and ( 001) planes for plagioclase relative to the structural reference frame (lineation X and normal to the foliation Z) in Fig. 4. The CPOs of amphibole are clear in all samples. They are very strong in sample 16lzwx2 and weak in samples #521, 16lzwx1, and 16lzwx3; the remaining samples show amphibole CPOs of moderate intensity (Fig. 4). The fabric strength (J-index) of amphibole varies from 1.6 to 5.4. Three different CPOs symmetries are observed. The strong amphibole CPOs in sample 16lzwx2 is characterized by an L-type fabric: the [001] axes form a high concentration approximately parallel to lineation and the [100] and the [010] axes form girdles approximately normal to the lineation, with a high concentration of the [100] axes normal to the foliation. Amphibole in samples #431-1 and #431-2 has an S-type fabric: the [100] axes form a high concentration approximately normal to foliation and the [010] and the [001] axes form girdles approximately parallel to the foliation. Amphibole in other three samples show weak CPOs with an orthorhombic symmetry, that is an SL-type fabric: the [001] axes forms a maximum parallel to the lineation, with some dispersion in the foliation plane, and the [100] axes have a weak maximum normal to the foliation, with dispersion in the plane normal to the lineation. In samples #521 and 16lzwx3, a weak maximum of the [010] axes normal to the foliation is also observed.

The CPOs of plagioclase is much weaker than that of amphibole in most samples. It is close to random in samples 16lzwx1 and #431-2 (Table 3 andFig Intracrystalline deformation indicators (GOS map, GOS distribution histograms, and misorientation angle distributions) of the amphiboles in three selected samples (16lzwx2, #521, and #431-1) are presented in Fig. 5. The exhumed amphibolite 16lzwx2 has a non-negligible proportion of amphibole grains with high GOS values (grains colored in green to yellow; GOS > 5 • , in Fig. 5a), exhibits relatively large frequencies for small misorientation angles (18% of the correlated misorientations are <10 • , Fig. 5g). The uncorrelated misorientation angle distribution remarkably deviates from that of a random amphibole CPOs (Fig. 5g), which is consistent with the strong CPOs (Fig. 3b). In contrast, the xenoliths #521 and, especially, #431-1 are characterized by lower GOS values, indicating weaker intracrystalline deformation in amphibole (Fig. 5b,c,e,f). The uncorrelated misorientation angle distribution of sample #521 approximates that of a uniform amphibole CPOs (Fig. 5h), which is consistent with the weak CPOs (Fig. 3a). However, this sample shows a non-negligible proportion of low angle grain boundaries (~10% of the correlated misorientations are <10 • , Fig. 5h). Despite of the low intragranular deformation, the amphiboles in sample #431-1 show a clear CPOs (the uncorrelated misorientation angle distribution differs significantly from the one for a uniform CPOs, Fig. 5i). Both xenoliths display a peak at 180 • in the misorientation angle distribution, which indicates twinning, not visible optically, in amphibole.

Seismic properties

The calculated seismic properties are illustrated in Fig. 6. The seismic properties of 100% amphibole aggregates were calculated using the CPOs measured in each sample. To evaluate the effect of plagioclase on seismic properties of amphibolite, the seismic properties of samples 16lzwx2 and #431-1, which have high plagioclase contents, were also calculated with various amp/pl ratios (Fig. 7). Variation of the calculated P wave velocity (Vp) and the polarization anisotropy of S waves (AVs) as a function of the orientation of the foliation and lineation is displayed in Fig. 8. All seismic data are summarized in Table 3. P, pressure estimated with the Al content using the amphibole geothermometer of [START_REF] Johnson | Experimental calibration of the aluminum-inhornblende geobarometer with application to Long Valley caldera (California) volcanic rocks[END_REF]. T, temperature estimated with the amp-pl geobarometer of [START_REF] Holland | Non-ideal interactions in calcic amphiboles and their bearing on amphibole-plagioclase thermometry[END_REF]. n.d., not detected.

1 σ, standard deviation of multiple measurements (5-15 points for each mineral).

Mineral abbreviations: amp, amphibole; pl, plagioclase, grt, garnet; Alm, almandine; Grs, grossular; Prp, pyrope; Sps, spessartine; An, anorthite; Ab, albite.

Table 2

Major element compositions (wt%) of amphibole, plagioclase, and garnet and equilibrium conditions.

The seismic anisotropies vary significantly among the samples, consistent with the variation in strength of the amphibole CPOs. An increase in the plagioclase content decreases both the anisotropies and the velocities of the P-and S-waves (Table 3 and Fig. 6). If only the amphibole is considered, the calculated maximum P-wave anisotropies (AVp) range between 6.6% and 17.9%, the maximum S-wave polarization anisotropies (AVs) between 4.0% and 13.9%, the maximum P-wave velocities (Vp) between 7.3 and 7.8 km/s, the minimum Vp between 6.30 and 6.90 km/s, the Vp/Vs1 ratios between 1.62 and 1.84. Sample #431-1 and sample #521 have the strongest and the weakest seismic anisotropies, respectively.

The seismic property (Vp, AVs, and Vs1) patterns of amphibole are consistent with the amphibole CPOs patterns (Fig. 4), even for samples composed of 50% amphibole and 50% plagioclase (Fig. 6). Samples 16lzwx1-3 and #521 show a similar pattern, which has an orthorhombic symmetry (Fig. 6). The direction of maximum Vp is approximately parallel to the lineation (the [001] amp maximum) and the minimum Vp is observed for propagation normal to the foliation (the (100) amp maximum); the strongest AVs is observed for propagation directions approximately normal to the lineation and parallel to the foliation (Fig. 6a-e). Samples #431-1 and #431-2 show a different pattern, which has a dominantly hexagonal symmetry with a slow revolution axis; high Vp and strong AVs are observed for all propagation directions within or at a low angle to the foliation plane (Fig. 6g-j, Fig. 8a,b). Apparent isotropy or weak anisotropy is observed for S waves propagating normal to the foliation in samples with an S-type or an L-or LS-type amphibole CPOs, respectively. The polarization direction of the fast S-wave (Vs1) is parallel to the foliation for S waves propagating within or at a low angle to the foliation in all samples (Fig. 6). Fast polarization directions are parallel to the lineation for all other propagation directions in samples with an L-or LS-type amphibole CPOs (Fig. 6a-f). However, in samples with an S-type amphibole CPOs, the fast polarization directions are at a high angle to the lineation for propagation directions oblique to both the foliation and the lineation (Fig. 6g-j).

The seismic properties of samples 16lzwx2 and #431-1 with various amp/pl ratios are shown in Fig. 6b-d & g-i and Fig. 7. The results show that: 1) The AVp and the maximum AVs decrease progressively with the increase of plagioclase content, indicating that amphibole controls the Fig. 4. Crystallographic preferred orientations of amphibole (a-f) and plagioclase (g-k). Pole figures are plotted as lower hemisphere equal area projections. Color bars (multiples of uniform distribution units) are set as 0-11 for amphibole, 0-2.5 for plagioclase. The insert indicates the orientation of the structural reference frame (S: foliation; L: lineation). J-amp: J-index of amphibole; J-pl: J-index of plagioclase. For samples 431-1, and 431-2, where the lineation could not be identified in the hand sample, the CPOs were rotated to present the maximum [001] of amphibole parallel to the lineation of other samples. For sample 521, where neither the foliation nor the lineation was identified, the CPOs were rotated to present the [001] and (010) maxima of amphibole parallel to the lineation and foliation of other samples, respectively. seismic anisotropy. It is consistent with the stronger anisotropy of the amphiboles and the almost random CPOs of plagioclase in these rocks;

2) The maximum Vp and Vs 1 only slightly increase with the increase of amphibole content, suggesting the similar seismic velocities of amphibole and plagioclase.

Discussion

Mechanisms producing the amphibole CPOs

Amphibole CPOs may form by dislocation creep under hightemperature conditions. In this case, dominant activation of the [001] (100) slip system and secondary activation of the [100](010) and 1/2 Table 3 Calculated seismic properties of amphibolites using the measured CPOs and different amp/pl ratios.

〈110〉{110} systems, corroborated by the TEM data [START_REF] Skrotzki | Defect structure and deformation mechanisms in naturally deformed hornblende[END_REF], produce a strong orientation of the [001] axes parallel to the lineation and the (010) planes parallel to the foliation [START_REF] Díaz Aspiroz | Development of lattice preferred orientation in clinoamphiboles deformed under low-pressure metamorphic conditions. A SEM/EBSD study of metabasites from the Aracena metamorphic belt (SW Spain)[END_REF][START_REF] Elyaszadeh | Different slip systems controlling crystallographic preferred orientation and intracrystalline deformation of amphibole in mylonites from the Neyriz mantle diapir, Iran[END_REF][START_REF] Skrotzki | Defect structure and deformation mechanisms in naturally deformed hornblende[END_REF]. The CPOs of the studied samples could be the result of deformation by dislocation creep. The L or LStype fabrics are associated with simple or pure shear while the S-type fabrics are related to either axial compression or transpressive deformation [START_REF] Tommasi | Upper mantle tectonics: three-dimensional deformation, olivine crystallographic fabrics and seismic properties[END_REF]. Samples 16lzwx1 and 2 show significant intracrystalline misorientation in amphibole (Table 1 and Fig. 5a), consistent with deformation by a dislocation creep. Hence, the bimodal amphibole grain size distribution in sample 16lzwx1 (Fig. 3a) is likely the result of dynamic recrystallization. Recent experimental studies suggest rigid body rotation in fine-grained amphibolite could also contribute to the observed L-type amphibole CPOs [START_REF] Ko | Crystal preferred orientation of an amphibole experimentally deformed by simple shear[END_REF][START_REF] Kim | New crystal preferred orientation of amphibole experimentally found in simple shear[END_REF]. In samples #431-1 and #431-2, the low intracrystalline misorientation and well-developed growth twins indicate that amphibole grains were not significantly strained internally. Thus, dislocation creep cannot be the process responsible for the formation of the amphibole CPOs in these samples. Amphiboles in samples #521 and 16lzwx3 have GOS intermediate values between those of the two groups (Table1) and subhedral shapes (Fig. 2c), indicating that the microstructures and CPOs of these rocks record the activities of more than one process.

Observations of naturally and experimentally deformed amphibolerich rocks indicate indeed that a variety of mechanisms may produce strong CPOs of amphibole, including rigid body rotation by cataclastic flow under low-temperature and high strain rate conditions [START_REF] Ko | Crystal preferred orientation of an amphibole experimentally deformed by simple shear[END_REF][START_REF] Kim | New crystal preferred orientation of amphibole experimentally found in simple shear[END_REF], diffusion creep in fine-grained aggregates [START_REF] Getsinger | Amphibole fabric formation during diffusion creep and the rheology of shear zones[END_REF], dissolution-precipitation creep when deformation occurs under melt/fluid-rich conditions (Giuntoli et al., 2018;[START_REF] Stokes | Deformation of amphibolites via dissolution-precipitation creep in the middle and lower crust[END_REF], rotation of non-equant amphibole grains in a deforming magma chamber [START_REF] Berger | Deformation mechanisms and reaction of hornblende: examples from the Bergell tonalite (Central Alps)[END_REF], oriented growth during partial crystallization of low melt fractions in a deforming rock [START_REF] Tommasi | Hydrous melts weaken the mantle, crystallization of pargasite and phlogopite does not: Insights from a petrostructural study of the Finero peridotites, southern Alps[END_REF], and topotactic growth during hydration reactions [START_REF] Padrón-Navarta | Oriented growth of garnet by topotactic reactions and epitaxy in high-pressure, mafic garnet granulite formed by dehydration melting of metastable hornblende-gabbronorite (Jijal Complex, Kohistan Complex, north Pakistan)[END_REF]. All these processes, except topotactic growth, have a common characteristic: they produce correlated amphibole CPOs and SPO (though the SPO may be erased by later textural re-equilibration and does not necessarily develop in the case of metamorphic reactions), but do not involve intracrystalline deformation. However, in magma chambers and fluid-rich settings, compaction may result in the superimposition of small amounts of solidstate deformation in the already crystallized parts of the system, leading to moderate intragranular deformation. Processes involving amphibole growth from fluids or melts may result in compositional zoning if the fluid composition changes during the growth [START_REF] Giuntoli | Replacement reactions and deformation by dissolution and precipitation processes in amphibolites[END_REF], but not if the fluid composition remains constant. Given the fact that anisometric growth of amphibole results in prismatic grains elongated parallel to the [001] and flattened parallel to the normal to (100), the CPOs produced by processes depending on the anisotropic grain growth of amphibole bears the same relation to the flow direction and plane as that produced by dislocation creep.

The textures of samples #431-1, #431-2, and #521 suggest that magmatic flow is likely the dominant mechanism for forming the amphibole CPOs. This inference is consistent with the petrological and geochemical data, which indicated that these amphibole-rich xenoliths represent residual cumulates of fractional crystallization of basaltic magma [START_REF] Xu | Constructing the early Mesozoic Gangdese crust in southern Tibet by hornblende-dominated magmatic differentiation[END_REF]. In a deforming magma mush, the symmetry of the CPOs will depend on the deformation regime, similarly to CPOs evolution by dislocation creep, but it also depends on the shape of the crystals [START_REF] Ildefonse | Effect of mechanical interactions on the development of shape preferred orientations: a two-dimensional experimental approach[END_REF]. The L-or LS-type fabrics will develop in simple or pure shear (2D flow) for prismatic grains, whereas the S-type fabric will form in axial compressions (such as compaction and transpression) and in pure or simple shear of platy grains. Development of the S-type fabrics as observed in samples #431-1 and #431-2 (Fig. 4e,f) may result from either axial compression, a strong component of compression normal to the shear plane, or pure or simple shear affecting platy hornblende grains with a weak elongation parallel to the [001]. The higher GOS values of amphibole in sample #521 may result from small amounts of solid-state deformation in the already crystallized parts of the system, maybe in response to compaction. Variation in the intensity and symmetry of the amphibole CPOs among the three samples may result from the sampling of different parts of the magma chamber. This interpretation is consistent with the variations in modal composition. The weak CPOs of plagioclase in these samples are coherent with the interstitial shapes of this phase, which implies late crystallization.

Despite the clear SPOs displayed by all three 16lzwx samples, the variation in the intensity of intragranular misorientation in amphibole among the amphibolites from the NBC complex might be explained by added contributions of oriented growth, dissolution-crystallization creep and dislocation creep. Given that these rocks are mafic migmatites [START_REF] Kang | Prolonged partial melting of garnet amphibolite from the Eastern Himalayan Syntaxis: Implications for the tectonic evolution of large hot orogens[END_REF], it is likely that the amphibole CPOs started to develop during the high-pressure granulite facies peak metamorphism (partial melting) and continued to evolve during the retrogression (exhumation). Note that sample 16lzwx2 show the strongest CPOs and SPO and intermediate GOS values in amphibole due to no garnet porphyroblast perturbation to the orientation of the foliation. When present as a major constituent of the rock, as in 16lzwx2, plagioclase shows a weak CPOs, which is consistent with the amphibole CPOs, pointing to a coherent deformation. In samples 16lzwx1 and 16lzwx3, plagioclase is a minor phase mainly occurring around garnet and locally forming pressure shadows (Fig. 3c). It may explain the poor orientation of plagioclase and the obliquity between the plagioclase CPOs and the foliation in sample 16lzwx3.

Implication for crustal seismic anisotropy in southern Tibet

The PmS delay time of the crust is up to 0.8 s with typical values of 0.45-0.75 s in southern Tibet [START_REF] Kong | Azimuthal anisotropy and mantle flow underneath the southeastern Tibetan Plateau and northern Indochina Peninsula revealed by shear wave splitting analyses[END_REF][START_REF] Sun | Crustal structure and deformation under the Longmenshan and its surroundings revealed by receiver function data[END_REF][START_REF] Sun | Crustal structure and deformation of the SE Tibetan plateau revealed by receiver function data[END_REF][START_REF] Wu | Deformation of crust and upper mantle in Central Tibet caused by the northward subduction and slab tearing of the Indian lithosphere: New evidence based on shear wave splitting measurements[END_REF][START_REF] Wu | Complex seismic anisotropy beneath western Tibet and its geodynamic implications[END_REF][START_REF] Zheng | Crustal azimuthal anisotropy beneath the southeastern Tibetan Plateau and its geodynamic implications[END_REF]. The average delay time of shear wave splitting associated with aligned crackinduced anisotropy in the brittle upper crust (~15 km) is generally <0.15 s [START_REF] Crampin | The fracture criticality of crustal rocks[END_REF]. Thus, 2/3 of the Pms delay times measured in southern Tibet are likely generated in the mid-to-lower crust. The delay time is sensitive to both the thickness of the anisotropic layer and the strength of the anisotropy along the propagation direction. The strength of the anisotropy depends on the orientation of the CPOs relative to the propagation direction, which is within 10 • from vertical for the PmS wave (Fig. 8). We have calculated the delay times associated with rocks composed of 100% amphibole or 50% amphibole and 50% plagioclase based on the CPOs data of amphibole and plagioclase in samples 16lzwx2 and #431-2 using the following equation [START_REF] Silver | Seismic anisotropy beneath the continents: probing the depths of[END_REF]:

δt = ( AV S / V s ) × L (1)
where δt is the shear wave delay time, AV S is the anisotropy for the considered propagation direction, V s is the shear wave average velocity, and L is the thickness of the anisotropic layer. The maximum anisotropy will be observed if the anisotropic layer in the deep crust has a vertical foliation and a horizontal lineation, that is a strike-slip fabric (dip angle = 90 • in Fig. 8f), or a vertical planar fabric (the dip angles for sample #431-1 in Fig. 8b or the dip angle = 90 • in Fig. 8d). In this case, 2/3 of the PmS delay time observed in southern Tibet may be explained by an amphibole-rich (~50% in vol.) layer with a thickness of 14-22 km (Fig. 9a). Significant delay times are nevertheless produced for all propagation directions at a high angle to the lineation in the case of the LS-type fabrics (Fig. 8b). Propagation directions at a high angle to the foliation result in low delay times in the case of the S-type and the LStype fabrics (a dip angle of 0 • in Fig. 8d or f). The observed delay times cannot be explained by amphibole CPOs associated with flat-lying foliations. Low delay times will also be produced for dipping lineation (the dip angle of >30 • for sample 16lzwx2 in Fig. 8d). Seismic radial anisotropy, which describes the discrepancies between the horizontally (Vs h ) and vertically (Vs v ) polarized shear wave velocities, is typically inferred from observations of a period-dependent discrepancy between the phase or group speeds of Rayleigh and Love waves [START_REF] Ritzwoller | Global surface wave diffraction tomorgraphy[END_REF][START_REF] Shapiro | Thinning and flow of Tibetan crust constrained by seismic anisotropy[END_REF]. Strong positive radial anisotropy (Vs h > Vs v , 4.8 ± 1.4% at 25-35 km depth) was observed in southeastern Tibet below the study area (Fig. 1b), which was interpreted as the result of lateral flow of the Tibetan crust [START_REF] Agius | Complex, multilayered azimuthal anisotropy beneath Tibet: evidence for co-existing channel flow and pure-shear crustal thickening[END_REF][START_REF] Hu | Direct surface wave radial anisotropy tomography in the crust of the Eastern Himalayan Syntaxis[END_REF][START_REF] Shapiro | Thinning and flow of Tibetan crust constrained by seismic anisotropy[END_REF][START_REF] Xie | Crustal radial anisotropy across Eastern Tibet and the Western Yangtze Craton[END_REF]. Negative radial anisotropy is, in contrast, observed north of the study zone [START_REF] Hu | Direct surface wave radial anisotropy tomography in the crust of the Eastern Himalayan Syntaxis[END_REF]. The S-type amphibole CPOs can produce strong positive radial anisotropy (9% for amphibolite of 50%amp-50%pl, Fig. 6i), independent of the propagation direction, when the foliation is subhorizontal. In the case of an LS-type amphibole CPOs, the radial anisotropy would also be positive, but its intensity will depend on the propagation direction (back azimuth), being minimum for propagation parallel to the lineation (Fig. 6d). Vertical foliations will, for both the LSand S-type amphibole CPOs, produce negative radial anisotropy, which will strongly vary as a function of the propagation direction. However, in the case of an LS-type CPOs, if the foliation is vertical and the lineation is horizontal, the radial anisotropy will nevertheless remain low for all back azimuths (Fig. 6d).

In conclusion, the strong crustal anisotropy observed in southern Tibet north of the study zone, using both PmS and surface waves, might be explained by a 15-22 km thick layer composed of ≥50% vol. of amphibole with an S-type CPOs and a steeply dipping foliation, as in a transpressive shear zone or a synkinematic magmatic body intruding a strike-slip zone. Such a fabric will produce fast PmS polarization directions parallel to the trend of the foliation. This is consistent with the PmS polarizations mainly parallel to the trend of the multiple dextral strike-slip faults accommodating the lateral escape of southeast Asia east of the study zone (Fig. 1). However, it is difficult to produce within a single deformed layer both a strong PmS anisotropy and a strong positive radial anisotropy, as the latter requires flat-lying foliations and the polarization anisotropy detected by S-waves propagating normal to the foliation is weak (for the LS-type CPOs, Fig. 6d) to null (for the S-type CPOs, Fig. 8k). A dip-angle variation of the foliation with depth is thus required.

Finally, amphibolites are common rocks, but the volumes needed to explain the observed delay times are significantly larger than those usually observed in the exposed middle to lower crust sections. Locally high Vp/Vs ratios (>1.8) in the crust of southern Tibet suggest the presence of mafic rocks, but lower values, typical of felsic rocks, predominate [START_REF] Hu | Direct surface wave radial anisotropy tomography in the crust of the Eastern Himalayan Syntaxis[END_REF]. As a consequence, CPOs of other minerals, in particular micas in the shallower sections of the crust, probably also contributes to the observed anisotropy.

Conclusion

In this study, the microstructure and seismic properties of amphibole-rich rocks from southern Tibet are presented. These samples display clear amphibole CPOs characterized by pronounced concentrations of the [001] axes parallel to the lineation and the (100) planes parallel to the foliation, respectively. However, the intensity and symmetry of the amphibole CPOs vary between samples, with the L-and LStype fabrics predominating among the samples from the NBC complex and the S-type fabrics among the xenoliths of the Gangdese batholith. Analyses of the microstructures and geochemical data suggest that the amphibole CPOs in the xenoliths formed in response to magmatic flow, whereas those in the NBC amphibolites likely formed by an association of oriented crystallization, dislocation creep, and dissolutionprecipitation creep.

All samples are characterized by strong seismic anisotropies, which tend to decrease with increasing plagioclase content. An amphibole-rich layer (15-22 km) with a steeply-dipping foliation and an S-type amphibole CPOs might account for both the strong azimuthal PmS anisotropy and the negative radial anisotropy observed to the north and east of the study zone. However, the required amphibolite volumes are too large, indicating that other minerals, in particular micas, probably also contributed to the observed anisotropy. An amphibole-rich layer with horizontal foliation at depths of 35 km may, on the other hand, explain the strong positive radial anisotropy in southern Tibet. These observations are best explained by the S-type amphibole CPOs, like those observed in the amphibolite xenoliths, which we interpret as produced by the magmatic flow.

The present observations corroborate that strong amphibole CPOs and crustal anisotropy may be formed by a variety of processes, ranging from magmatic deformation to solid-state creep. These different processes produce a similar relation between amphibole CPOs and flow direction and plane. However, in the case of magmatic flow, the S-type fabric may form even in simple shear if the amphibole grains are more flattened than elongated.

Fig. 1 .

 1 Fig. 1. a) Geomorphic map of the Tibetan-Himalayan orogenic belt. White bars represent crustal shear wave azimuthal anisotropies (fast wave polarization directions and delay times) measured using the PmS phase. Data sources of the PmS anisotropy: southwestern Tibet[START_REF] Wu | Complex seismic anisotropy beneath western Tibet and its geodynamic implications[END_REF], central south Tibet[START_REF] Wu | Deformation of crust and upper mantle in Central Tibet caused by the northward subduction and slab tearing of the Indian lithosphere: New evidence based on shear wave splitting measurements[END_REF], and east Tibet[START_REF] Kong | Azimuthal anisotropy and mantle flow underneath the southeastern Tibetan Plateau and northern Indochina Peninsula revealed by shear wave splitting analyses[END_REF][START_REF] Sun | Crustal structure and deformation of the SE Tibetan plateau revealed by receiver function data[END_REF][START_REF] Sun | Crustal structure and deformation under the Longmenshan and its surroundings revealed by receiver function data[END_REF][START_REF] Zheng | Crustal azimuthal anisotropy beneath the southeastern Tibetan Plateau and its geodynamic implications[END_REF]. The color bar at the bottom is elevation in kilometers. The white box denotes the location of Figure1b. b) Radial anisotropy in southeastern Tibet at 35 km depth (modified from Hu et al., 2020). The color bar at the bottom shows the range of radial anisotropy. Red and blue colors refer to positive and negative radial anisotropies, respectively. The black box marks approximately the location of the Eastern Himalayan Syntaxis (EHS). Stars indicate locations of studied samples: 16lzwx1-3 are the samples collected in the NBC complex while #521, #431-1, and #431-2 are collections from the xenoliths. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .

 2 Fig. 2. Photomicrographs illustrating typical microstructures of amphibole-rich rocks from southern Tibet. a) Nearly round-shaped garnet (Grt) porphyroblasts enclosed in a matrix mainly composed of amphibole (Hbl) and minor fine-grained plagioclase (Pl) and quartz (Qtz). b) Distinct foliation defined by the amphibole SPO and the shapes of plagioclase aggregates. c) and d) Details of the microstructure illustrating the lack of evidence of intracrystalline deformation in both amphibole and plagioclase in sample 161zwx3 and the xenoliths. Minor quartz occurs as isolated irregularly shaped grains. In (d) plagioclase has polysynthetic twins and quartz has a weak undulatory extinction. (a-c) are plane-polarized light images and (d) is a cross-polarized light image. Ilm, ilmenite. Scale bars: a) 500 μm, b) 1000 μm, c) 300 μm, d) 400 μm.

Fig. 3 .

 3 Fig. 3. EBSD phase maps of amphibole-rich rocks. Scale bars: 2 mm in (a, b, c) and 5 mm in (d, e, f).

  . 4). Sample 16lzwx2 has the best-defined plagioclase CPOs characterized by alignment of the [100] maximum subparallel to the lineation and of the [001] maximum normal to the foliation. All samples show a poor concentration of the [001] axes at a high angle to the foliation and of the [100] axes at a low angle to the lineation, respectively. In sample 16lzwx3, the plagioclase CPOs is oblique (~20-25 • ) to the macroscopic lineation and foliation. In samples 16lzwx1 and #431-2, the distributions of the [001] and [010] axes are similar while the [100] axes are dispersed in a wide girdle at a low angle to the foliation.

Fig. 5 .

 5 Fig. 5. Amphibole GOS maps (a-c), distributions of GOS (d-f), and misorientation angle distributions (g-i) from three selected samples (16lzwx2, #521, and #431-1). The black lines in GOS maps are high angle (>10 • ) boundaries. The colored scale bar at the left represents GOS values of 0 • to 10 • from bottom to top. Misorientation angle distributions are presented as neighboring points (correlated), randomly selected points (uncorrelated), and the random CPOs (uniform) of amphibole is also shown for references. The misorientation peaks at 180 • in samples #521 and #431-1 likely indicate twinning in amphibole. Scale bars: 2 mm (a), 5 mm (b, c).

Fig. 6 .

 6 Fig. 6. Calculated seismic properties of 100% amphibole aggregates for all samples and amp-pl aggregates with various amp/pl ratios for sample 16lzwx2 (b-d) and #431-1 (g-i). Vp: P wave velocity; AVs: S wave anisotropy; Vs1: fast shear wave velocity; AVp: P wave anisotropy; Max. AVs: maximum S wave anisotropy; Max. Vs1: maximum fast shear wave velocity. AVs = 200*(Vs1 -Vs2)/(Vs1 + Vs2); AVp = 200 * (Max.Vp -Min.Vp)/(Max.Vp + Min.Vp). The black bars in AVs and Vs1 plots indicate the fast wave polarization directions. The reference frame is the same as in Fig. 4.

Fig. 7 .

 7 Fig. 7. Relationship between mineral proportion and seismic properties using the CPOs of amphibole and plagioclase in the samples 16lzwx2 and #431-1.

Fig. 8 .

 8 Fig. 8. Variation of calculated P wave velocity and AVs as a function of the dip angle of the XY plane (a, b), XZ plane (c, d), and YZ plane (e, f). X marks the lineation and Z the normal to the foliation.

Fig. 9 .

 9 Fig. 9. Relationship between delay time of S wave and thickness of anisotropic layers of 100% amphibole aggregates and amphibolite composed of 50% amp and 50% plagioclase. a) The delay times calculated using the maximum AVs value. b), c), d) The delay times calculated using the AVs values for propagation directions of the PmS parallel to X, Y, and Z, respectively. The shaded area shows the range of delay times observed from the middle-lower crust in southern Tibet. It is assumed that 2/3 of the delay times are contributed by the middle-lower crust.

Table 1

 1 Modal composition and microstructure parameters derived from EBSD mapping.

	Sample	Modal composition (%)	Index rate	Amp

  AVs1, maximum propagation anisotropy for the fast S wave; AVs2, maximum propagation anisotropy for the slow S wave; Max. Vp, Max. Vs1, Max. Vs2, maximum P wave velocities, fast S wave, and slow S wave; Min. Vp, Min. Vs1, Min. Vs2, minimum P wave velocities, fast S wave, and slow S wave.

	Sample	Composition	AVp	Max. Vp	Min.	Max.	Max.	Max.	Min.	Min.	Max.	Min.	Max. Vp/	Min. Vp/
			(%)	(km/s)	Vp	AVs	AVs1	Vs1	Vs1	AVs2	Vs2	Vs2	Vs1
					(km/s)	(%)	(%)	(km/s)	(km/s)	(%)	(km/s)	(km/s)
	16lzwx1	100%amp	12.39	7.58	6.69	8.26	5.97	4.29	4.04	5.34	4.15	3.93	1.79	1.64
	16lzwx2	100%amp	15.96	7.79	6.64	9.79	6.53	4.34	4.07	7.54	4.17	3.87	1.84	1.62
		75%amp + 25%	13.71	7.57	6.60	9.35	6.20	4.22	3.97	6.88	4.08	3.81	1.83	1.64
		pl
		50%amp + 50%	11.04	7.32	6.55	8.84	5.89	4.10	3.86	6.24	3.97	3.73	1.82	1.66
		pl
	16lzwx3	100%amp	8.45	7.40	6.80	5.96	4.25	4.24	4.06	3.51	4.12	3.98	1.76	1.66
	#521	100%amp	6.58	7.34	6.90	4.04	3.25	4.19	4.05	3.08	4.13	4.01	1.77	1.65
	#431-1	100%amp	17.89	7.59	6.30	13.88	13.56	4.45	3.89	3.41	4.01	3.87	1.71	1.63
		75%amp + 25%	14.57	7.45	6.40	11.50	11.39	4.33	3.86	2.62	3.96	3.86	1.72	1.66
		pl
		50%amp + 50%	10.87	7.29	6.50	9.03	8.83	4.19	3.84	1.80	3.90	3.83	1.75	1.70
		pl
	#431-2	100%amp	17.21	7.57	6.40	13.11	12.93	4.44	3.90	3.03	4.01	3.89	1.71	1.63

Note: AVp, maximum P wave propagation anisotropy; AVs, maximum S wave polarization anisotropy;

Acknowledgments

We dedicate this article to Prof. Adolphe Nicolas, an exceptional researcher who contribute to our understanding of how deformation occurs in the Earth for over 50 years. We are grateful to Fabrice Barou and Wenlong Liu for guidance on EBSD measurement. Deng Xiao is appreciated for help on microprobe analysis. Huizi Jian is thanked for processing seismic data. Pengyuan Han is deeply thanked for sample (#521, #431-1, and #431-2) collection from southern Tibet. Constructive reviews of Dr. Haemyeong Jung, Dr. Loic Labrousse, and editor Dr. Philippe Agard helped us to improve the paper. This work was completed while Xiong Wang was a PhD student at China University of Geosciences (Wuhan) and Université de Montpellier, supported by the National Natural Science Foundation of China (Grants 41425012 and 41590623) and the MOST special fund from the State Key Laboratory of GPMR. This paper was written while Xiong Wang was a postdoctoral fellow at Southern University of Science and Technology, funded by the National Natural Science Foundation of China (Grant 41974098).

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi. org /10.1016/j.tecto.2021.228869.


Mauro Zucchelli

Samuel Deslauriers-Gauthier

Rachid Deriche

Brain Tissue Microstructure Characterization Using dMRI Based Autoencoder Neural-Networks

Keywords: Spherical Harmonics, Rotation Invariant Features, Neural Networks, Diffusion MRI, Autoencoder

In recent years, multi-compartmental models have been widely used to try to characterize brain tissue microstructure from Diffusion Magnetic Resonance Imaging (dMRI) data. One of the main drawbacks of this approach is that the number of microstructural features needs to be decided a priori and it is embedded in the model definition. However, the number of microstructural features which is possible to obtain from dMRI data given the acquisition scheme is still not clear. In this work, we aim at characterizing brain tissue using autoencoder neural networks in combination with rotation-invariant features. By changing the number of neurons in the autoencoder latent-space, we can effectively control the number of microstructural features that we obtained from the data. By plotting the autoencoder reconstruction error to the number of features we were able to find the optimal trade-off between data fidelity and the number of microstructural features. Our results show how this number is impacted by the number of shells and the bvalues used to sample the dMRI signal. We also show how our technique paves the way to a richer characterization of the brain tissue microstructure in-vivo.

Introduction

Diffusion Magnetic Resonance Imaging (dMRI) has been widely used in research to probe the brain tissue microstructure in the living human brain. To estimate microstructural features from the dMRI signal, it is first necessary to select a biophysical model which can capture the desired features [START_REF] Novikov | Quantifying brain microstructure with diffusion mri: Theory and parameter estimation[END_REF]. The choice of the model is limited by the type of dMRI data available, e.g. the type of diffusion sequence, the number of gradients (or shells), the number of echo times, etc [START_REF] Lampinen | Neurite density imaging versus imaging of microscopic anisotropy in diffusion mri: A model comparison using spherical tensor encoding[END_REF]. Although this approach has been successful, it has several limitations. The biggest limitation is that the choice of the model completely decides the number and type of microstructural feature which are obtainable from the data [START_REF] Novikov | On modeling[END_REF]. A model-fitting will always provide a result but the accuracy of this result is not simple to verify. Nowadays, several software tools [START_REF] Fick | The dmipy toolbox: Diffusion mri multicompartment modeling and microstructure recovery made easy[END_REF][START_REF] Garyfallidis | Dipy, a library for the analysis of diffusion mri data[END_REF] easily permit researchers to fit multiple models on the same data. However, the comparison and validation of such models is still an open question 1 .

Neural networks (NN) have become widely popular in medical imaging for their astonishing classification given a sufficiently large and accurate training set [START_REF] Liu | A comparison of deep learning performance against healthcare professionals in detecting diseases from medical imaging: a systematic review and meta-analysis[END_REF]. One of the most important features of NN is that they can be considered a model-free approach. Normally, no information on the characteristics of the data that we want to classify is given to the NN, because they are able to learn the classes representation directly from the training data.

NN have been already applied to dMRI data in order to obtain microstructural features [START_REF] Golkov | q-space deep learning: Twelve-fold shorter and model-free diffusion mri scans[END_REF][START_REF] Zucchelli | Investigating the effect of dmri signal representation on fully-connected neural networks brain tissue microstructure estimation[END_REF]. In these works, a multi-compartment model has been chosen for generating the microstructural features used as labels for the NN training set. This approach limits the representation power of the NN to the model used for labeling the training data, not fully exploiting the model-free characteristic of the NN.

Autoencoders are a particular family of NN, principally used for denoising and feature extraction [START_REF] Vasilev | q-space novelty detection with variational autoencoders[END_REF]. An autoencoder NN is composed of two parts, an encoder, and a decoder. The encoder that takes as input the data and outputs a certain number of features (also known as latent-space), corresponding to the number of neurons in the last layer of the encoder. The decoder takes as input these features and outputs the original data (see Fig. 1). An autoencoder NN is trained by trying to minimize the error between the data and the output of the decoder. In order to prevent the network from learning the identity operation, the number of features produced by the encoder is strictly lower than the size of the input data. This enables the autoencoder NN to codify each input in a small number of features which can then be used for classification or other tasks. Autoencoders have been recently used in dMRI by Vasilev and colleagues to detect abnormal voxels in multiple-sclerosis patients' data [START_REF] Vasilev | q-space novelty detection with variational autoencoders[END_REF]. One of the main limitations of using autoencoders with dMRI data is that the dMRI signal is strongly directional-dependent and even a slight change in the orientation in the underlying brain fiber bundles will result in a completely different signal profile. This means that if the dMRI signal samples are used for training the autoencoder at least some of its features will be used to capture the orientation profile. Tissue orientation can be easily retrieved from dMRI data using the popular spherical deconvolution approaches [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion mri: non-negativity constrained super-resolved spherical deconvolution[END_REF], therefore it is generally not included in the scope of microstructural estimation technique. In fact, Rotation Invariant Features (RIF) can be used as an input microstructural models, instead of the row signal samples [START_REF] Novikov | Rotationallyinvariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion mri[END_REF][START_REF] Zucchelli | A computational framework for generating rotation invariant features and its application in diffusion mri[END_REF].

In this article, we aim to combine high order RIF, derived from the 4 th order Spherical Harmonics expansion of the dMRI signal, with autoencoders. Our purpose is to characterize the number of microstructural features which is possible to obtain from dMRI pulse-gradient spin-echo (PGSE) data, regardless of the tissue orientation. We also investigate how this number changes to the number of shells used to sample the data. Our results show that this technique is indeed feasible to model dMRI data and show also that increasing the number of shells and the b-value provide a richer characterization of the brain tissue microstructure in-vivo.

Materials and methods

Rotation-invariant features of the dMRI signal

The dMRI signal at each b-value can be represented as a spherical function f (u), where the unit vector u corresponds to the gradient direction. Considering the diffusion signal is real-valued and antipodal-symmetric, we can expand it into a truncated Spherical Harmonics (SH) series

f (u) = L l=0,even l m=-l c lm Y m l (u) (1)
where c lm are the real-valued SH coefficients Y m l are the real SH [START_REF] Homeier | Some properties of the coupling coefficients of real spherical harmonics and their relation to gaunt coefficients[END_REF] of degree l and order m, and L is the maximum degree. In [START_REF] Zucchelli | A computational framework for generating rotation invariant features and its application in diffusion mri[END_REF], the authors defined a set of algebraic independent RIF. Each invariant can be calculated as

I l [f ] = l1 m1=-l1 • • • l d m d =-l d c l1m1 • • •c l d m d G(l 1 , m 1 | • • • |l d , m d ) (2)
where G are the generalized Gaunt coefficients [START_REF] Homeier | Some properties of the coupling coefficients of real spherical harmonics and their relation to gaunt coefficients[END_REF], namely the integral of d SH and the vector l

is equal to [l 1 , l 2 • • • , l d ].
Eq. ( 2) provides a rotation invariant indices for any set of l. In this work, we focused on the 12 algebraic independent RIF obtained from the SH expansion at L = 4 of the diffusion signal [START_REF] Zucchelli | A computational framework for generating rotation invariant features and its application in diffusion mri[END_REF].

Being invariants, these RIF depends only on the "shape" of the signal and not on the orientation. For example, the value of the RIF change between crossing fibers and single fiber voxels, if two voxels present different intra-axonal signal fractions, or if the water diffusivity changes between the voxels [START_REF] Novikov | Rotationallyinvariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion mri[END_REF]. Since multiple factors contribute to the invariants numerical values they cannot be considered biomarkers per se, but by removing the orientation-dependency they greatly improve the estimation of brain tissue microstructure. In particular, neural network-based microstructural feature estimation considerably benefits from this form of rotation invariant signal representation [START_REF] Zucchelli | Investigating the effect of dmri signal representation on fully-connected neural networks brain tissue microstructure estimation[END_REF].

In-vivo dataset

In this work, we consider three subjects of Human Connectome Project (HCP) diffusion MRI young-adults dataset [START_REF] Sotiropoulos | Advances in diffusion MRI acquisition and processing in the human connectome project[END_REF]. For each subject, 288 dMRI volumes has been acquired: 18 volumes at b-value 0 s/mm

Autoencoder design

Figure 1 shows a graphical representation of the autoencoder used in this work. We design our autoencoder as a fully connected NN for taking as input the 12 RIF calculated from each shell of the HCP dataset. Therefore, the number of input (and output) channels corresponds to 12 for the single-shell b=1000 s/mm2 data, 24 for the two shells up to b=2000 s/mm 2 , and 36 for the complete dataset with three shells up to b=3000 s/mm 2 . Following the input layer, we add a first hidden layer composed of 50 neurons combined with a Rectified Linear Unit (ReLU) functions. After this layer, the data is passed to another layer of neurons plus ReLU forming the latent-space layer. The output of the latentspace layer represents the set of microstructural features that characterize the dMRI signal. The cardinality of the latent-space layer is the main variable of interest for our autoencoder design. This number defines the amount of features that the autoencoder uses to represent the input. In this work, we created several neural networks with latent-space ranging from 3 neurons up to the number of the input channels. In this last case, the network should be able to simulate the identity transform perfectly recovering the input, therefore we use it only as a reference limit-case. These first three elements form the encoder. The decoder takes as input features output of the latent space and passes it to the second hidden layer composed of 50 neurons. As in the case of the first hidden layer, each neuron is followed by a ReLU function. The output of the second hidden layer is forwarded to the neurons of the output layer having the same number as the input of the network. Each autoencoder is optimized trying to minimize the Mean Square Error (MSE) between the input RIF at the different shell and the output of the network. Since each RIF have a different range of values, we scale each invariant according to a factor of 1, 10, or 100 to bring them in the same numerical range. This normalization is necessary in order to avoid that the biggest invariant (I 0 , the mean of the signal ) dominates the optimization reducing the influence of the other RIF. We trained the network using each of the three HCP subjects, independently, splitting the voxel randomly. Sixty percent of each brain voxels were used for training the network, and forty percent for testing. To test the robustness of the method, we used 10 different random seeds for dividing the data into training and testing, and for the network initialization and optimization. All the autoencoders implemented in this work are written in python using the Pytorch 2 library and trained using the Adam optimizer with a learning rate of 5 × 10 -5 and a weight decay of 1 × 10 -6 . The training set is split in batches of 1000 elements and we trained the networks for 300 epochs.

All of the autoencoder hyper-parameters (including the number of layers, the learning rate, and the weight decay), were optimized by minimizing the MSE of the invariants on one of the HCP subjects.

Results and Discussion

In this work, we used the RIF extracted from the different shells of the dMRI signal as input, and the autoencoder latent-space represents the microstructural features. However, it is unclear how many microstructural features is it possible to extract from PGSE multi-shell dMRI data, in particular, as a function of the number of shells and the b-values used in the acquisitions. Figure 2 shows the log MSE between the original and reconstructed RIF for the testing set voxels of three HCP subjects (in red, green, and blue) as a function of the number of neurons in the latent-space. Each transparent line represents one random-seed instance of the same network, and the solid line the average of the ten networks results. Figure 2 (left) shows the results for the autoencoders trained using the 12 RIF at b-value 1000 s/mm 2 . At this b-value, the difference in the reconstruction error between the subjects is in the same range of the changes obtained training the same subject with a different network initialization, in particular, when the number of features of the latent-space is lower than 9. We calculated the elbowpoint of the average of the curves for each subject, and we mark it on the graph with black crosses. The elbow-point represents an empirical way to find the "optimal" trade-off between reconstruction error and the number of features. In our experiments, two subjects (160123 and 473952) have reached an optimal number of features with 6 neurons in the latent-space, and one subject (174437) with 9 neurons. These features contain information on both the inner properties of the axons and extracellular space (like the axonal density and water diffusivity) and information regarding the axons architecture (crossing bundles, fibers fanning, etc). This number of features is higher than the information usually retrieved from single-shell data, which are generally modeled as the three eigenvalues of the diffusion tensor. However, in this case, we used SH of order 4 for modeling the signal, while the diffusion tensor can be considered an order 2 representation. Therefore, the increase in the number of features can be explained by the increase in the order of the signal representation. Figure 2 (center) shows the same graph for the two shells data at b-value 1000 and 2000 s/mm 2 . In this case, we have 24 RIF as input (12 for each shell). The autoencoders show a different reconstruction error for each subject when the latent-space grows bigger than 15 neurons. Before this point, the difference between the MSE of the different subjects is in the range of variance induced by the change of the random seed.

The optimal number of features, calculated by finding the elbow-points of the curves, increase to 13 for subjects 160123 and 473952 and 14 for subject 174437. This means that there is a minimum increase of 5 features while moving from single-shell low-b-value data, to two shells data. Training the network with all the three shells, up to b-value 3000 s/mm 2 further increase this difference (see Figure 2 right). In this case, we have 36 input channels considering 12 RIF per shell. Using the b-value 3000 s/mm 2 data we observe a clear difference between the reconstruction error of the three subjects. As in the previous cases, this difference increases with the increase of the number of features of the latent-space.

It is necessary to underline that each subject has been used for training its autoencoder independently to the other two subjects. This means that there is something in the b-value 3000 s/mm 2 RIF which led the autoencoder to have a systematic different reconstruction error between the subject. This can be due to noise or another hidden acquisition artifact that may be completely uncorrelated to tissue microstructure. Although the errors are different, the trend of the errors is very close to each other. In fact, all three subjects have an optimal number of features at 15 latent-space neurons. The increase to the two-shells data is reduced to only one additional microstructural feature. This means that adding a third shell at b-value 3000 s/mm 2 while adding some subject-specificity and reduce the variance of the results does not increase the number of microstructural features obtainable from the data sensitively. Although we do not have the data, we can imagine that with a further increase in the number of shells we will probably reach a plateau in the number of microstructural features. This result is in line with other studies whose highlight that in order to be sensitive to the full brain tissue microstructure it is necessary to include other dMRI acquisition data, for example, Spherical Encoding [START_REF] Lampinen | Neurite density imaging versus imaging of microscopic anisotropy in diffusion mri: A model comparison using spherical tensor encoding[END_REF], shells acquired using multiple diffusion-times [START_REF] Novikov | Quantifying brain microstructure with diffusion mri: Theory and parameter estimation[END_REF], or multiple echo-time [START_REF] Veraart | Te dependent diffusion imaging (teddi) distinguishes between compartmental t2 relaxation times[END_REF].

Figure 3 shows three examples of the quality of the autoencoder RIF reconstruction for the three shells data with 15 neurons in the latent space (HCP subject 160123). The three considered RIF are the signal mean invariant I 0 , the degree 2 power spectrum invariant I 22 , and one of the new invariant developed by using the original framework presented in [START_REF] Zucchelli | A computational framework for generating rotation invariant features and its application in diffusion mri[END_REF]: invariant I 224 . As it is possible to see these three RIF were recovered almost perfectly by the autoencoder. Although we only show three examples this level of quality is consistent for all the 36 RIF in the three shells.

Figure 4 shows the 15 microstructural features obtained for the HCP subject 160123 using the 3-shells RIF. As it is possible to see, the different features highlight different parts of the brain tissues. Most of these microstructural features present a high-contrast between white matter and gray matter. In particular, some single-fiber white-matter areas like the corticospinal tract and corpus-callosum present the higher contrast in most of the features, similarly with classical fractional anisotropy maps. We must remind the reader that the RIF are very sensitive to the underlying fiber geometry. Therefore, it is not surprising that most of the autoencoder-based microstructural features can provide high contrast in the different white-matter regions. The work of Vasilev and colleagues [START_REF] Vasilev | q-space novelty detection with variational autoencoders[END_REF]] shows that the autoencoder-based features of the dMRI signal can be used to discriminate the voxels affected by multiple sclerosis. Although we based this work only on healthy subjects, we are positive that by combining autoencoders and rotation invariant features the resulting microstructural features might even be more precise in detecting pathological conditions.

Conclusions

In this work, we explore the use of the combination of autoencoder neural networks and RIF for the estimation of tissue microstructural features in PGSE dMRI data. There are two main advantages of using this technique. First, the use of RIF permits to eliminate the influence of tissue orientation in the microstructural features derived from the autoencoder latent-space. Second, the autoencoder is a model-free technique. This is particularly important in clinical data where the choice of the wrong model could lead to an incorrect characterization of the microstructural changes induced by pathology. One of the drawbacks of the technique is the reduced interpretability of the microstructural features. While "axonal density" and "extra-axonal diffusivity" are directly linked to the physical properties of the brain tissue, the autoencoder-derived features are more difficult to label. However, the main objective of these features is to help to identify a pathological condition in the brain tissue and be effectively used as biomarkers. Therefore, our future work will be aimed to test our technique on clinical data like multiple sclerosis and Alzheimer's disease.

2 ,

 2 90 volumes at b-values 1000 s/mm 2 , 90 volumes at b-values 2000 s/mm 2 , and 90 volumes at b-values 3000 s/mm 2 . Each volume is composed of 145 × 174 × 145 voxels with a resolution of 1.25 × 1.25 × 1.25 mm 3 . Before calculating the SH expansion, we normalized each volume dividing the diffusion signal by the mean of the b-value 0 s/mm 2 volumes.

Fig. 1 .

 1 Fig. 1. Graphical representation of the autoencoder design used in this work. The input layer counts 12 × number of b-values channels, each hidden layer is composed of 50 neurons, and the output layer has the same size of the input. For each data we trained several autoencoder with the number of neurons in the latent space ranging from 3 to the size of the input channel.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Log MSE of the difference between the input RIF and the network RIF as a function of the size of the latent-space, for 1 shell (b-value = 1000 s/mm 2 , on the left), two shells (up to b-value = 2000 s/mm 2 , at the center), and three shells (up to b-value = 3000 s/mm 2 , on the right). Each transparent line corresponds to a different random network initialization and the solid line represents the average for each of the three HCP subjects (in red, blue, and green). The black crosses correspond to the elbow point of the log curves.

Fig. 4 .

 4 Fig. 4. The 15 autoencoder microstructural features obtained from the 3-shells data for HCP subject 160123.

see the recent MEMENTO challenge: https://my.vanderbilt.edu/memento/

https://pytorch.org/

5 Acknowledgments This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002. This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 694665 : CoBCoM -Computational Brain Connectivity Mapping). Data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657)