N Gravielle

H Bundaleski

P Khemliche

Roncin

C A Ríos Rubiano

G A Bocan

M S Gravielle

N Bundaleski

H Khemliche

P Roncin

Ab initio potential for the He-Ag(110) interaction investigated using grazing-incidence fast-atom diffraction

Keywords: number(s): 79.20.Rf, 34.20.-b, 34.35.+a

Ab initio potential for the He-Ag(110) interaction investigated using grazing-incidence fast-atom diffraction

I. INTRODUCTION

An accurate interaction potential for atom-surface systems is central for achieving a correct description of several processes and phenomena in surface science. This better understanding is in turn most relevant for technological applications, for example in the chemical industry. For the particular case of He-Ag systems, early semi-empirical potentials, as well as more recent ones based on ab initio calculations, have been proposed and probed by means of diverse physical properties, such as adsorption coefficients or scattering cross sections [1][2][3][4]. Regarding the phenomenon chosen to test the potential, two important aspects need to be considered: the energy range it can probe and its sensitivity to the characteristics of the model potential. In this respect, grazing-incidence fast-atom diffraction (GIFAD) on crystal surfaces has recently emerged as a powerful surface analysis method that allows one to probe surface interactions with an exceptional sensitivity in the 10 meV to few eV energy range [5][6][7][8][9][10][11].

A simple way of understanding the unexpected presence of interference signatures in the distributions of scattered atoms, after swift and grazingly impinging on a surface along low-indexed crystallographic directions [6], is to separately consider the fast movement of the projectile parallel to the incidence channel and its slow motion in the corresponding perpendicular plane. It was found that this particular geometrical condition helps to prevent the quantum decoherence originated by surface thermal vibrations [12,13], making it possible to observe GIFAD effects for a wide variety of materials. In fact, GIFAD patterns have been observed for insulator surfaces [5,6], where electronic excitations are strongly suppressed due to the presence of a wide band gap, but also for semiconductors [7], metals [14,15], as well as adsorbate-covered metal surfaces [16,17].

In this work we present both experimental data and theoretical results for projectile momentum distributions associated with grazing scattering of fast He atoms on a Ag(110) surface. This system corresponds to the first and simplest metallic case for which GIFAD effects were experimentally observed [14]. Here, we make use of the exceptional sensitivity of GIFAD patterns to the shape of the interaction potential to thoroughly check the quality of an ab initio He-Ag(110) potential.

We describe the He-Ag(110) interaction potential with a full three-dimensional (3D) potential energy surface (PES) that takes into account the projectile's three degrees of freedom. Such a 3D PES is obtained from an accurate density-functional theory (DFT) calculation derived with the QUANTUM ESPRESSO code [18], which is specially adequate to perform selfconsistent calculations on metallic systems. We calculate a dense 3D grid of ab initio energies and obtain the complete 3D PES by means of a sophisticated interpolation technique [19].

The scattering process is evaluated within the surface eikonal (SE) approximation, which is a distorted wave method that makes use of the eikonal wave function to represent the elastic collision with the surface [20]. On the other hand, the motion of the fast projectile is classically described by considering axially channeled trajectories for different initial conditions. The SE approach is a semiclassical approximation that includes a clear description of the main mechanisms of the process. It has been successfully applied to investigate fast-atom diffraction from insulator surfaces [21][22][23]. In contrast to other approaches [24,25], the corrugation of the complete 3D PES is taken into account without averaging the projectile-surface potential along the incidence direction. SE projectile distributions using the DFT potential are compared with the experimental spectra for helium incidence along three different directions of the Ag(110) surface: [110],

[001], and [112]. The paper is organized as follows. The experimental method and the theoretical model, including details of the potential evaluation, are summarized in Secs. II and III, respectively. Results are presented and discussed in Sec. IV, and in Sec. V we outline our conclusions. Atomic units (a.u.) are used unless otherwise stated.

II. EXPERIMENTAL METHOD

The experimental setup has been described in detail elsewhere [14]. A primary 3 He + ion beam is extracted from an ECR ion source biased at the desired voltage, it is mass selected by a 180 • magnet before being neutralized in a charge exchange chamber. The beam enters the UHV chamber, crosses the interaction region, and is imaged by a positionsensitive detector (PSD) located half a meter downstream on the opposite side of the chamber. A set of movable apertures is used to limit the atomic beam divergence down to the desired angular resolution as monitored on the PSD. Then, the Ag(110) surface is gradually inserted into the beam. The direct beam progressively disappears while the scattered particles appear around the direction corresponding to specular reflection on the surface. The Ag crystal is then rotated around the surface axis until the beam is aligned along a low-index crystallographic direction. The Ag(110) surface is prepared by cycles of large angle 500 eV Ar + sputtering, annealing at 600 • K and grazing angle sputtering by 5 keV Ar + ions. The Ag(110) crystal is then inserted almost parallel to the 3 He beam. The projectile total energy is tuned between 300 eV and 1 keV whereas the incidence angle is adjusted between 0.5 and 2 degrees, producing energies normal to the surface in the sub-eV range.

III. THEORETICAL MODEL

A. Surface eikonal approximation

Within the surface eikonal model the transition matrix per unit area A reads [20]

T (SE) if = 1 A A d R os a if ( R os ), (1)
where R os determines the initial position of the projectile on the surface plane and

a if ( R os ) = 1 (2π ) 3 +∞ -∞ dt |v z ( R P )| × exp[-i Q. R P -iη( R P )] V SP ( R P ) (2)
is the transition amplitude associated with the classical path R P ≡ R P ( R os ,t), with v z ( R P ) the component of the projectile velocity perpendicular to the surface plane. In Eq. (2) Q = K f -K i denotes the projectile momentum transfer, with K i(f ) the initial (final) projectile momentum satisfying energy conservation (i.e., K f = K i ). The phase η is the eikonal-Maslov phase, which is defined along the projectile path as [22]

η( R P ) = t -∞ dt V SP ( R P (t )) + φ M , (3)
where V SP is the projectile-surface interaction and φ M = νπ/2 is the Maslov correction term, with ν the Maslov index as defined in Ref. [26]. Note that the surface eikonal approach takes into account the quantum interference between transition amplitudes corresponding to different axially channeled projectile trajectories. Details on calculations are given in Refs. [20,22]. From Eq. ( 1), the transferred momentum distribution is obtained as

dP (SE) dQ tr = (2π ) 4 m 2 P dθ f cos θ f K 2 f s -Q 2 tr T (SE) if 2 , ( 4
)
where m P is the projectile mass, K f s = K f cos θ f is the final momentum parallel to the surface, and Q tr = K f cos θ f sin ϕ f is the transverse momentum transfer, that is, the component of Q perpendicular to the incidence channel on the surface plane, with θ f and ϕ f the exit polar and azimuthal angles, respectively.

B. Projectile-surface interaction

The interaction energy of the He atom with the Ag(110) surface is described with a full adiabatic 3D PES that depends on the atomic position R P = (X,Y,Z). The PES is constructed from a grid of 252 ab initio energies over which an interpolation is performed.

All ab initio data are obtained from the DFT-based QUANTUM ESPRESSO code [18], which uses a plane-wave basis set to expand the system wave functions and is particularly efficient to model metallic surfaces. The values of relevant QUANTUM ESPRESSO parameters are chosen so that ab initio energies are calculated to a prescribed accuracy (differences <5 meV with respect to the converged result). The exchange-correlation energy is calculated within the generalized gradient approximation (GGA), using the Perdew-Burke-Ernzerhof energy functional (PBE) [27]. The electron-core interaction is described with ultrasoft pseudopotentials [28] (specifically, He.pbe-mt_fhi.UPF and Ag.pbe-d-rrkjus.UPF from www.quantum-espresso.org). The energy cutoff in the plane-wave expansion is 35.0 Ryd for the wave functions and 245.0 Ryd for the charge density and potential; the fractional occupancies are determined through the broadening approach of Marzari-Vanderbilt [29] with σ = 0.01 Ryd.; and the Brillouin-zone integration is performed with a 10 × 7 × 1 Monkhorst-Pack grid of special k points, ensuring a k grid equally dense in all directions. The Ag lattice constant, obtained from a bulk calculation, is a = 7.865 a.u., this value being slightly higher than the experimental one (7.709 a.u.).

The Ag(110) surface is modeled by means of the supercellslab scheme. A four-layer slab is used with a 2 × 2 cell in the plane parallel to the surface (atomic coverage of 0.25) and a supercell of length 25.03 a.u. along the normal to the surface (ẑ axis). The latter length guarantees that the energy for He-Ag(110), with the He atom midway between slabs, is XY independent and thus provides a reasonable representation of the asymptotic region. This configuration, for which there is hardly any interaction with the surface, is chosen as the energy reference. In order to get the surface equilibrium geometry, the interlayer distance is relaxed from its bulk value d 0 = 2.781 a.u. The two bottom layers are kept fixed during the relaxation process, which continues until the energy difference between consecutive iterations is less than 10 -4 a.u. and all components of all forces are less than 10 -3 a.u. Geometry corrections due to relaxation amount to -9.14% and +4.11% for the first and second interlayer distances respectively, in accord with experimental results [30]. Once relaxed, the slab is kept frozen for the calculations that follow.

Given the closed-shell electronic structure of the He atom (1s 2 ), a non-spin-polarized calculation is enough for an adequate description of the ground state. A 3D energy grid is obtained. It consists of 42 equidistant points Z in the range -2.0 a.u. Z 8.25 a.u., with Z = 0 corresponding to the topmost surface layer, and 6 sites [(X,Y ) values] uniformly spread within the unit cell, as indicated in Fig. 1. Once the 252-point energy grid is derived, the corrugation reducing procedure [19] together with spline interpolation over the ab initio data are used to build the PES. Further details on the calculations are given in the Appendix. In Figs. 2(a) and 2(b) we show equipotential curves for two different distances to the surface, Z = 3.0 a.u. and Z = 4.0 a.u., respectively, corresponding to the characteristic region of He closest approach to the surface for the incidence conditions considered in this work.

IV. RESULTS

In this article, experimental momentum distributions of 3 He atoms scattered from a Ag(110) surface under grazingincidence conditions are employed in order to probe the DFT potential energy surface. The evaluation of the SE transition matrix from Eq. ( 1) involves the integration on the initial position R os of the classical trajectory on a given area A. This integral was calculated with the Monte Carlo technique using more than 2 × 10 5 classical trajectories, with random R os values given by a Gaussian distribution, which avoids the unphysical interferences that arise in the momentum spectra when uniform random positions are used.

We start by considering incidence along the [1 10] channel, which presents the strongest corrugation and, consequently, the richest diffraction patterns. In order to provide an overall understanding of the GIFAD scenario, in Fig. 3 we compare experimental and theoretical diffraction charts for 3 He atoms impinging on the silver surface with an energy E = K 2 i /(2m P ) = 0.5 keV. Such diffraction charts display the intensity of the projectile distribution as a function of the transverse transferred momentum Q tr and the normal impact energy E ⊥ = E sin 2 θ i , where θ i denotes the incidence angle measured with respect to the surface plane. Note that different values of E ⊥ allow one to probe potential contours across the incidence direction for different distances to the surface. Therefore, the general good agreement between the experimental and simulated two-dimensional diffraction charts of Fig. 3 is a signature of the quality of the present DFT potential in the 0.05-0.3 eV range of perpendicular energies.

A closer examination of Fig. 3 reveals a series of maxima and minima structures in the SE diffraction chart, which are symmetrically placed with respect to the incidence direction, for which Q tr = 0. Similar patterns are also present in the experimental distribution, which displays a slight momentum asymmetry attributed to a small misalignment of the incident beam with respect to the low-indexed crystallographic direction. As it is the case with all diffraction techniques, these interference structures stem from two different contributions, the structure factor, responsible for the presence of the Bragg peaks, and the form factor, responsible for the intensity variations. Both mechanisms can be separately analyzed within the SE model by factorizing the transition matrix as [23] T . The T (SE) 1 (Q tr ) factor is the result of the interference between projectile trajectories whose initial positions are separated by a distance smaller than d and presents an oscillatory structure as a function of Q tr , which is associated with the projectile deflection by the surface potential. When a large number of diffraction orders is observed, these oscillatory structures are sometimes called supernumerary rainbows [25]. Thus, Bragg peak positions provide information about the crystallographic structure only, but their intensities are extremely sensitive to the shape of the PES across the incidence channel. What is more, it was found that small changes in the surface potential can modify strongly the supernumerary rainbow factor, even transforming a GIFAD maximum into a minimum [21,24].

(SE) if = T (SE) 1 (Q tr ) S n tr (Q tr ), (5)
In Fig. 4 we plot experimental and theoretical GIFAD patterns for 0.5 keV He atoms impinging along the [001] direction. Perpendicular energies in the range from 0.20 to 0.32 eV are considered. In the figure, SE differential probabilities obtained by considering an integration region A equal to 8 × 8 reduced unit cells are displayed with a blue dashed line. These results present narrow Bragg maxima, whose positions are indicated with vertical dashed lines in the figure. Note that the width of these Bragg peaks is governed by n tr (in our case, n tr = 8), and it gets smaller as n tr increases. On the other hand, since the factor T (SE) 1 determines the number of observed Bragg maxima, the smaller corrugation of this crystallographic channel results in narrower spectra than the ones of Fig. 3 for the same normal energies.

In order to simulate the experimental conditions and establish a comparison with experimental data, in Fig. 4 we also plot SE differential probabilities convoluted with a Lorentz function (blue solid line), where the parameters of the line broadening were obtained from the observed linewidths [14]. This convolution takes into account not only the experimental divergence of the incident beam but also the broadening introduced by both thermal vibrations of lattice atoms and inelastic processes, which contribute to deteriorate the coherence [31]. We find a fairly good agreement between convoluted SE results and experimental intensities, with the exception of the Bragg maxima of order m = ±2 for E ⊥ = 0.32 eV, whose intensities are overestimated by the SE curve. This fact is related to the presence of sharp rainbow peaks in the T (SE) 1 function, which are originated by the classical description of the projectile motion [32]. Such rainbow peaks affect the outermost Bragg maxima when they are close to each other, as observed in Fig. 4(d), but this deficiency does not have any influence on the SE distribution at smaller deflection angles.

In Fig. 5, experimental momentum spectra for 1 keV He atoms impinging on a silver surface along the [001] direction are compared with SE distributions obtained both including and not including the convolution. As reported in previous works [21,24], the intensity of the central peak, corresponding to the Bragg maximum of order m = 0, is extremely sensitive to the corrugation amplitude of the surface potential. In a simple Hard-Wall model this trough to peak amplitude is called h and the intensity of the Bragg maximum of order m = 0 oscillates with the geometrical phase = K i⊥ h, where K i⊥ = K i sin θ i is the initial projectile momentum perpendicular to the surface [33]. This means that the corrugation amplitude h is probed with an accuracy, which is a fraction of the perpendicular wavelength. In the figure we observe how a small variation in the E ⊥ value strongly modifies the relative intensity of the zeroth-order Bragg peak, which switches from being almost an absolute maximum for E ⊥ = 0.21 eV to a minimum for E ⊥ = 0.29 eV. This behavior is well reproduced by the simulations based on the DFT potential, indicating that the ab initio model gives an appropriate description of the He-Ag(110) interaction in this normal energy range.

To complete the analysis of the surface potential we consider an He beam impinging along the [112] channel with E = 1 keV. In Fig. 6 experimental momentum distributions for normal energies ranging from 0.30 to 0.44 eV are compared with convoluted and unconvoluted SE spectra. For E ⊥ = 0.44 eV [Fig. 6(c)] we find that the SE momentum spectrum, including experimental and inherent uncertainties through convolution, is in very good agreement with the experimental data in the whole Q tr range. But, since the Q tr position of the classical rainbow peak depends on E ⊥ , when the normal energy decreases and the classical rainbow maximum becomes close to the outermost Bragg position, as observed from the unconvoluted values of Fig. 6(a), the corresponding Bragg intensity starts to be overestimated. However, as explained above, this failure is a characteristic of the semiclassical methods and it is not related to the potential model. Furthermore, taking into account that the number of observed Bragg maxima depends strongly on the shape of the equipotential curves across the incidence channel [24], the general agreement between the experimental and simulated diffraction spectra of Fig. 6 can be considered as an indication of the performance of the present DFT potential in such a perpendicular energy range.

Finally, in Fig. 7 we plot the effective corrugation z of the potential across the different channels, [110], [001], and [112], as a function of the perpendicular energy E ⊥ . The corrugation z is defined as the normal distance between the maximum and the minimum of an equipotential surface that is obtained by averaging the 3D PES along the incidence direction. Experimentally derived data, displayed with full symbols in the figure, were deduced from Bragg peak intensities under different angles of incidence by using the Hard-Wall model [33]. In all the cases considered the theoretical curve is close to the experimentally derived data, showing an increase of the corrugation with the normal energy, in opposition to what happens for the LiF(001) surface along the [110] direction [22].

V. CONCLUSION

We have addressed, both experimentally and theoretically, GIFAD patterns for He atoms scattered from a Ag(110) surface in order to probe an ab initio surface potential that was obtained from DFT by making use of the QUANTUM ESPRESSO code. Momentum distributions of scattered projectiles were derived from the SE approximation, which takes into account the quantum interference in terms of the coherent superposition of transition amplitudes for different projectile paths with the same final momentum.

SE results for incidence along three different crystallographic directions of the Ag(110) surface are compared with experimental data considering a wide range of energies associated with the motion normal to the surface plane. Several aspects of the experimental momentum distributions are well reproduced by the theoretical simulations. Taking into account the extreme sensitivity of fast atom diffraction patterns, we thus conclude that the present ab initio DFT surface potential provides an appropriate description of the He-Ag(110) interaction for perpendicular energies in the E ⊥ = 0.1 eV-0.5 eV range. In addition, we have investigated the effective corrugation across the different incidence channels, as a function of the perpendicular energy, finding good agreement with the experimentally derived corrugation.

APPENDIX: ON THE HE-AG(110) INTERACTION POTENTIAL PRECISION

In this Appendix we present supplementary information on the interaction potential used in this work. Figure 8 shows the ab initio energies as a function of Z, for the six XY grid sites (Fig. 1). As it was pointed out in the main text, these calculations were performed for a set of parameters carefully chosen so that the energy differences were well converged within a 5 meV margin. Mark the energy's XY asymptotic independence (see Fig. 8 inset), which is indicative of an enough amount of vacuum between slabs.

Once the XY Z grid is obtained, the corrugation reducing procedure together with spline interpolation are used to build up the PES as a function of the He atom XY Z coordinates. Within this scheme, the quality of the interpolation will be The quality of the interpolation is checked by contrasting interpolated to test ab initio calculations for off-grid trajectories along X, Y , and Z. The results are shown in Fig. 9. On analyzing this figure we must keep in mind that the region closest to the surface, being the one with the highest corrugation, is also the most challenging region for the interpolation. Therefore, in the middle and bottom panels we have considered trajectories along the [1 10] and [001] directions for a fixed in-grid Z = 2.5 a.u. ∼ 1.3 Å. The differences observed (right column) are all below 10 meV. Regarding this result, note that for every dynamic simulation performed in this work the closest distance to the surface reached by the He atoms was Z min 2.8 a.u., where the interpolation will perform even better than at the tested distance. In the top panel of Fig. 9 we have considered a trajectory along Z for an off-grid site, corresponding to the midpoint between a TOP and HOLLOW site. We see here that in the Z region of interest for the dynamics, the differences between interpolated and ab initio results are around a few meVs.

Finally, we reckon it worth mentioning that the DFT calculation 5 meV convergence margin or the interpolation error do not by themselves represent how well the PES will compare to experimental data. There are several other factors that might affect the PES precision. The forced periodicity of the supercell-slab scheme, the approximation used to include exchange correlation energy in the DFT calculation or the neglect of Van der Waals interactions, although not expected to play an important role in the process under study, are other sources of error to keep in mind.

FIG. 1 .

 1 FIG. 1. (Color online) Depiction of the atomic positions of the Ag(110) surface, with a the lattice constant. The six (XY ) sites shown in the figure correspond to the ones used for the PES calculation.Circles: TOP and HOLLOW sites standing for atoms of the first and second layers, respectively. Squares: BRIDGE_1 and BRIDGE_2 sites corresponding to the middle points between first and second TOP atomic neighbors, respectively. Triangles: middle points between the mentioned sites. The incidence directions of He atoms are also indicated.

FIG. 2 .

 2 FIG. 2. (Color online) Equipotential curves for the interaction between the He atom and the Ag(110) surface, as a function of the coordinates X and Y , along the [1 10] and [001] directions, respectively. Two different distances to the surface are considered: (a) Z = 3.0 a.u., and (b) Z = 4.0 a.u.

FIG. 3 . 1 (

 31 FIG. 3. (Color online) Two-dimensional diffraction chart displaying (a) experimental intensities and (b) SE differential momentum probabilities, as a function of the normal energy E ⊥ and the transverse momentum transfer Q tr , for 0.5 keV 3 He atoms impinging along the [110] direction. where T (SE) 1 (Q tr ) is derived from Eq. (1) by evaluating the R os integral over one reduced unit cell. For incidence along the [1 10] and [001] channels, the second factor reads S n tr (Q tr ) = sin(n tr β)/ (n tr sin β), where n tr denotes the number of reduced unit cells along the transverse direction that are included in the integration region A, and β = Q tr d/2, with d the spacial periodicity of the channel. The S n tr (Q tr ) factor is due to interference between projectile trajectories whose initial positions are separated by a distance just equal to d and gives rise to Bragg peaks placed at Q tr = m2π/d, with m an integer number. From Eq. (5), the intensities of such Bragg maxima are completely determined by the envelope factor, T (SE) 1

FIG. 4 .FIG. 5 .

 45 FIG. 4. (Color online) Momentum distributions, as a function of the transverse momentum transfer Q tr , for 0.5 keV 3 He atoms impinging on Ag(110) along the [001] direction. The normal incidence energy is: (a) 0.20 eV, (b) 0.24 eV, (c) 0.28 eV, and (d) 0.32 eV. Empty circles: experimental data. Solid blue line: SE differential probability convoluted to include inherent uncertainties. Dashed blue line: SE differential probabilities without convolution, as explained in the text. The vertical dashed lines show Bragg peak positions and the arrow indicates the position of the classical rainbow maximum.

FIG. 6 .

 6 FIG. 6. (Color online) Similar to Fig. 4 for 1 keV 3 He atoms impinging along the [1 12] direction with normal incidence energy (a) 0.30 eV, (b) 0.37 eV, and (c) 0.44 eV.

FIG. 7 .

 7 FIG. 7. (Color online) Effective corrugation z of the potential across the low-indexed crystallographic channels [1 10], [001], and [1 12], as a function of the perpendicular energy. Full symbols, experimentally derived results in Hard-Wall approximation; solid curves, values derived from the present DFT potential.

FIG. 8 .

 8 FIG. 8. Energy vs Z for the six XY sites used to make the interpolation grid. Z = 0 corresponds to the position of the topmost atomic layer of the Ag(110) surface. A detail of the asymptotic region is shown in the inset.

FIG. 9 .

 9 FIG. 9. (Color online) Interpolation checks. Full circles, ab initio data used in the interpolation; red crosses, ab initio data used to test the interpolation procedure. Top: Energy (left) and difference between interpolated and ab initio energy (right), as a function of Z, for the off-grid top-hollow site, located at mid-distance between a TOP and a HOLLOW site. Middle: Energy (left) and difference between interpolated and ab initio energy (right), as a function of X, for Z = 2.5 a.u. ∼ 1.3 Å and Y = 0. Bottom: Energy (left) and difference between interpolated and ab initio energy (right), as a function of Y , for Z = 2.5 a.u. ∼ 1.3 Å and X = 0.

ACKNOWLEDGMENTS C.R.R. and M.S.G. acknowledge financial support by CONICET, UBA, and ANPCyT of Argentina. M.S.G. would like to express gratitude to the Université Paris-Sud for support during a portion of this work. G.A.B. acknowledges financial support by ANPCyT and is also thankful to H. F. Busnengo, J. D. Führ, and M. L. Martiarena regarding the PES calculation. The experimental work was supported by the ANR under Contract No. ANR-07-BLAN-0160-01.


B Lalmi

H Khemliche

A Momeni

P Soulisse

P Roncin

High resolution imaging of superficial mosaicity in single crystals using Grazing Incidence Fast Atom Diffraction

Keywords: Diffraction, Fast Atom Diffraction, mosaic structure, atom scattering, particle-surface interactions, ionic insulators

A new table top technique is used to simultaneously analyse the local morphology of crystalline surfaces as well as the misalignment of large scale domains at the topmost surface layer. The approach is based on fast atom diffraction at grazing incidence (GIFAD); the diffraction pattern yields the topology of the surface electronic density, similar to that imaged by Atomic Force Microscopy. If superficial mosaicity is present, diffraction patterns arising from each mosaic domain can be distinguished, providing high sensitivity on the property of each domain. Taking NaCl(001) as an example, we observe a discrete tilt angle distribution of mosaic domains following an arithmetic progression with a 0.025°±0.005 difference; a twist mosaic angle of 0.09°±0.01 is also observed.

Introduction

Actual single crystals always contain defects that induce local distortions to the periodic arrangement.

These defects can take several forms: point defects (interstitial atoms, vacancies), linear (dislocations), two dimensional (stacking faults, antiphase boundaries), and three dimensional (amorphous cluster, precipitates) [1][2][3][START_REF] West | Solid State Chemistry and its Applications[END_REF]. These defects also induce local modifications in the electronic structure, which can affect considerably the physico-chemical properties [START_REF] Schulze | [END_REF]. Mosaicity is one of the fundamental defects most often present in single crystals. The concept was first introduced to explain the decrease of intensity in the X-ray reflectivity studies [6][7][8]. Mosaic crystals are considered to be formed by a large number of small perfect crystallites slightly disoriented with respect to each other. This is due to an arrangement of dislocations in a pattern separating the crystal into blocks [9,[START_REF] Kosevich | Dislocations in Theory of Elasticity[END_REF]. The spread of misalignment angles of the crystal blocks gives the degree of mosaicity of the material. Generally, mosaicity can be described in terms of two crystallographic parameters: the out-of plane misorientation called tilt mosaic and the in-plan misorientation called twist mosaic [START_REF] Chierchia | [END_REF]12].

In X-ray diffraction, tilt and twist mosaic are given by the FWHM of symmetric and asymmetric -scan rocking curve respectively [13][14][15][16]. However, additional structural imperfections and instrumental broadening also affect the linewidth [START_REF] Guinier | Theorie et Technique de la Radiocristallographie[END_REF][START_REF] Fatemi | [END_REF][START_REF] Fewster | X-ray and Neutron Diffraction: Theory and Applications[END_REF][START_REF] Hordon | [END_REF], imposing measurement for several diffraction orders. Such treatment requires a comparatively high precision and the determined values are often marred by considerable uncertainty. Furthermore, because of the penetration of X-rays in solids, the results are insensitive to the superficial layer, even in grazing geometry which probes a typical depth around 100 Å. Local microscopy such as STM or AFM are excellent tools to characterize such crystalline distortion but only in a comparatively reduced scale because atomic resolution cannot be reached on scales larger than few 100 nm so that superficial mosaicity, is still largely unknown.

In this paper, we present a quantitative characterization of superficial mosaicity, using an exclusively surface sensitive technique, Grazing Incidence Fast Atom Diffraction (GIFAD). This technique is based on the coherent scattering of helium atoms at energies in the keV range and at grazing incidence angle around 1°. In this geometry, the incident atom energy normal to the surface is below 1 eV, thus precluding penetration below the surface.

Experimental details

All experiments were performed in an ultrahigh vacuum chamber with a base pressure in the 10 -10 Torr range.

He + ions are extracted from an ion source at 450 eV energy; they are subsequently neutralized by charge exchange in a He cell. Then two variable apertures are used to reduce the beam divergence and the beam size to less than 1 mrad and 100µm respectively. In the specific experimental conditions used here, the divergence of the direct beam is 0.3 mrad (0.017°). The beam scattered from the sample, and a small fraction of the direct beam, reach a detector located at L=865 mm. The latter is of made of microchannel plates and a phosphor screen; images are recorded by a high resolution CCD camera.

In addition to providing the lattice parameter, the intensity distribution within the diffraction pattern reveals the profile of the valence electronic density, thus providing a map of the surface topography comparable to that obtained in AFM but by via reciprocal space technique i.e. with interferometric accuracy, see e.g. [21][22][23]. The exact beam position, required for an accurate determination of the incident angle, is recorded with the sample surface in a fully retracted position. The commercial NaCl(001) sample, 10 x 10 mm in size, was air-cleaved and then rapidly introduced in the vacuum chamber; it was further annealed in UHV at 600°C.

Results and discussions

Figure 1 displays diffraction images obtained under an incidence angle close to 0.9 degree and correspond to a primary beam aligned along [-110] and [100], the two main surface channeling directions. It is worth recalling that GIFAD does not allow for exchange of reciprocal lattice vectors along the beam direction. This is due to the fact that the coupling strength is too weak and the associated energy too high [22,24]. As a result, only one Laue circle intercepting the diffraction spots and the direct beam can be observed. Each image contains two types of information, one related to the macroscopic orientation of the surface and the other connected to the microscopic atomic structure.

The peak spacing yields reciprocal vectors of 1.57 ± 0.02 Å -1 and 2.22 ± 0.02 Å -1 which correspond well to the projected lattice parameter along the [-110] and [100] directions of the NaCl (001) surface respectively. The geometric description of each image can be summarised with three points, the beam position B, the specular spot position S and the center C of the Laue circle. The line BS gives the direction of the incidence plane which contains the surface normal while its perpendicular bisector is the projection of the surface plane and contains C.

When the beam is perfectly aligned along one of the sample crystallographic axes, C is located on BS.

The intensity of the spots is given by the form factor that, within the Hard-Wall description [25], is simply the modulus of a Fourier-like transform of the corrugation function. Here, the observed intensity is well fitted by a simple cosine corrugation with a half amplitude of 0.29 Å, consistent with the value measured with thermal helium scattering [26]. More important, the measured intensity distribution is also extremely sensitive to the angle of incidence (the radius of the Laue circle) and to the misalignment angle between the beam and the crystallographic axis, so that the analysis provides additional evidence of any in plane misalignment.

Despite the small beam size, the grazing geometry results in a beam length on the sample of several millimetres. At 0.9°incidence, this full length is typically 6 mm, though the actual length is smaller because the sample surface is not fully inserted into the beam, thus a small fraction of the incident beam reaches the detector and provides an angular reference. GIFAD patterns of the same sample but at a different location on the sample surface, recorded under an incidence angle close to 0.5°, are shown in fig. 2. At this angle, the projected length of the incident He beam onto the surface is nearly twice larger than it is at 0.9°.

The image in fig. 2 The image in fig. 2.b, recorded with the incident beam along the [100] direction, does not contain any sign of tilted domains with rotation vector parallel to the <110> direction.

The sensitivity of the diffraction pattern to mosaic structures of the surface is illustrated in fig. 3. The presence of tilted domains, whether the rotation axis is normal or parallel to the beam, is depicted in fig. 3a and3b respectively. These schemes describe well the diffraction images of fig. 2a and2b; along the [-110] direction, the fit to the data in fig. 2a where D is the diameter of the Laue circle. Thus rotation angles below one degree would be difficult to resolve.

In fact the quasi-gaussian broadening, within the Laue circles, of the diffraction spots observed in fig. 2 could well originate from such a mosaicity. Assuming that images in fig. 1 originate from a single domain and that the extra broadening in fig. 2 is due to rotation of the domains around the incident beam direction, an angular distribution with 1.2°FWHM is deduced.

A third situation can be encountered when the atomic beam interacts with twisted mosaic domains; this is shown in fig. 3c. Here the Laue circles arising from the twisted domains are displaced by a distance ξ d proportional to the twist angle. An example of this type of mosaicity is given in fig. 3d, it was recorded on a cleaved KBr(001) crystal with the beam along the [110] direction.

In the framework of quasi-elastic scattering, which is the case in GIFAD [25], and if we assume that the surface is perfectly flat without defects, diffracted intensity distributions presented in fig. 2 to 4 would not be authorized according to the energy conservation law. Thus such structured intensity distributions could only be related to topographic features of the surface.

Features similar to those shown in fig. 2a, with comparable tilt angle values, have been observed by high resolution X-Ray diffraction on cleaved NaCl(001) [27] and RbI single crystals [28]. Discrete mosaic distributions have also been reported on thin pseudomorphic ZnSe films on GaAs(001), above a critical thickness of 100 nm [29,30]. In the later studies, the tilt angle of the mosaic domains, derived from the positions of satellites spots in SPA-LEED measurements, is found to be 0.2°with respect to the macroscopic plan, which is only slightly larger than the resolution limit in SPA-LEED. With GIFAD, tilt angles as small as 0.01°can be measured, limitation comes only from the divergence of the primary beam. The latter can always be improved, but at the expense of a reduced atom flux that would increase the acquisition; the data shown here were typically acquired in 10-100 seconds.

The highly regular distribution of tilted mosaic domains described in fig. 2a is present on a large fraction of the cleaved NaCl(001) sample surface in the form of patches of mm size, meaning that some macroscopic regions of the surface are nearly free from such defects. In addition to this non-uniform patch distribution, we also observe a strong anisotropy in tilt mosaicity. Along the [-110] and [110] directions, the diffraction images bear the signature of tilted domains that have a finite component of their in-plane rotation vector normal to theses directions. Diffraction images measured along the [100] direction instead show that the domains have no rotation component in this direction. Consequently, we derive that the tilt vector should be parallel to the <010> direction.

On the origin of the observed mosaicity, some information may be gained by comparing results form samples having different shapes. Using the same sample preparation procedure, we clearly observe that round samples yield more mosaicity than square ones. Finally polished KBr(001) samples, used as a X-ray monochromator window, showed the lowest mosaicity. It is worth mentioning that all these crystals were purchased from the same supplier. It is interesting to note that CaF 2 (111) samples, cleaved with similar procedures, do not exhibit any measurable mosaicity.

Finally, this uniform distribution of the tilt mosaic domains within patches appears to be a robust characteristic of some ionic crystals and points to some kind of self-organisation of these domains.

An example of a peculiar organization mode is shown in fig. 4. Here, the probe beam is off axis with respect to the main crystallographic directions, exactly 7°away from the [0-10] direction. The surface corrugation (in terms of valence electron density) as averaged by the beam, is very flat so that only the zeroth-order is observed. As in fig. 2a, a series of equally spaced diffraction spots can be observed, each spot corresponding to specular reflection from a single mosaic domain. A projection along the vertical axis shows two regions, each described by series of peaks equally spaced by 0. 027°; however these two regions (marked by I and II in fig. 4) are phase shifted by a half period. Variations induced by moving the beam across the sample surface indicates a main domain size of a few 100 microns, this is consistent with the number of illuminated domains at a given incidence angle.

Conclusions

In conclusion, we have demonstrated that Grazing Incidence Fast Atom Diffraction (GIFAD) can provide a microscopic as well as long range structural analysis of the top most surface layer. In addition to quantitative information on surface corrugation, it can provide clear signature of surface topographic imperfections.

Qualitatively, and quantitatively, we have reported a first experimental characterization of the superficial mosaicity on a NaCl(001) single crystal. Two types of superficial mosaic structures have been observed; (i) twist mosaic defined by macroscopic blocks lying parallel to the (100) surface but slightly rotated one from another along the [100] direction, and (ii) tilt mosaic domains with a misalignment angle being a multiple of an elementary angle equal to 0.025°.

Because GIFAD is a particularly fast technique than can be operated in situ during thin film growth, it may become a valuable tool to follow in real time the formation dynamics of mosaic domains in constrained multilayered systems.

Figure captions

  .a was obtained along the [-110] axial surface channelling, it shows an unexpected superposition of diffraction patterns. One can clearly distinguish a series of diffraction spots, forming well defined Laue circles whose centers are aligned in the specular plane (see schematic representation fig.3.a). The different diameters are associated with different effective angles of incidence indicating the presence of different domains within the surface illuminated by the atomic beam. The fact that they appear at different incidence angles indicates that these domains are tilted with a component of the rotation vector located in the surface plane and perpendicular to the beam direction, i.e. parallel to the <100> direction. Similar patterns are observed when the beam is oriented along the [110] axis. The profile through the line in fig.2a is shown in the insert; the multipeak structure is very well described by a regular series of peaks spaced by 0.025°±0.005. We observe similar features on cleaved KBr(001) samples, with a characteristic tilt angle comparable to the value derived for NaCl(001).

  When the incident beam is parallel to the rotation angle of the tilted domains (see fig.3b.), Laue circles arising from each domain are simply rotated around the beam position by an amount γ equal to the angle between blocks. Then the separation between subsequent spots is given by D γ ×

Figure 1 :

 1 Figure 1:GIFAD patterns recorded with 450 eV4 He atoms on NaCl(001), under incidence angles of (a) 0.87°along [-110] and (b) 0.85°along[100]. The Laue circle passing though the spots is centered on C, intercepts the zeroth order on S and the direct beam on B at the lower part of the image. Note that a large fraction of the direct beam is cut to protect the detector.

Figure 2 :

 2 Figure 2: GIFAD patterns recorded with 450 eV 4 He atoms on NaCl(001), under incidence angles of 0.5°along the [-110] (a) and [100] (b) directions. The insert in (a) shows a line profile and a fit using equally spaced Gaussians.

Figure 3 :

 3 Figure 3:Sensitivity of the diffraction pattern to mosaic defects. Tilt mosaic with a rotation vector in the surface plan, perpendicular (a) and parallel (b) to the beam. Twist mosaic with a rotation vector normal to the surface (c), an experimental example is shown in (d).

Figure 4 :

 4 Figure 4: Diffraction profile measured off axis (7°way from [0-10]); only the zeroth order diffraction spot can be observed (see text). The intensity distribution along the vertical direction indicates two regions of equally spaced peaks, arising from individual tilt mosaic domains separated by 0.055°/2); these two regions are offset by half spacing.

Figure 1 Figure 4 (

 14 Figure 1 (figure-1.eps)

Acknowledgment

This work was supported by the ANR under contract number ANR-07-BLAN-0160-01.


Mohammad Veisi

Farhad Sobouti

Sébastien Chevrot

Madjid Abbasi

Esmaeil Shabanian

Upper mantle structure under the Zagros collision zone; insights from 3D teleseismic P-wave tomography

Keywords: Traveltime tomography, Zagros Mountains, Iranian Plateau, Continental collision, Slab detachment, Lithospheric delamination

We map the lateral variations of P-wave velocity beneath the Zagros collision zone in Iran down to 800 km depth by regional travel time tomography. We invert 32,293 relative travel time residuals from five temporary seismic networks across the collision zone in NW and central Zagros and northern Iran, and a large number of permanent stations distributed over the Iranian Plateau. A salient feature of our model is a ~225 km thick lithosphere beneath the Zagros, almost twice as thick as in the rest of Iran. Slab detachment at a depth range of 250-300 km from the base of the subducted continental Arabia is clearly distinguishable in central Zagros. In NW Zagros detachment appears to be at an earlier stage. Beneath the central Iran Plateau, we observe

J o u r n a l P r e -p r o o f

Journal Pre-proof shallow low-velocity regions down to 200 km depth, and smaller patches with significantly higher velocity in the 200-400 km depth range. This pattern could suggest a post-collisional foundering of the mantle lithosphere in the upper plate. The subducted slab seems to penetrate into the lower mantle in a segmented fashion, with no evidence for a stagnant slab above the 660 km interface.

Introduction

The Iranian Plateau has formed as a result of the closure of the Neo-Tethys Ocean and the Arabia-Eurasia continental collision. Some of the more recent estimates put the onset of collision in the Oligocene [START_REF] Egan | Subsidence and uplift mechanisms within the South Caspian Basin: insights from the onshore and offshore Azerbaijan region[END_REF][START_REF] Agard | Zagros orogeny: a subduction-dominated process[END_REF][START_REF] Mcquarrie | Retro-deforming the Arabia-Eurasia collision zone: Age of collision versus magnitude of continental subduction[END_REF][START_REF] Pirouz | Early Neogene foreland of the Zagros, implications for the initial closure of the Neo-Tethys and kinematics of crustal shortening[END_REF]. The Zagros Mountains on the southwestern margins of the plateau are the current locus of the collision in Iran (Fig. 1), accommodating a third of the total convergence of the Arabian platform under central Iran. The Main Zagros Thrust (MZT), the surface expression of the intracontinental suture (e.g., [START_REF] Stöcklin | Structural history and tectonics of Iran: a review[END_REF][START_REF] Paul | Seismological evidence for crustal-scale thrusting in the Zagros mountain belt (Iran)[END_REF]2010), defines the northeast boundary of the Zagros range. Coeval with the deformation of the Zagros, mountain building processes and volcanism have been active in the Alborz Mountains in the north (e.g., [START_REF] Guest | Late Cenozoic shortening in the westcentral Alborz Mountains, northern Iran, by combined conjugate strike-slip and thinskinned deformation[END_REF]2007;[START_REF] Ballato | Middle to late Miocene Middle Eastern climate from stable oxygen and carbon isotope data, southern Alborz Mountains, N Iran[END_REF], while central Iran has experienced relatively recent plateau formation (5-10 Ma) and post-collisional volcanism [START_REF] Omrani | Arcmagmatism and subduction history beneath the Zagros Mountains, Iran: a new report of adakites and geodynamic consequences[END_REF][START_REF] Hatzfeld | Comparisons of the kinematics and deep structures of the Zagros and Himalaya and of the Iranian and Tibetan plateaus and geodynamic implications[END_REF][START_REF] Chiu | Zircon U-Pb age constraints from Iran on the magmatic evolution related to Neo-Tethyan subduction and Zagros orogeny[END_REF].

A thorough understanding of the Zagros collision zone and its Late Cenozoic history requires detailed tomographic images of upper mantle structures. Several questions, in particular, are J o u r n a l P r e -p r o o f important in this regard; (1) What is the amount of lithospheric thickening since the onset of collision [START_REF] Paul | Seismological evidence for crustal-scale thrusting in the Zagros mountain belt (Iran)[END_REF] Has this lithospheric thickening been limited to the margins of the upper plate, or did it propagate far into the plate interior in central Iran? (3) Has the collision resulted in an effective detachment of the oceanic slab from the underthrusting Arabian continental margin? (4) If the Iranian Plateau is the result of a late-stage lithospheric foundering/delamination, then, is there any evidence of that process left anywhere in the upper mantle?

Many body and surface wave tomographic studies from the global scale (e.g., [START_REF] Van Der Meer | Atlas of the underworld: Slab remnants in the mantle, their sinking history, and a new outlook on lower mantle viscosity[END_REF]) to the regional scale (e.g., [START_REF] Chang | Joint inversion for threedimensional S velocity mantle structure along the Tethyan margin[END_REF]Koulakov, 2011, Maggi and[START_REF] Maggi | Surface waveform tomography of the Turkish-Iranian plateau[END_REF][START_REF] Al-Lazki | Pn tomographic imaging of mantle lid velocity and anisotropy at the junction of the Arabian, J o u r n a l P r e -p r o o f Journal Pre-proof Eurasian and African plates[END_REF] have mapped the upper mantle seismic velocities and revealed the large-scale structure of the Tethyan subduction zone beneath the Zagros and Iranian Plateau.

These studies have consistently found a high-velocity and thick lithosphere beneath the Zagros, whereas central Iran is characterized by a low-velocity and thin lithosphere.

Local 2-D tomographic imaging of Iran was first performed by [START_REF] Kaviani | A strong seismic velocity contrast in the shallow mantle across the Zagros collision zone (Iran)[END_REF] and many other studies have followed since (e.g., [START_REF] Paul | Seismic imaging of the lithospheric structure of the Zagros mountain belt (Iran)[END_REF][START_REF] Motaghi | High resolution image of uppermost mantle beneath NE Iran continental collision zone[END_REF][START_REF] Rahmani | Deep velocity image of the north Zagros collision zone (Iran) from regional and teleseismic tomography[END_REF][START_REF] Mahmoodabadi | Shear wave velocity structure of the upper-mantle beneath the northern Zagros collision zone revealed by nonlinear teleseismic tomography and Bayesian Monte-Carlo joint inversion of surface wave dispersion and teleseismic P-wave coda[END_REF]. Several 3-D studies also focused on different regions of Iran (e.g., [START_REF] Shomali | Lithospheric structure beneath the Zagros collision zone resolved by non-linear teleseismic tomography[END_REF][START_REF] Shad Manaman | New constraints on upper-mantle S-velocity structure and crustal thickness of the Iranian plateau using partitioned waveform inversion[END_REF]Mahmoodabdi et al., 2019). These previous studies, which suffered from a rather limited station coverage, have confirmed the geometry of the underthrusted Arabian margin beneath central Iran and of the subducted oceanic slab. However, the available global tomographic models have poor resolution in the shallow upper mantle beneath Iran, mainly because of their coarse spatial discretization. On the other hand, regional-scale 2-D models, owing to their limited aperture (<600 km), have a poor resolution at depth, in particular in the mantle transition zone. In this study, we combine several datasets from permanent and temporary seismic networks installed in Iran to improve both the J o u r n a l P r e -p r o o f aperture and ray coverage. The new extended dataset allows us to significantly improve the spatial resolution in both the shallow and deep parts of the mantle, providing key insights into the interaction of the two converging plates in the subduction zone and the fate of the subducted slab.

Data

Our dataset consists of teleseismic P-wave arrival times from 202 temporary, and 102 permanent stations in Iran (Fig. 1 and Table 1).

Data preparation and inversion method

We picked 32,293 first P-wave arrival times from 1588 teleseismic events with magnitudes greater than 5.5 and epicentral distances between 28 o and 95 o (Fig. 2a). For picking the relative P travel times we used the method developed by [START_REF] Chevrot | Optimal measurement of relative and absolute delay times by simulated annealing[END_REF]. The advantage of this method is that we can obtain robust and accurate travel time measurements even on very noisy data. The J o u r n a l P r e -p r o o f travel times need to be corrected for crustal thickness, station elevation, and ellipticity of the Earth. For crustal thickness correction, we used the Moho depth map from [START_REF] Kaviani | Crustal and uppermost mantle shear wave velocity structure beneath the Middle East from surface wave tomography[END_REF].

We improved that map by including new crustal thicknesses estimated from P-wave receiver functions obtained from the temporary stations. We gridded the surface of the model with 15×15 km blocks. Inside each block, we determined the crustal thickness by interpolating from the nearest points available from the receiver functions and the Moho map. The station elevation corrections were calculated by computing the travel times between the surface of the WGS84 ellipsoid and the station elevations. We also corrected the travel times for the ellipticity of the Earth. We then computed average event travel times that were removed from the travel time residuals. These relative P-wave travel times were then inverted following the method of [START_REF] Aki | Determination of the three-dimensional seismic structure of the lithosphere[END_REF], known as the ACH method.

Model parameterization in any inversion must follow some rules. The cell size in the uppermost layer of the model should not be smaller than the smallest distance between the stations, the depth to the bottom of the model should be no more than 2/3 of the largest station distances, and the cells should be large enough to contain sufficient ray crossings to prevent smearing in the output model [START_REF] Evans | Teleseismic velocity tomography using the ACH method: theory and application to continental scale[END_REF]. We found that the choice of 25-km cubic cells was a good compromise and our model has been parameterized with 33 layers from the surface down to 825 km depth.

We added damping and smoothing constraints to regularize the inversion. The tomographic model is given by the solution of 𝐽 𝑇 𝑊 𝑜𝑏𝑠 𝛿𝛥𝑡 = (𝐽 𝑇 𝑊 𝑜𝑏𝑠 𝐽 + 𝑊 𝑑𝑎𝑚𝑝 + 𝑊 𝑠𝑚𝑜𝑜𝑡ℎ )𝛥𝑚, Equation 1J o u r n a l P r e -p r o o f where 𝛿𝛥𝑡 is the vector of relative traveltimes, 𝐽 is the Jacobian matrix, 𝑊 𝑜𝑏𝑠 is the weight matrix of the observations, 𝑊 𝑑𝑎𝑚𝑝 is the damping matrix, 𝑊 𝑠𝑚𝑜𝑜𝑡ℎ is the smoothing matrix and 𝛥𝑚 is the vector of velocity perturbations. The Jacobian matrix for each event is given by 𝐽 𝑖𝑗 = (𝐺 𝑖𝑗 -𝐺 ̅ 𝑖𝑗 ), Equation 2where each element of 𝐺 𝑖𝑗 is the travel time of the ith ray in the jth cell, and the elements of 𝐺 ̅ 𝑖𝑗 are the averages of 𝐺 𝑖𝑗 over all the stations that have recorded the event. Each element of the Jacobian is the travel time of the ith ray in the jth cell from which the mean travel time of all rays of the event passing through the cell is subtracted. The concept behind this modification is that 𝐺 ̅ 𝑖𝑗 𝛥𝑚 is equivalent to the mean travel time of each event, therefore, the modified Jacobian matrix has the effect of removing the average ray contribution. This eliminates the necessity for arbitrarily choosing one mean velocity as a reference. Recently, [START_REF] Maupin | Combining asynchronous data sets in regional body-wave tomography[END_REF] has shown that using demeaned event sensitivity kernels also removes the long wavelength biases arising from asynchronous station deployments. Therefore, our inversion scheme is well adapted to the joint inversion of travel time data coming from different permanent and temporary sub-arrays that may not have recorded the same number of events.

We regularize the inversion by adding damping and smoothing constraints. To smooth the solution, we damp the norm of the Laplacian of the model. The procedure for finding the optimum damping and smoothing coefficients is as follows. 1) We start the tomography with two initial damping and smoothing factors. 2) We run the inversion code iteratively inside two loops.

The outer loop consists of 10 iterations over the damping factor, and the inner loop consists of 20 iterations over the smoothing factor. In this step 10 L-curves, one for each damping factor are J o u r n a l P r e -p r o o f produced (see Figure S1 for selected L-curve). 3) We repeat steps 1 and 2 four times with different initial factors and increments to ensure that the best initial values are found. We pick the best L-curve from a set of 40. To choose the best damping and smoothing factors, we impose a RMS reduction of at least 50 percent, and a smooth final model not dominated by small-scale anomalies. The inversion achieved a final RMS reduction of 56.8%. The LSQR (Least Squares using the QR factorization) method of [START_REF] Paige | LSQR: An algorithm for sparse linear equations and sparse least squares[END_REF] was used to solve the equations. This method is very fast and employs sparse matrix techniques.

We performed a synthetic checkerboard test to identify the parts of the model that are well resolved and to determine the amount of vertical smearing in regions of low ray density. In this test, the same complete ray distribution as in our real dataset has been used. The starting synthetic velocity model is made of 75-km cubic velocity anomalies organized in 3 horizontal layers in the shallow part of the model (< 450 km) and 2 layers of 100-km cubic anomalies in the deeper parts. The input velocity anomalies had ±5% perturbations. The synthetic blocks were placed at 50-125 km, 200-275 km, 350-425 km, 500-600 km, and 675-775 km depth. We also added 10% random noise to the synthetic relative residuals. The synthetic test achieved an RMS reduction of 62%. Fig. 3 shows the results of the test in layers 4, 10, 16, 23, and 30 of the model.

In the shallow parts of the model (i.e., layer 4) the velocity variations are best resolved in regions of good station coverage. In the mid-upper mantle depths (i.e., layers 10 and 16) where ray crossing reaches its maximum, good resolving power is achieved in many parts of the model. In the deeper parts (i.e., layers 23 and 30), the ray density is distributed more evenly among the cells, and thus the inversion returns a more uniform pattern of velocity anomalies. Fig. 4 shows the results of the synthetic test along 7 cross-sections from northwest to southeast. In general, the resolution is good in the upper 400 km. In the regions far from the temporary arrays (e.g., J o u r n a l P r e -p r o o f sections d, e, and f) the near-surface anomalies show somewhat lower amplitudes with respect to the deeper features. Beneath the CIGSIP array and neighboring regions (section b), the resolution is good from the top to the bottom of the model. However, vertical smearing in the lower half of the model becomes significant in the northwestern part of the model (sections a, b and c). The first anomaly, marked as ALH (Arabian Lithosphere High) is a shallow high-velocity feature beneath the Zagros and SSZ extending down to a depth of about 270-300 km. It also extends beneath central Iran by as much as 75 km north of the MZT (see also Fig. 5d). In the eastern half of the model, the anomaly is poorly resolved at shallower depth where station coverage is sparse (e.g., Fig. 7b andc). We interpret the ALH anomaly as the Arabian margin underthrusted beneath central Iran. The second anomaly marked as TSH (Tethyan Slab High), is a highvelocity region located north of the MZT under SSZ and central Iran. In most of the cross sections it is located below 400 km depth and extends to the bottom of the upper mantle, but in NW, as well as SE (Fig. 6k andl) it penetrates into the lower mantle as well. The resolution test (e.g., Fig. 4b-e) indicates that these features in the upper mantle are well resolved. We interpret the high-velocity TSH as the subducted Neotethyan slab. The third anomaly, marked as CIH (Central Iran High), is a high-velocity body beneath central Iran, observed at 200 km depth in the western part of the model. Towards the eastern regions and beneath central and eastern Alborz it gradually deepens to about 450 km depth (Fig. 7c). In the northern end of the profiles crossing the Alborz (e.g., Fig. 6a andb), a high-velocity shallow zone characterizes the transition between the Alborz Mountains and the South Caspian Basin, which suggests that the nature of the

Tomographic Results

Fig

Discussion

Seismic velocities in the Earth are influenced by temperature, anelasticity, partial melting, water content, and rock composition. Although the effect of anelasticity and rock composition can be largely ignored in the upper mantle where the oceanic slab is descending [START_REF] Goes | Shallow mantle temperatures under Europe from P and S wave tomography[END_REF], isolating their relative contribution to the velocity field is challenging. In the following, we therefore assume that velocity perturbations mainly result from temperature anomalies.

Several scenarios have been proposed to explain the geodynamic evolution of the Arabia-Eurasia plate system. Among them, models suggested by [START_REF] Dercourt | Atlas Tethys, paleoenvironmental maps[END_REF] and [START_REF] Şengör | Paleotectonics of Asia: fragments of a synthesis[END_REF] and another one by Stampfli andBorel (2002, 2004) have found wider acceptance. In the first model, except for a short-lived intra-oceanic subduction near the southern margin of the Neo-Tethys that resulted in the emplacement of the Oman ophiolites on the Arabian margin, all of the subduction has taken place at the northern margin of Neo-Tethys beneath Anatolia and central Iran. In the model proposed by [START_REF] Stampfli | The TRANSMED transects in space and time: constraints on the paleotectonic evolution of the Mediterranean domain[END_REF] the subduction of the Neo-Tethys ridge beneath the Iranian blocks (Lut and SSZ) in the Mesozoic was followed by the opening of an intra-oceanic back-arc basin east of Arabia (the so-called Semail Ocean). During the Cenozoic, the Semail Ocean subducted beneath Eurasia. [START_REF] Hafkenscheid | Subduction history of the Tethyan region derived from seismic tomography and tectonic reconstructions[END_REF] According to this reconstruction, the high-velocity bodies presently observed in the mantle transition zone are mainly the Arabian margin lithosphere that has been pushed down since the Oligocene break-off and during the course of continental collision. [START_REF] Agard | Zagros orogeny: a subduction-dominated process[END_REF] presented a synthesis of various geological data to reconstruct the deformation history in the Zagros collision zone from 150 Ma to present. Their reconstruction is broadly similar to that of [START_REF] Hafkenscheid | Subduction history of the Tethyan region derived from seismic tomography and tectonic reconstructions[END_REF]. They outlined a major slab break-off episode in Late Paleocene, and a more recent partial slab tear during the Miocene (~ 10 Ma) based on adakite formation in parts of the UDMA in central Iran [START_REF] Omrani | Arcmagmatism and subduction history beneath the Zagros Mountains, Iran: a new report of adakites and geodynamic consequences[END_REF].

The structure of the Zagros lithosphere

Our model shows high velocities in the Zagros lithosphere and low velocities north of MZT in central Iran. A similar structure has been mapped in previous studies (e.g., [START_REF] Asudeh | Pn velocities beneath Iran[END_REF][START_REF] Bijwaard | Closing the gap between regional and global travel time tomography[END_REF][START_REF] Kaviani | A strong seismic velocity contrast in the shallow mantle across the Zagros collision zone (Iran)[END_REF][START_REF] Shomali | Lithospheric structure beneath the Zagros collision zone resolved by non-linear teleseismic tomography[END_REF][START_REF] Koulakov | High-frequency P and S velocity anomalies in the upper mantle beneath Asia from inversion of worldwide traveltime data[END_REF][START_REF] Priestley | The Zagros core: deformation of the continental lithospheric mantle[END_REF][START_REF] Mahmoodabadi | Upper-mantle velocity structure beneath the Zagros collision zone, Central Iran and Alborz from nonlinear teleseismic tomography[END_REF]. We interpret these anomalies as evidence for a colder and thicker lithosphere beneath the Zagros versus a thinner and warmer lithosphere in the surrounding regions. [START_REF] Priestley | The Zagros core: deformation of the continental lithospheric mantle[END_REF] used temperature modeling to derive the depth of the lithosphere-asthenosphere boundary from shear wave velocity profiles. They obtained a value of ~120 km for the thickness of the lithosphere everywhere except in the Zagros, for which their estimate was twice as much. [START_REF] Priestley | The Zagros core: deformation of the continental lithospheric mantle[END_REF] with our cross-section in Fig. 6c. We also identify a very thick lithosphere in the Zagros (anomaly ALH), but our model does not support a detachment of the Arabian lithosphere. [START_REF] Mahmoodabadi | Shear wave velocity structure of the upper-mantle beneath the northern Zagros collision zone revealed by nonlinear teleseismic tomography and Bayesian Monte-Carlo joint inversion of surface wave dispersion and teleseismic P-wave coda[END_REF] also resolved a shallow low-velocity feature in NW Zagros south of the SSZ, and explained it as a narrow low-velocity corridor formed as a consequence of delamination of the Arabian lithosphere. In general, we observe a thick lithosphere along the strike of the Zagros throughout the range (e.g., Fig. 5b,c). However, important variations in the geometry of the frontal edge of the advancing Arabian plate are observed from NW to SE. In NW Zagros, the Arabian lithosphere achieves a maximum of 120 km of underthrusting beneath the SSZ at the location of the Zagros03 seismic line (Fig. 5d). In contrast, in central Zagros They argued that the abnormal thickness of the lithosphere could be the result of internal deformation of the Arabian lithosphere as it collides against the relatively strong central Iran lithosphere. Several previous tomographic studies have reported a similar non-cylindrical structure of the Zagros lithosphere along the collision zone. The P and S-wave imaging of the upper mantle in Iran by [START_REF] Alinaghi | Seismic tomographic imaging of P-and S-waves velocity perturbations in the upper mantle beneath Iran[END_REF] and [START_REF] Koulakov | High-frequency P and S velocity anomalies in the upper mantle beneath Asia from inversion of worldwide traveltime data[END_REF] show that central Zagros is underlain by a low velocity anomaly in the shallow upper mantle that is not observed beneath NW and SE Zagros, in good agreement with [START_REF] Chang | Joint inversion for threedimensional S velocity mantle structure along the Tethyan margin[END_REF] and [START_REF] Priestley | The Zagros core: deformation of the continental lithospheric mantle[END_REF].

The structure of the subducted slab

A prominent feature of our tomographic model is the combined ALH + TSH high-velocity region beneath the collision zone that extends from the surface down to the transition zone. In the southeastern part of the study region, it penetrates the lower mantle. The apparent dip of the structure is about 60 o -70 o towards northeast throughout the collision front. We interpret this structure as the underthrusting Arabian margin, still partially connected to the last fragments of the subducted oceanic slab, with a detachment around 200-300 km depth. The degree of detachment varies along the collision zone; whereas slab detachment seems to have remained at an early stage beneath NW Zagros, it seems more mature beneath central Zagros. The depth of detachment is just beneath the thickened Arabian lithosphere, and since the top of the detached slab is still observed at a shallow level, it might be a relatively recent event developed in the latter stages of collision. [START_REF] Omrani | Arcmagmatism and subduction history beneath the Zagros Mountains, Iran: a new report of adakites and geodynamic consequences[END_REF] documented a post-Late Miocene formation of

J o u r n a l P r e -p r o o f

Journal Pre-proof volcanic rocks with adakitic signature in UDMA in central Iran, and based on their spatially limited distribution along the arc, they inferred that they must have formed as a result of slab melting during a slab break-off event. [START_REF] Agard | Zagros orogeny: a subduction-dominated process[END_REF] expanded on this idea and used other geological evidence to suggest that a late partial break-off/tear caused the detachment of the Neo-Tethyan slab from the continental lithosphere at around 10 Ma. They also argued that the extent of break-off must be greater in central Zagros than in the northwest. A number of previous tomographic studies have provided evidence for slab detachment in the Zagros. [START_REF] Van Der Meer | Atlas of the underworld: Slab remnants in the mantle, their sinking history, and a new outlook on lower mantle viscosity[END_REF] showed a tomographic image beneath northern Zagros where a detached slab segment with its top just below the Zagros root descends into the lower mantle. In western Zagros, Mahmoodabadi et al. ( 2019) evidenced a high-velocity anomaly at 300 km depth disconnected from the Zagros lithosphere. They interpreted it as a remnant of the Neo-Tethyan oceanic slab detached from the leading continental edge. [START_REF] Rahmani | Deep velocity image of the north Zagros collision zone (Iran) from regional and teleseismic tomography[END_REF] in their tomography study on the CIGSIP array obtained a similar image. They found a pronounced high-velocity anomaly at depths greater than 350 km north of the Zagros suture and interpreted it as a remnant of the Neo-Tethyan slab. However, their tomographic model does not show a separation between the continental lithosphere at the surface and the remnant slab at depth. They performed a synthetic modeling and showed that the apparent continuity between the two high-velocity regions could be the result of smearing. Nevertheless, they stopped short of explicitly arguing for a slab detachment event.

The available evidence from seismic anisotropy is in good agreement with a slab detachment event beneath the Zagros, as well as to the generally-complex lithospheric thickness variations between the Zagros and central Iran as discussed in the previous section. Recent shear wave splitting studies by [START_REF] Sadeghi-Bagherabadi | Seismic anisotropy and its geodynamic implications in Iran, the easternmost part of the Tethyan Belt[END_REF], [START_REF] Kaviani | Mantle-flow diversion beneath the Iranian plateau induced by Zagros' lithospheric keel[END_REF] evidenced slab detachments (i.e. profiles CR6 to CR12, and CR37 and region to its SE in Fig. 6 and7). The evidence from seismic anisotropy corroborates our evidence for slab detachment and furthermore, it indicates that the zone of detachment possibly has a significant role in the organization of the mantle flow field beneath the collision zone.

Most cross sections in NW Zagros image a Tethyan slab descending into the lower mantle directly. This feature is not as clear in the southern part of the model (e.g. profile CR37 in Fig. 7). [START_REF] Kaviani | Mantle transition zone thickness beneath the Middle East: Evidence for segmented Tethyan slabs, delaminated lithosphere, and lower mantle upwelling[END_REF] studied the variations of the mantle transition zone thickness beneath the Middle East. They identified a significant depression of the 660 boundary and a modest uplift of the 410 in NW Zagros, corresponding to our western profiles. They interpreted this feature as J o u r n a l P r e -p r o o f Journal Pre-proof resulting from the presence of a slab segment in western Zagros. In south central Zagros they observed the uplift of both 410 and 660 discontinuities, which precludes the presence of stagnant slab in that part of the mantle.

To evaluate the reliability of our tomographic model and assess the amount of smearing that may mask the distances between the prominent anomalous regions, we performed additional synthetic resolution experiments to test if the major anomalies observed in our tomographic model are adequately resolved. Specifically, we modelled the ALH, TSH, and CIH anomalies in the region encompassing the CIGSIP and Zagros03 arrays where there is good station coverage and raycrossing density. Fig. 8 shows the results of the synthetic test. The real-data anomalies were modelled by uniform horizontal prisms placed across the western part of the model. The ALH anomaly had a polygonal shape with its base at 200 km depth. The TSH and CIH had rectangular sections. We varied the position of the top of the TSH prism, from 25 km to 175 km below the base of ALH. Fig. 8 illustrates the extent of vertical smearing in our inversions. In the 200-400 depth range, where ray crossing is maximum, horizontal boundaries have been displaced by as much as 50 km along the vertical direction. Vertical boundaries are much better resolved. In one model run we assumed an unbroken slab down to 650 km depth (Fig. 8d). The inversion returned a continuous unbroken high-velocity region from surface to bottom, a very different pattern from what is observed in the real-data tomography. Anomalous bodies vertically separated by more than 100 km can be resolved. We also observe that the geometry of the resolved synthetic anomalies changes by only a moderate amount between the various cross-sections. Overall, the structures beneath profiles CR 9 to 14 are the best resolved.

van Hunen and Allen (2011) quantitatively investigated slab break-off subsequent to continental subduction. They found that the mechanical strength of the slab has a major role in the timing of J o u r n a l P r e -p r o o f Journal Pre-proof slab break-off after the onset of collision. For an old and strong subducting slab this time is 20-25 Ma, while for a young and weak slab it can be as early as 10 Ma. [START_REF] Van Hunen | Continental collision and slab break-off: A comparison of 3-D numerical models with observations[END_REF] applied their findings to the Arabia-Eurasia collision zone by assuming an age of 35 Ma for the initial collision and a Permian age of 200 Ma for the subducting Neo-Tethyan slab. They deduced that slab break-off could not have occurred before 15-10 Ma. Better estimates now put the onset of collision in the Zagros at around 27-25 Ma [START_REF] Egan | Subsidence and uplift mechanisms within the South Caspian Basin: insights from the onshore and offshore Azerbaijan region[END_REF]Agard et al., 2011 and references therein). Following van Hunen and Allen (2011), slab detachment in the Zagros could have occurred as recently as 5 Ma ago, significantly more recently than what the geological evidence suggests. This discrepancy can be resolved by assuming that the subducting slab was weaker than what is expected for its age. One can also argue that the actual slab break-off in the Zagros has been a faster process. In fact, numerical modeling by [START_REF] Duretz | Numerical modelling of spontaneous slab breakoff and subsequent topographic response[END_REF][START_REF] Duretz | Dynamics of slab detachment[END_REF] incorporating a combination of viscous and Peierls creep mechanisms predicted that an 80 Ma old oceanic lithosphere subducting at a rate of 5 cm/yr can break off significantly sooner, at about 11 Ma after the onset of collision and at a depth of about 300 km. In terms of time and depth of break-off, the predictions by [START_REF] Duretz | Numerical modelling of spontaneous slab breakoff and subsequent topographic response[END_REF][START_REF] Duretz | Dynamics of slab detachment[END_REF] are in good agreement with the available geological evidence and our tomography results in the Zagros. Since our tomographic images evidence a break off developed unevenly along the collision zone, this would suggest that the process is still ongoing.

The upper mantle structure beneath central Iran and the Alborz

In our model we observe slow velocities in the shallow upper mantle (<175 km) beneath central Iran and the Alborz (anomaly CIL). Voluminous high-velocity bodies are observed between 200 to 350 km depth (anomaly CIH), beneath the southern flank of the Alborz Mountains and adjacent regions in central Iran (e.g., Fig. 5e).

J o u r n a l P r e -p r o o f

Journal Pre-proof

Previous studies have reported slow propagation of body and surface waves and high attenuation rates across the Turkish-Iranian Plateau (e.g., [START_REF] Al-Lazki | Pn tomographic imaging of mantle lid velocity and anisotropy at the junction of the Arabian, J o u r n a l P r e -p r o o f Journal Pre-proof Eurasian and African plates[END_REF][START_REF] Maggi | Surface waveform tomography of the Turkish-Iranian plateau[END_REF][START_REF] Al-Damegh | Regional seismic wave propagation (Lg and Sn) and Pn attenuation in the Arabian Plate and surrounding regions[END_REF][START_REF] Kaviani | A strong seismic velocity contrast in the shallow mantle across the Zagros collision zone (Iran)[END_REF][START_REF] Mahmoodabadi | Upper-mantle velocity structure beneath the Zagros collision zone, Central Iran and Alborz from nonlinear teleseismic tomography[END_REF]. These observations have led many authors to suggest that central Iran has a thin lithosphere and a warm upper mantle. The LAB depth beneath the southern side of the Alborz Mountains calculated by [START_REF] Rastgoo | Deep structure of the Alborz Mountains by joint inversion of P receiver functions and dispersion curves[END_REF] from receiver functions and dispersion curves is around 150 km. Central

Iran is an elevated plateau with a surface relief of 1000-1500 m built during the course of the collision. Despite its relatively high topography, it is only moderately thickened compared to the Zagros. The measured Moho depths from receiver functions in central Iran range from 42 km north of the Zagros suture [START_REF] Paul | Seismological evidence for crustal-scale thrusting in the Zagros mountain belt (Iran)[END_REF] to 48 km just south of the Alborz Mountains [START_REF] Radjaee | Variation of Moho depth in the central part of the Alborz Mountains, northern Iran[END_REF]. [START_REF] Hatzfeld | Comparisons of the kinematics and deep structures of the Zagros and Himalaya and of the Iranian and Tibetan plateaus and geodynamic implications[END_REF], through a detailed analysis of available geologic, present-day kinematics of deformation, and deep structure data, discussed the possibility of convective removal of the mantle lid in Tibet and in central Iran. They argued that the low seismic velocities in the uppermost mantle just below the Moho depth beneath both plateaus are consistent with Late-Cenozoic removal of the mantle lithosphere. Several recent seismic tomographies in Tibet have produced strong evidence for lithospheric foundering/delamination as a result of continental collision [START_REF] Ren | Finite frequency tomography in southeastern Tibet: evidence for the causal relationship between mantle lithosphere delamination and the north-south trending rifts[END_REF][START_REF] Chen | Lithospheric foundering and underthrusting imaged beneath Tibet[END_REF][START_REF] Huang | P and S wave tomography beneath the SE Tibetan Plateau: Evidence for lithospheric delamination[END_REF]. The seismic structure of these models consists of a shallow lowvelocity zone in the uppermost mantle sitting over a deeper high-velocity body that extends as deep as the mantle transition zone. In the central Iranian Plateau, the seismic structure of the upper mantle is much less clear than in Tibet and [START_REF] Hatzfeld | Comparisons of the kinematics and deep structures of the Zagros and Himalaya and of the Iranian and Tibetan plateaus and geodynamic implications[END_REF] were not able to present unequivocal evidence from seismic tomography for convective removal of the mantle lithosphere. Seismological studies remain so far inconclusive on the hypothesis of mantle J o u r n a l P r e -p r o o f Journal Pre-proof delamination in central Iran. The low-resolution surface wave study by [START_REF] Maggi | Surface waveform tomography of the Turkish-Iranian plateau[END_REF] in western Asia mapped a low-velocity upper mantle beneath central Iran which was interpreted as an indicator of a thin delaminated lithosphere and a warm shallow mantle. The regional P wave tomography by [START_REF] Mahmoodabadi | Upper-mantle velocity structure beneath the Zagros collision zone, Central Iran and Alborz from nonlinear teleseismic tomography[END_REF] has revealed a similar low-velocity shallow structure, but the authors invoked the subduction model by [START_REF] Verdel | A Paleogene extensional arc flare up in Iran[END_REF] to explain their observations. According to [START_REF] Verdel | A Paleogene extensional arc flare up in Iran[END_REF], flat subduction beneath the Eurasian margin in the Cretaceous was followed by slab roll-back, extension, and magmatism in central Iran and asthenospheric upwelling during the Eocene, resulting in a pre-collisional thin lithosphere and warm upper mantle on the Eurasian side. Petrological evidence for postcollisional delamination in northern Iran comes from the region of Damavand, a Quaternary volcano in central Alborz. [START_REF] Davidson | The geology of Damavand volcano, Alborz Mountains, northern Iran[END_REF] and [START_REF] Liotard | Origin of the absarokite-banakite association of the Damavand volcano (Iran): Trace elements and Sr, Nd, Pb isotope constraints[END_REF] argued that the geochemistry of the Damavand magmatism points to lithospheric delamination in the Alborz. [START_REF] Shabanian | Structural control on volcanism in intraplate post collisional settings: Late Cenozoic to Quaternary examples of Iran and Eastern Turkey[END_REF] analyzed structural controls on the Alborz volcanism in a post-collisional setting and proposed that the volcanic activity of Damavand has been the result of magma generation due to a rising asthenosphere following lithospheric delamination in a transtensional environment. [START_REF] Rastgoo | Deep structure of the Alborz Mountains by joint inversion of P receiver functions and dispersion curves[END_REF] used their 2-D seismic profile along the southern flank of the Alborz Mountains to verify the structural model proposed by [START_REF] Shabanian | Structural control on volcanism in intraplate post collisional settings: Late Cenozoic to Quaternary examples of Iran and Eastern Turkey[END_REF]. In western and central Alborz their model revealed a low-velocity layer in the 50-100 km depth range and a higher-velocity zone at 100-150 km depth. The spatial distribution of these velocity zones was interpreted as the result of post-collisional delamination of the lower part of the western Alborz lithosphere.

The previous seismological studies proposing delamination in central Iran, have only provided evidence for a thin "delaminated" lithosphere, without revealing anything that can be taken as J o u r n a l P r e -p r o o f proof for the presence of the sunken root itself. In several of our profiles and depth plots, the high-velocity anomaly CIH can be spotted at different depths between 200 and 350 km. It is the best observed in the northern margin of central Iran and the southern side of the Alborz range.

The CIH anomaly progressively deepens from northwest to southeast (compare successive profiles in Fig. 6 and7), suggesting a foundering body at different stages of descent. We suggest that the CIH body might offer tangible evidence for the sunken root. In terms of depth interval and vertical configuration, the CIH anomaly is similar to what seismic studies mentioned above have described in Tibet. [START_REF] François | Buildup of a dynamically supported orogenic plateau: Numerical modeling of the Zagros/Central Iran case study[END_REF] through numerical modeling of continental subduction, attempted to test the mechanisms for the buildup of dynamic topography in the central Iranian Plateau. They showed that a recent slab break-off (5-10 Ma) after the onset of collision, could have indirectly initiated a mantle flow field beneath the overriding plate, resulting in the delamination of the upper plate lithospheric mantle. A late-stage isostatic readjustment of the uplifted plateau would follow without undergoing significant crustal thickening. Their successful models were able to produce small-scale convective instabilities beneath the upper plate with wavelengths comparable to the size of the fast anomalies we have observed in our tomographic model. Effective detachment of an oceanic slab from continental lithosphere is often followed by the upwelling of the asthenosphere through slab gaps, eduction, slab retreat, and eventually volcanism and exhumation of ultra-high-pressure rocks, which are not currently observed in central Iran. On the other hand, the lack of deep seismicity, especially deep extensional events in the Zagros, indicates that the subducted slab does not exert gravitational pull on the base of the continental plate. To reconcile these apparently contradictory observations, we propose that slab break-off in NW Zagros is very young.

Conclusions

Table1: Summary of network and data information used in this study

Network name

No

  The temporary stations belong to several arrays deployed at different times in various parts of the collision zone. The Zagros01 and Zagros03 arrays consist of two transects deployed in central and western Zagros for a few months in 2001 and 2003 respectively. The CIGSIP array operated for one year (2013-2014) along three parallel lines in western Iran covering the Zagros and Alborz mountains. The IASBS-CAM-1 (2008-2012) and IASBS-CAM-2 (2014-2016) arrays were installed in NW Iran and western Alborz, respectively, and operated for 6 to 30 months. We also use the data of the permanent stations of the Iranian national networks IRSC and INSN recorded between 2012 and 2019 to increase the resolution of our output model in the regions between the temporary arrays. With this distribution of seismic stations our tomographic study covers a significant part of the collision zone: NW and central Zagros, western and central Alborz, and a significant part of central Iran.

  Fig. 2b shows the average station residuals. Stations in the Zagros Mountains and the Sanandaj-

Fig. 5

 5 Fig.5shows 12 horizontal plots of the velocity anomalies between 25 and 825 km depth. At

  J o u r n a l P r e -p r o o fCaspian lithosphere is different from the continental regions surrounding it. The fourth major anomalous zone in the model is the shallow low-velocity anomaly CIL (Central Iran Low) beneath central Iran and the Alborz. This anomaly is confined to the upper 150-200 km and reveals a thinner central Iran lithosphere bounded by the thick Zagros and South Caspian lithospheres.

  used the above models to predict the present positions of Neo-Tethyan lithospheric domains subducted J o u r n a l P r e -p r o o f since the Mesozoic. They compared their predictions with results from global mantle P-wave tomographies and concluded that most of the subducted Neo-Tethyan lithosphere is presently residing in the lower mantle between 800 and 2000 km depth. They further argued for the occurrence of an Early Oligocene slab detachment under the northern Zagros suture zone.

  argued that the thickening of the Zagros lithosphere results from shortening.Motaghi et al. (2017a,b) inverted receiver functions from J o u r n a l P r e -p r o o f seismic profiles in western and southern Zagros to derive shear wave velocity profiles, from which they inferred a deep LAB and thick lithosphere beneath the Zagros and SSZ (> 200 km).They interpreted this thickening as resulting from the underthrusting of the Arabian lithosphere beneath central Iran.[START_REF] Mohammadi | New constraints on lithospheric thickness of the Iranian plateau using converted waves[END_REF] used P-to-S and S-to-P converted phases to constrain LAB depths of 130 and 150 km beneath western Zagros and SSZ, respectively, and significantly shallower depths of 80-90 km in central Iran and the Alborz. They also mapped a sharp lithospheric thickness variation between the SSZ and UDMA in central Iran. The P-wave tomography study by Mahmoodabadi et al. (2019) in western Zagros (coinciding with the region of Zagros03 array) also found a thick high-velocity lithospheric mantle beneath the Zagros and low-velocities beneath central Iran and Alborz. The 2-D P-wave tomography of Rahmani et al. (2019) along the CIGSIP profile in NW Zagros showed a low-velocity wedge at a depth of 150 km near the frontal edge of the Arabian plate beneath the suture zone. They interpreted it as the beginning of the detachment of the lower part of the Arabian lithosphere from its upper part. The profile of Rahmani et al. (2019) coincides

  J o u r n a l P r e -p r o o f (southeast of Zagros03 line in Fig 5) our model resolves low-amplitude velocity anomalies on both sides of the suture zone, with no clear evidence for significant underthrusting beneath central Iran. In SE Zagros, Motaghi et al. (2017b) observed evidence for lithospheric buckling with no indication of significant underthrusting of the Arabian lithosphere beneath central Iran.

  and Arvin et al. J o u r n a l P r e -p r o o f (2021) have revealed a complex mantle flow pattern beneath Iran compared to a simpler flow regime beneath the Anatolian Plateau and Arabian Plate, where lithospheric thickness varies relatively smoothly. Kaviani et al. (2021) hypothesized that the southwest mantle flow field is deflected beneath the Zagros by the lithospheric keel, producing a NW-SE toroidal flow beneath central Iran parallel to the strike of the collision zone. Similarly, Sadeghi-Bagherabadi et al. (2018) invoked a combination of the Zagros keel effect and an asthenospheric edge-driven convective flow beneath central Iran to explain the NW-SE fast splitting directions under the southwestern margin of central Iran. In NW Iran and the Alborz region where the lithosphere is thinner, the anisotropic fast-axis directions suggest a simple mantle flow field away from the collision zone. The inferred circular pattern around the Zagros keel implies sharp boundaries for the Zagros lithosphere, in good agreement with our tomography results. It also implies a gap in the slab to explain the mantle flow transfer to central Iran. This gap may have result from the slab detachment event. The proposed keel deflection and associated circular flow is mainly concentrated in NW and South central Zagros, coinciding with locations where we have

Figure 9

 9 Figure 9 summarizes our findings of the lithospheric and upper mantle structure beneath the

Figure 1 .

 1 Figure 1. Map showing the seismic networks used in this study as well as the major tectonic units: Zagros Fold and Thrust Belt (ZFTB), Main Zagros Thrust (MZT), Alborz Mountains, Sanandaj-Sirjan Zone (SSZ), Urmia-Dokhtar Magmatic Arc (UDMA), the Central Iran Microblock, and Makran Subduction Zone. The positions of the crosssections, numbered from 1 to 40, are also shown with black solid lines. J o u r n a l P r e -p r o o f

Figure 2 .Figure 3 .

 23 Figure 2. (a) Map showing the epicenters of teleseismic events used in this study. Epicentral distances are

Figure 4 .

 4 Figure 4. Results of the spike test inversion along 7 cross-sections. Location and surface

Acknowledgments

This research was funded by several IASBS internal grants. M.V. benefited from support provided by the Ministry of Science, Research and Technology of Iran and Laboratoire Géosciences Environnement Toulouse -GET. We acknowledge ISC and INSN for use of their permanent station data. The Zagros01 and Zagros03 arrays were installed by LGIT, Grenoble and IIEES, Tehran. The CIGSIP array was supported by the Strategic Priority Research Program (B) (Grant no. XDB03010802) and the International Partnership Program (GIH21776) of the Chinese Academy of Sciences. The IASBS-CAM arrays were jointly operated by IASBS and J o u r n a l P r e -p r o o f Journal Pre-proof

Univ. of Cambridge. Maps were prepared using the Generic Mapping Tools (Wessel & Smith 1998). We thank associate editor Ling Chen, Zhouchuan Huang, and two other anonymous reviewers for their valuable comments and critique of the manuscript.

 [START_REF] Karasözen | Seismotectonics of the Zagros (Iran) from orogen-wide, calibrated earthquake relocations[END_REF].

J o u r n a l P r e -p r o o f

Journal Pre-proof

J o u r n a l P r e -p r o o f

Journal Pre-proof Declaration of interests ☒ 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. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

J o u r n a l P r e -p r o o f

Journal Pre-proof

We map the lateral variations of P-wave velocity beneath the Zagros collision zone in Iran down to 800 km depth by regional travel time tomography. Our tomographic study documents variations of lithospheric thickness between the Zagros and central Iran and a post-collisional slab detachment around 250 km depth. If we assume that the detachment occurred right beneath the base of the subducted Arabian continental margin, then the maximum thickness of the continental root must be no greater than 225 km. This would imply a doubling of the Zagros lithosphere as a result of collision. The absence of high-velocity anomalies far north from the MZT at the lithospheric depth range indicates that the lithosphere thickening has not propagated into the interior of central Iran. The high-velocity anomalies at 200-300 km depth beneath central Iran suggest the presence of fragments of delaminated mantle lithosphere. They also suggest post-collisional slab break-off in central Zagros that seems to be at an earlier stage in NW Zagros compared to central Zagros. From the shallow depth of the top of the detached slab, we infer that the break-off is recent

(5-10 Ma), in good agreement with geological records. Effective detachment of an oceanic slab from continental lithosphere is often followed by the upwelling of the asthenosphere through slab windows, eduction, slab retreat, and eventually volcanism and exhumation of ultra-high-pressure rocks, which are not currently observed in central Iran. On the other hand, lack of deep seismicity, especially deep extensional events in the Zagros, indicates that the subducted slab does not exert gravitational pull on the base of the continental plate to reconcile these apparently contradictory observations, we propose that slab break-off in NW Zagros is very young. • Evidence for a thicker lithosphere beneath the Zagros versus thinner lithosphere in the surrounding regions.

• Slab detachment near the base of the thickened Zagros lithosphere.

• Seismological evidence for mantle lithospheric delamination beneath central Iran.

J o u r n a l P r e -p r o o f