Nadaya Cubas
email: nadaya.cubas@sorbonne-universite.fr

Philippe Agard

Roxane Tissandier

Earthquake ruptures and topography of the Chilean margin controlled by plate interface deformation

What controls the location and segmentation of mega-earthquakes in subduction zones is a long-standing problem in earth sciences. Prediction of earthquake rupture extent mostly relies on interplate coupling models based on Global Navigation Satellite Systems providing patterns of slip deficit between tectonic plates. We here investigate if and how the strongly and weakly coupled patches revealed by these models relate to the distribution of deformation along the plate interface, i.e. basal erosion and/or underplating. From a mechanical analysis of the topography applied along the Chilean subduction zone, we show that extensive plate interface deformation takes place along most of the margin. We show that basal erosion occurs preferentially at 15 km depth while underplating does at 35 ± 10 and 60 ± 5 km depth, in agreement with P-T conditions of recovered underplated material, expected pore pressures, and spatial distribution of marine terraces and uplift rates. South of the Juan Fernandez Ridge, large sediment input favors shallow accretion and underplating of subducted sediments, while along northern Chile, extensive basal erosion provides material for the underplating. We then show that, along the accretionary margin, the two last major earthquakes are limited along their down-dip end by underplating while, along the erosive margin, they are surrounded by both basal erosion and underplating. Segments with heterogeneously distributed deformation largely coincide with lateral earthquake terminations. We therefore propose that long-lived plate interface deformation promotes stress build-up and leads to earthquake nucleation. Earthquakes then propagate along fault planes shielded from this long-lived permanent deformation, and are finally stopped by segments of heterogeneously distributed deformation. Slip deficit patterns and earthquake segmentation therefore reflect the along-dip and along-strike distribution of the plate interface deformation.

Topography acts as a mirror of distributed plate interface deformation and should be more systematically studied to improve the prediction of earthquake ruptures.

Introduction

Predicting the extent of subduction zone seismic ruptures mostly relies on geodetic observations of interseismic coupling (e.g., [START_REF] Moreno | Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone[END_REF], notwithstanding spatial resolution issues [START_REF] Loveless | Spatial correlation of interseismic coupling and coseismic rupture extent of the 2011 Mw= 9.0 Tohoku-oki earthquake[END_REF]. Some earthquakes have however ruptured low coupled zones [START_REF] Noda | Stable creeping fault segments can become destructive as a result of dynamic weakening[END_REF] or only part of highly locked patches (e.g., [START_REF] Konca | Partial rupture of a locked patch of the Sumatra megathrust during the 2007 earthquake sequence[END_REF]. This was the case for the 2014 Mw 8.4 Iquique earthquake, which only ruptured the northern end of the North Chile gap (last broken in 1877, Ruiz et al., 2014a), or the Tocopilla Mw 7.8 earthquake, which struck the down-dip and less coupled part of the same 1 gap seven years before [START_REF] Delouis | Recent crustal deformation in the Antofagasta region (northern Chile) and the subduction process[END_REF][START_REF] Metois | Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38-18 S)[END_REF] (Fig. 1). Predicting the location and extent of earthquakes clearly requires a better understanding of earthquake mechanics.

Subduction earthquake propagation has so far been either related to megathrust frictional and pore fluid pressure properties [START_REF] Perfettini | Seismic and aseismic slip on the Central Peru megathrust[END_REF][START_REF] Kaneko | Towards inferring earthquake patterns from geodetic observations of interseismic coupling[END_REF][START_REF] Moreno | Locking of the Chile subduction zone controlled by fluid pressure before the 2010 earthquake[END_REF], inherited stress states [START_REF] Konca | Partial rupture of a locked patch of the Sumatra megathrust during the 2007 earthquake sequence[END_REF][START_REF] Kaneko | Towards inferring earthquake patterns from geodetic observations of interseismic coupling[END_REF] or seafloor roughness [START_REF] Kodaira | Subducted seamount imaged in the rupture zone of the 1946 Nankaido earthquake[END_REF][START_REF] Wang | Invited review paper: Fault creep caused by subduction of rough seafloor relief[END_REF][START_REF] Bletery | Mega-earthquakes rupture flat megathrusts[END_REF]. Megathrusts, commonly regarded as interfingered seismic and aseismic patches with contrasting frictional properties, are conventionally sub-divided into four domains with depth (A-D, [START_REF] Lay | Depth-varying rupture properties of subduction zone megathrust faults[END_REF]. The highly coupled domain-B hosts mega-earthquakes, while the partly coupled domain-C is characterized by moderate slip earthquakes generating significant high frequency seismic radiation [START_REF] Lay | Depth-varying rupture properties of subduction zone megathrust faults[END_REF].

The distribution of these patches might persist over time scales up to million years and control forearc morphology [START_REF] Song | Large trench-parallel gravity variations predict seismogenic behavior in subduction zones[END_REF]. A spatial correlation between the coastline and the down-dip limit of the strongly locked zone has indeed been evidenced along the northern chilean margin [START_REF] Béjar-Pizarro | Andean structural control on interseismic coupling in the North Chile subduction zone[END_REF][START_REF] Metois | Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38-18 S)[END_REF]. Coastal uplift rates inferred from marine terraces, between 0.1 and 0.3 mm/a over the last 400 ka [START_REF] Saillard | From the seismic cycle to long-term deformation: linking seismic coupling and Quaternary coastal geomorphology along the Andean megathrust[END_REF], have been linked to the change of mechanical coupling between domains B and C [START_REF] Saillard | From the seismic cycle to long-term deformation: linking seismic coupling and Quaternary coastal geomorphology along the Andean megathrust[END_REF] or to underplating [START_REF] Delouis | Recent crustal deformation in the Antofagasta region (northern Chile) and the subduction process[END_REF][START_REF] Adam | Crustal dynamics and active fault mechanics during subduction erosion. Application of frictional wedge analysis on to the North Chilean Forearc[END_REF][START_REF] Clift | Slow rates of subduction erosion and coastal underplating along the Andean margin of Chile and Peru[END_REF]. Numerical investigations recently suggested that underplating, through transient stripping of the slab-top near the base of the forearc crust, may generate periodic Myr-long uplift sequences and a trackable 100 m-high topographic signal [START_REF] Menant | Transient stripping of subducting slabs controls periodic forearc uplift[END_REF].

Bathymetric features are considered as potential seismic barriers, based on the distribution of incoming plate roughness [START_REF] Wang | Invited review paper: Fault creep caused by subduction of rough seafloor relief[END_REF]Lallemand et al., 2018) or gravimetric anomalies [START_REF] Bassett | Gravity anomalies, crustal structure, and seismicity at subduction zones: 2. Interrelationships between fore-arc structure and seismogenic behavior[END_REF]. Subduction of bathymetric features induces large fracture networks, promoting removal of material from the bottom of the overriding plate (basal erosion, [START_REF] Ranero | Subduction erosion along the Middle America convergent margin[END_REF] and potentially characteristic topographic features [START_REF] Dominguez | Upper plate deformation associated with seamount subduction[END_REF][START_REF] Collot | Origin of a crustal splay fault and its relation to the seismogenic zone and underplating at the erosional north Ecuador-south Colombia oceanic margin[END_REF][START_REF] Ruh | Crustal deformation dynamics and stress evolution during seamount subduction: High-resolution 3-D numerical modeling[END_REF]. Hetereogeneities in the structure of the damaged plate interface and stress field, as a result, would favor aseismic slip and impede the propagation of large ruptures [START_REF] Wang | Invited review paper: Fault creep caused by subduction of rough seafloor relief[END_REF]. Very few seismic surveys, however, have imaged seamounts surrounding a large rupture [START_REF] Kodaira | Subducted seamount imaged in the rupture zone of the 1946 Nankaido earthquake[END_REF][START_REF] Geersen | Subducting seamounts control interplate coupling and seismic rupture in the 2014 Iquique earthquake area[END_REF].

We herein explore an alternative explanation: rather than variations of frictional properties or of seafloor roughness, earthquake segmentation might relate to the distribution of deformation along the plate interface [START_REF] Wang | Invited review paper: Fault creep caused by subduction of rough seafloor relief[END_REF], as observed along strike-slip faults (Wesnousky, 2006). Both underplating and basal erosion require the redistribution of deformation ultimately leading to the migration of the plate boundary [START_REF] Vannucchi | Toward a dynamic concept of the subduction channel at erosive convergent margins with implications for interplate material transfer[END_REF][START_REF] Agard | The subduction plate interface: rock record and mechanical coupling (from long to short timescales)[END_REF]. We here show that the location of such distributed deformation along the plate interface can be captured from a simple mechanical analysis of the topography, since underplating and basal erosion both impact forearc morphology. This method is applied over a 2000km long strech along the Chilean margin, known to transition from erosive to accretionary from north to south, as incoming sediment thickness increases [START_REF] Von Huene | Subduction erosion and basal friction along the sediment-starved convergent margin off Antofagasta, Chile[END_REF][START_REF] Clift | Controls on tectonic accretion versus erosion in subduction zones: Implications for the origin and recycling of the continental crust[END_REF][START_REF] Clift | Slow rates of subduction erosion and coastal underplating along the Andean margin of Chile and Peru[END_REF] (Fig. 1a). This margin is one of the best instrumented subduction zones, along which five mega-earthquakes, including a domain-C event, have occurred during the last 25 years (Fig. 1c). We first map out the location of distributed deformation zones associated with basal ero-sion and underplating and compare them with the megathrust segmentation proposed so far. We then discuss how distributed deformation along the plate interface impacts and therefore can be used to constrain earthquake extent and mechanics.

Methodology

We applied the critical taper theory (CTT, Davis et al., 1983, Fig. 2). The theory predicts the topographic slope of a wedge on the verge of failure as well as dips of its internal faults, from the frictional and pore fluid pressure properties of the wedge and megathrust [START_REF] Dahlen | Noncohesive critical Coulomb wedges: An exact solution[END_REF]. We used the solution of [START_REF] Dahlen | Noncohesive critical Coulomb wedges: An exact solution[END_REF] for a non cohesive wedge, describing the critical taper, formed by the topographic slope α and the slab dip β, as a function of the angle Ψ B formed by the maximum principal stress σ 1 and the base of the wedge and the angle Ψ 0 formed by σ 1 and the top of the wedge. The solution for the lower compressional branch is:

(α + β) c = Ψ B -Ψ 0 (1)
with

Ψ B = 1 2 arcsin( sinφ b sinφ b ) - 1 2 φ b , (2)
Ψ 0 = 1 2 arcsin( sinα sinφ int ) - 1 2 α . (3)
The angles φ int and φ b are the internal and basal coefficients of friction defined as µ int = tanφ int and µ b = tanφ b , and

φ b = arctan 1 -λ b 1 -λ tanφ b , (4)
α = arctan 1 -ρ w /ρ 1 -λ tanα . (5)
The internal and basal Hubbert-Rubbey fluid pressure ratios λ and λ b are defined in [START_REF] Davis | Mechanics of fold-and-thrust belts and accretionary wedges[END_REF] as:

λ = P -ρ w gD |σ z | -ρ w gD , (6)
λ b = P b -ρ w gD |σ z | -ρ w gD , (7)
where ρ and ρ w are the wedge material and water densities and D is the water depth. The solution is exact if λ = λ b and the approximation is valid for small tapers as used in this study [START_REF] Wang | A note on pore fluid pressure ratios in the Coulomb wedge theory[END_REF].

The relationship between the topographic slope α and the slab dip β forms an envelope separating different mechanical states (Fig. 2). The wedge is at critical state if the taper formed by α and β follows a critical envelope. In that case, activation of the megathrust requires internal faulting. In contrast, if the taper exceeds this critical limit, the wedge enters in a stable domain (in the sense of the CTT), where the only possible active fault is the megathrust.

Dynamic effective frictions of seismogenic zones have been shown to reach extremely low values (µ ∼ 0.01 -0.03, [START_REF] Fulton | Low coseismic friction on the Tohoku-Oki fault determined from temperature measurements[END_REF][START_REF] Gao | Strength of stick-slip and creeping subduction megathrusts from heat flow observations[END_REF]. As a consequence, a forearc above a seismogenic megathrust, with standard α and β, will systematically fall in the CTT stable domain (Fig. 2). In contrast, aseismic megathrusts are characterized by larger effective friction (µ ∼ 0.1 -0.15, [START_REF] Cubas | Megathrust friction determined from mechanical analysis of the forearc in the Maule earthquake area[END_REF][START_REF] Gao | Strength of stick-slip and creeping subduction megathrusts from heat flow observations[END_REF], allowing the wedge to reach critical conditions [START_REF] Cubas | Megathrust friction determined from mechanical analysis of the forearc in the Maule earthquake area[END_REF] (Fig. 2). Furthermore, at critical state, if the basal effective friction of the wedge approaches the internal effective friction, the dip of internal faults within the wedge will decrease and become parallel to the plate interface. A small difference of friction will thus lead to deformation either by basal accretion (in the lower plate) or by basal erosion (in the upper plate), instead of standard accretion characterized by thrust faults reaching the surface [START_REF] Dahlen | Noncohesive critical Coulomb wedges: An exact solution[END_REF] (Fig. 2). If the basal effective friction reaches the internal one, a highly fractured forearc is even expected.

We here used inversion to retrieve the effective basal friction, the internal friction and internal pore pressure ratio from the α -β co-variation. To do so, we followed the procedure developed by [START_REF] Cubas | Megathrust friction determined from mechanical analysis of the forearc in the Maule earthquake area[END_REF] and detailed in Supplementary Materials. We used ETOPO 1 [START_REF] Amante | ETOPO1 arc-minute global relief model: procedures, data sources and analysis[END_REF] to derive the bathymetric and topographic slope, and slab 2.0 [START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF] for the slab dip. We built swath profiles perpendicular to the trench every 0.1 degree along trench to get average topographic slopes and slab dips with their standard deviations necessary to the inversion procedure (Suppl. Mat. Fig. S1).

We explored from 10 to 45 o for φ int , from 1 to 44.9 o for φ ef f b , and from 0.35 to 0. 975 for λ. Along the α -β profiles, we selected segments with a slope comprised between the two extreme theoretical envelops of the parameter set (Fig. 3, Suppl.

Mat. Fig. S2). For each of these segments, we retrieved probability density functions for the friction of the wedge, the pore pressure ratio of the wedge and the effective friction of the megathrust (Suppl. Mat. Fig. S3) and best misfit sets (Suppl. Mat. Fig. S4). We only kept segments with misfits lower than 0.1, and with values consistent with standard frictions: from 25 to 43 o for φ int (consistent with lab experiment values [START_REF] Byerlee | Friction of rocks[END_REF]), from 1 to 42.9 o for φ ef f b (same range as the internal friction but considering additional pore fluid pressure effect), and from 0.35 to 0. 975 for λ. Figure 3a-c shows the fit between the α-β profiles and the envelops corresponding to the best misfit of the identified critical segments. These segments (Fig. 3b-d) have then been reported on Fig. 1b.

The CTT inversion only provides the internal pore pressure, which can be considered as a lower bound for the megathrust pore pressure. Moreover, close effective internal and basal frictions implies close internal and megathrust pore pressure.

Results

Locating plate interface deformation

Since we seek to capture distribution of deformation along the plate interface, we compared the internal and megathrust effective frictions and mapped the difference (Fig. 1b). A one degree difference implies a dip difference < 10 o between the megathrust and the forward-verging thrusts. This representation allows inferring areas prone to either standard accretion, or to basal erosion/accretion, the latter being herein defined as distributed deformation along the plate interface.

The transition from erosion to accretion is particularly well captured by the method despite its relative simplicity (Fig. 1a-b). North of 26 o S (i.e., at the transition between erosive and accretionary subduction), the difference of effective frictions is systematically lower than 2 o , and mostly lower than 1 o , which means that deformation is essentially located along the plate interface. The difference increases notably south of the Juan Fernandez Ridge, and progresses inland towards the south, consistent with active faults mapped near the Arauco peninsula (Fig. 1d). The method also captures basal erosion linked to the few subducted seamounts imaged around the Iquique rupture [START_REF] Geersen | Subducting seamounts control interplate coupling and seismic rupture in the 2014 Iquique earthquake area[END_REF] (Fig. 1c).

To interpret the origin of plate interface deformation, we set the limit between accretion and interface deformation for a difference of friction angle of one degree and then searched for the depth distribution of each process (Fig. 4a-b). Results for different thresholds are also presented.

Accretion mostly occurs at depths shallower than 20 km whereas distributed deformation along the plate interface is found at every depth with three favorable peaks: a first one at 15 km depth, a second at 35±10 km depth and a third at 60±5 km depth.

Basal erosion due to bathymetric features is known to be efficient at shallow depth [START_REF] Von Huene | Generic model of subduction erosion[END_REF], but might also occur slightly deeper as testified by subsidence of forearc basins [START_REF] Clift | Slow rates of subduction erosion and coastal underplating along the Andean margin of Chile and Peru[END_REF]. Underplating has been observed at shallow depth [START_REF] Kimura | Seismic evidence for active underplating below the megathrust earthquake zone in Japan[END_REF][START_REF] Tréhu | Geologic controls on up-dip and along-strike propagation of slip during subduction zone earthquakes from a high-resolution seismic reflection survey across the northern limit of slip during the 2010 Mw 8.8 Maule earthquake, offshore Chile[END_REF][START_REF] Bangs | Basal accretion along the south central Chilean margin and its relationship to great earthquakes[END_REF], even for seamount moats [START_REF] Clarke | Seamount chain-subduction zone interactions: Implications for accretionary and erosive subduction zone behavior[END_REF], but mostly occurs below 30-40 km, i.e. the long-term coupling-decoupling transition [START_REF] Agard | The subduction plate interface: rock record and mechanical coupling (from long to short timescales)[END_REF]. The compilation of P-T conditions for recovered underplated material shows a similar distribution as CTT, with peaks at depths of 30±5 km and 50±5 km and a gap in between (Fig. 4b). Numerical modelling also suggests that erosion predominates over underplating at depths shallower than 20 km [START_REF] Menant | Transient stripping of subducting slabs controls periodic forearc uplift[END_REF]. As a consequence, we assume that the interface deformation documented here mostly relates to basal erosion for depths < 20 km, whereas underplating dominates at greater depths, in particular below the coast. This depth divide between basal erosion and underplating is strengthened by estimates of pore fluid pressure ratio (Fig. 4d-f).

If low pore fluid pressure is mainly found in accretion conditions, higher fluid pressure dominates where basal erosion occurs [START_REF] Von Huene | Subduction erosion and basal friction along the sediment-starved convergent margin off Antofagasta, Chile[END_REF][START_REF] Vannucchi | Toward a dynamic concept of the subduction channel at erosive convergent margins with implications for interplate material transfer[END_REF]. On the contrary, recovered underplated material from greater depth present very consistent P-T conditions of lithostatic pressure, which rules out any significant fluid overpressure [START_REF] Agard | The subduction plate interface: rock record and mechanical coupling (from long to short timescales)[END_REF]. The first peak for basal erosion at low pore fluid pressures in Fig. 4e could be related to the 20 km depth assumption and reveal the minor occurrence of underplating at shallower depth. The trimodal distribution deduced from CCT (one peak for basal erosion and two peaks for underplating) and differences in pore pressure remains for a 0.5 o or a 2 o difference of friction.

Most of the critical segments attesting to distributed interface deformation have lengths smaller than 7 km, and do not exceed 20 km (Fig. 4c). This is in good agreement with the length of graben and horst and seamounts observed on the subducting plate [START_REF] Von Huene | Subduction erosion and basal friction along the sediment-starved convergent margin off Antofagasta, Chile[END_REF][START_REF] Geersen | Subducting seamounts control interplate coupling and seismic rupture in the 2014 Iquique earthquake area[END_REF]. This is also consistent with seismic observations [START_REF] Kimura | Seismic evidence for active underplating below the megathrust earthquake zone in Japan[END_REF][START_REF] Tréhu | Geologic controls on up-dip and along-strike propagation of slip during subduction zone earthquakes from a high-resolution seismic reflection survey across the northern limit of slip during the 2010 Mw 8.8 Maule earthquake, offshore Chile[END_REF][START_REF] Bangs | Basal accretion along the south central Chilean margin and its relationship to great earthquakes[END_REF] and the long-term rock record showing generally ≤ 300-500 m thick exhumed slices from the top of the slab [START_REF] Agard | The subduction plate interface: rock record and mechanical coupling (from long to short timescales)[END_REF] extending 5-10 km downdip at most [START_REF] Plunder | How continuous and precise is the record of P-T paths? Insights from combined thermobarometry and thermodynamic modelling into subduction dynamics (Schistes Lustrés, W. Alps)[END_REF][START_REF] Plunder | Geodynamics of the Tavşanlı zone, western Turkey: Insights into subduction/obduction processes[END_REF]. Such tabular shape ratios, with seamounts probably near the higher bound [START_REF] Bonnet | No large earthquakes in fully exposed subducted seamount[END_REF], further support interface deformation parallel to the megathrust [START_REF] Fagereng | Melange rheology and seismic style[END_REF][START_REF] Rowe | The thickness of subduction plate boundary faults from the seafloor into the seismogenic zone[END_REF].

Relating plate interface deformation with long-term coastal uplift

In order to approach internal effective friction, the basal friction of the wedge needs to be relatively high (Fig. 2a). Deformation captured at depth by CTT therefore corresponds to underplating in the making, when deformation is still distributed: once tectonic slicing is achieved, basal friction will reach low values characteristic of mature faults, bringing the taper in the stable domain and impeding detection by CTT.

We here show that active, incipient underplating takes place along most of the Chilean margin (Fig. 1, 5a), whether or not strain accumulation ultimately turns these critical zones into tectonic slices [START_REF] Agard | The subduction plate interface: rock record and mechanical coupling (from long to short timescales)[END_REF][START_REF] Bonnet | No large earthquakes in fully exposed subducted seamount[END_REF][START_REF] Kimura | Transition of accretionary wedge structures around the up-dip limit of the seismogenic subduction zone[END_REF]. Along the Tocopilla and Maule segments, incipient underplating takes place below the slab -continental Moho intercept (orange curve on Fig. 1c), whereas it is located above and near the intercept between 25-33 o S. Along the northern erosive margin (Fig. 5a-b), the thin oceanic sediment cover together with eroded material and possibly pieces stripped from the subducting plate are being underplated. Southward, larger sediment input favors shallow accretion and underplating of subducted sediments (Fig. 5a-b), as locally observed by seismic surveys [START_REF] Tréhu | Geologic controls on up-dip and along-strike propagation of slip during subduction zone earthquakes from a high-resolution seismic reflection survey across the northern limit of slip during the 2010 Mw 8.8 Maule earthquake, offshore Chile[END_REF].

To evaluate the possible link between coastal uplift and regions of distributed deformation along the plate interface, we compare our spatial distribution of underplating with documented marine terraces and uplift rates [START_REF] Saillard | From the seismic cycle to long-term deformation: linking seismic coupling and Quaternary coastal geomorphology along the Andean megathrust[END_REF] (Fig. 5a-c). While the irregular distribution of marine terraces may relate to fragmentary preservation rather than discontinuous underplating, most of them coincide with segments of underplating (or with standard accretion near the Arauco peninsula).

The typical length of the segments showing distributed deformation along the interface is modest (∼10-20 km). This advocates for propagation of deformation only into the topmost part of the subducting plate, e.g. into the stack of oceanic/eroded material or along a hydrated basalt layer within the crust. Although the comparison with uplift rates inferred from morphometric analysis and modelling of landscape evolution [START_REF] Melnick | Rise of the central Andean coast by earthquakes straddling the Moho[END_REF] shows no particular pattern, underplating of relatively thin slices would be consistent with both the low uplift rates and the rock record [START_REF] Agard | The subduction plate interface: rock record and mechanical coupling (from long to short timescales)[END_REF].

Shallow normal faults identified along northern Chile have been related to potential underplating [START_REF] Adam | Crustal dynamics and active fault mechanics during subduction erosion. Application of frictional wedge analysis on to the North Chilean Forearc[END_REF][START_REF] Clift | Slow rates of subduction erosion and coastal underplating along the Andean margin of Chile and Peru[END_REF]. We here show that they are located above underplating areas, and that normal faults transition southward to thrusts where accretion prevails (Fig. 5a-c).

Large megathrust earthquakes are surrounded by plate interface deformation

We now compare zones of distributed deformation with areas struck by recent large earthquakes (Fig. 1c). The Mw 8.1 2014 Iquique earthquake is surrounded by patches of distributed deformation associated with basal erosion according to their depth and to seamounts detected in the area [START_REF] Geersen | Subducting seamounts control interplate coupling and seismic rupture in the 2014 Iquique earthquake area[END_REF]. The aftershock, located 50 km south of the main event, also falls in between patches of distributed deformation. The Mw 7.8 2007 Tocopilla earthquake is limited up-dip as well as to the south by interface deformation attributed to underplating. The Mw 8.1 1995 Antofagasta earthquake seems limited up-and down-dip by distributed deformation along the plate interface. The high slip areas of the Mw 8.2 2015 Illapel and the Mw 8.7 2010 Maule earthquakes are both limited down-dip by plate interface deformation. Their frontal extent partially propagates in areas at critical state but corresponding to standard accretion. A small portion of the earthquake ruptures do overlap with the plate interface deformation, but high slip areas are always clearly limited by plate interface deformation. Moreover, Fig. 6 shows the diversity of published co-seismic models, emphasizing the uncertainties related to the proposed slip areas, that need to be taken into account for the comparison.

To summarize, along the accretionary part of the margin, mega-earthquakes are delimited down-dip by underplating whereas along the erosive domain, the extent is controlled by both basal erosion and underplating.

We now investigate the relationships between the critical segments characterized by relatively high effective basal friction, hence expected to behave aseismically, both with the along-strike segmentation inferred for large historical earthquakes (Mw ≥ 7.5, Saillard et al., 2017, Fig. 5a-d) and with the estimates of coupling [START_REF] Metois | Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38-18 S)[END_REF]Klein et al., 2018, Fig. 5a-e).

Segments with limited plate interface deformation (<15 km) coincide with a significant number of rupture terminations (blue and grey overlays, respectively: Fig. 5a-d), particularly in the southern region. Between 29 and 24 o S, rupture terminations coincide with segments where plate interface deformation is reduced. The underplating segments also show a long wavelength distribution (4 to 8 o degrees latitude) corresponding to the three main seismic regions visible on Fig. 5a,d (respectively located north of Mejillones, between Mejillones and Punta Choros, and to the south of Punta Choros). No correlation exists with sediment input (Fig. 5b). The proportion of a swath profile at critical state is then compared to the average estimate of coupling along this profile (Fig. 5e). A positive correlation is evidenced at small scale (1.5 and 2 o degrees latitude, Suppl. Mat. fig. S7) between the segments with minimal underplating (blue overlays; Fig. 5a) and the regions of minimal coupling (red overlays; Fig. 5e). On a larger scale (4 and 8 o degrees latitude, Suppl. Mat. fig. S7), an anti-correlation is observed with coupling: zones of high coupling correspond to domains where underplating is somewhat more limited (e.g., between 31 and 38 o S).

As shown by Fig. 1d, areas marked by intense plate deformation are mostly located down-dip of the highly locked patches, in particular in the south. Continuous and voluminous sediment influx provides spatially and temporally stable mechanical conditions [START_REF] Olsen | Thick, strong sediment subduction along south-central Chile and its role in great earthquakes[END_REF]. This leads to stationary seismic asperities followed, once the Moho is crossed, by the regular underplating of very thin and relatively short slices (typically 300-500 m thick, 5-10 km long). This spatial and temporal stability promotes topographic build-up, explaining the spatial coincidence between strongly coupled patches and efficient underplating. In contrast, segments with limited underplating might reveal spatially and temporally heterogeneous conditions along the megathrust [START_REF] Olsen | Thick, strong sediment subduction along south-central Chile and its role in great earthquakes[END_REF], i.e. heterogeneous deformation impeding earthquake propagation and significant topographic build-up.

Noteworthily, the hypocenters of the last major events lie at the boundary with a patch of distributed deformation. Finally, we show that the high frequency radiations generated by Maule and Illapel earthquakes [START_REF] Wang | Frequency-dependent energy radiation and fault coupling for the 2010 Mw8. 8 Maule, Chile[END_REF][START_REF] Meng | Dual megathrust slip behaviors of the 2014 Iquique earthquake sequence[END_REF][START_REF] Meng | Double pincer movement: Encircling rupture splitting during the 2015 Mw 8.3 Illapel earthquake[END_REF] coincide with areas of extensive interplate deformation (Fig. 1c), along which spatially restricted fault planes might have been dynamically triggered.

Link between earthquake ruptures and plate interface deformation

We herein show two major results: (1) recent earthquakes are bounded by extensive plate interface deformation characterized, along the southern accretionary part of the margin by underplating at their down-dip edge and, along the northern erosive part of the margin by both basal erosion and underplating (Fig. 1); (2) along strike, segments characterized by minimal interplate deformation largely coincide with lateral earthquake terminations (circles; Fig. 5d). These observations demonstrate a close link between distributed interplate deformation and earthquake segmentation.

We therefore propose that plastic deformation and stress build-up associated with interplate deformation along distributed fault planes of limited extent (Fig. 7 -step 1) eventually leads to mega-earthquake nucleation (Fig. 7 -step 2). This is suggested by the location of hypocenters at boundaries with distributed deformation (Fig. 1). This interseismic deformation might be accompanied by some micro-seismicity, as observed prior to the Iquique rupture [START_REF] Schurr | Forming a Mogi doughnut in the years prior to and immediately before the 2014 M8. 1 Iquique, northern Chile, earthquake[END_REF] and in the upper plate in northern Chile [START_REF] Sippl | Seismicity structure of the northern Chile forearc from> 100,000 double-difference relocated hypocenters[END_REF]. If high pore fluid pressure prevails, this deformation could also be accompanied by transient seismic/slip events [START_REF] Liu | Spontaneous and triggered aseismic deformation transients in a subduction fault model[END_REF][START_REF] Kimura | Seismic evidence for active underplating below the megathrust earthquake zone in Japan[END_REF][START_REF] Collot | Subducted oceanic relief locks the shallow megathrust in central Ecuador[END_REF], which concords with observations for the Iquique earthquake (Ruiz et al., 2014a;[START_REF] Meng | Dual megathrust slip behaviors of the 2014 Iquique earthquake sequence[END_REF][START_REF] Socquet | An 8 month slow slip event triggers progressive nucleation of the 2014 Chile megathrust[END_REF] and others elsewhere (Tohoku, Kato et al. (2012); Guerrero, [START_REF] Radiguet | Triggering of the 2014 M w 7.3 Papanoa earthquake by a slow slip event in Guerrero, Mexico[END_REF]).

Once nucleated, large earthquakes propagate along well localized and smoothed rate-weakening fault planes [START_REF] Bletery | Mega-earthquakes rupture flat megathrusts[END_REF] limited by elongate zones of underplating which, in addition, inhibit further rupture propagation and slip (Fig. 1, 7 -step 3).

These earthquakes abut on regions of heterogeneously distributed deformation and stress concentrations (Fig. 7 -step 4), impeding the development of well localized slip zones [START_REF] Wang | Invited review paper: Fault creep caused by subduction of rough seafloor relief[END_REF], and again significant slip (Wesnousky, 2006). Some of these fractures might however slip seismically during large events and produce the observed high frequency radiations [START_REF] Meng | Double pincer movement: Encircling rupture splitting during the 2015 Mw 8.3 Illapel earthquake[END_REF][START_REF] Meng | Dual megathrust slip behaviors of the 2014 Iquique earthquake sequence[END_REF].

Down-dip segmentation and plate interface deformation

Within this framework, domain-C [START_REF] Lay | Depth-varying rupture properties of subduction zone megathrust faults[END_REF] would correspond to the region of incipient underplating, and domain-B to the smoothed rate-weakening fault zone. This is consistent with partial coupling of domain-C, where deformation is dominated by creep mechanisms, and elastic strain is accumulated along undeformed segments of the megathrust, producing moderate slip earthquakes with a higher recurrence time and significant coherent short-period seismic radiations. Domain-C might more specifically be delimited by the two underplating peaks (Fig. 4), the upper one corresponding to the slicing of oceanic sediments and/or eroded continental material, the second to the propagation of deformation into the hydrated basalt layer once the continental Moho is crossed [START_REF] Agard | The subduction plate interface: rock record and mechanical coupling (from long to short timescales)[END_REF], a process ultimately leading to the underplating of very thin and relatively short slices (typically 300-500 m thick, 5-10 km long, [START_REF] Agard | The subduction plate interface: rock record and mechanical coupling (from long to short timescales)[END_REF].

Data suggest that domain-C, by locally concentrating interplate permanent deformation due to specific geometry (i.e., morphological asperities) and/or mechanical behaviour (i.e., frictional properties, fluid content, porosity), shields the rate-weakening domain-B from long-lived deformation (Fig. 7). Contrary to the prevalent paradigm, domain-B should not be regarded as an 'asperity', i.e. a strongly coupled yet mechanically weak zone with high recurrence time earthquakes, but instead as a relatively weak zone devoid of permanent deformation and only storing elastic energy. In this new paradigm, loading along domain-C tightly controls slip deficit as well as rupture nucleation and extent along domain-B (Fig. 7), explaining the correlation observed between historical earthquakes, earthquake terminations and distributed plate interface deformation (Fig. 5 a,d). On the short-term scale, elastic deformation generates subsidence at the coast during domain-B earthquakes and uplift for domain-C events [START_REF] Melnick | Rise of the central Andean coast by earthquakes straddling the Moho[END_REF]. On the long-term, the thin underplated slices ultimately build a Myr topographic signal of slow uplift, particularly difficult to capture with geodetic measurements due to the interplay of elastic and plastic deformation [START_REF] Menant | Transient stripping of subducting slabs controls periodic forearc uplift[END_REF]. The Coastal Cordillera present along much of the Chilean margin could result from this long-term underplating process. The absence of upper plate deformation above the seismogenic zone in conjunction with topographic build-up above underplating areas account for the correlation between seismic behavior and gravity anomalies [START_REF] Song | Large trench-parallel gravity variations predict seismogenic behavior in subduction zones[END_REF].

Conclusions

We herein show, through a simple CTT analysis of topography, that, along the Chilean subduction zone, earthquakes are bounded by long-lived, extensive plate interface deformation occurring preferentially at 35 ± 10 and 60 ± 5 km depth. Along the Central Chilean margin, earthquakes are limited along their down-dip edge by underplating while, along Northern Chile, they are surrounded by distributed deformation related to both basal erosion and underplating. Plate interface deformation might control earthquake nucleation, extent and arrest, and its distribution provide a relatively reliable image of the seismic potential at depth. This topography analysis may also remedy the lack of preserved marine terraces or weakely constrained coupling models. This method could therefore be used more systematically to identify regions with precursory transient signals and to enhance predictions of future earthquakes extent.

Appendix See supplementary material for detailed method and supporting figures. [START_REF] Santibáñez | Crustal faults in the Chilean Andes: geological constraints and seismic potential[END_REF][START_REF] Tréhu | Geologic controls on up-dip and along-strike propagation of slip during subduction zone earthquakes from a high-resolution seismic reflection survey across the northern limit of slip during the 2010 Mw 8.8 Maule earthquake, offshore Chile[END_REF]. (b) Difference between internal and basal effective frictions of identified critical segments. Large differences lead to standard accretion with thrust faults reaching the surface, while small differences lead to the development of faults parallel to the plate interface allowing for either basal erosion or underplating. (c) Interpreted patches with distributed deformation along the plate interface related to either basal erosion or underplating compared to last major events: Antofagasta 1995 [START_REF] Chlieh | Crustal deformation and fault slip during the seismic cycle in the North Chile subduction zone, from GPS and InSAR observations[END_REF]), Tocopilla 2007[START_REF] Béjar-Pizarro | Andean structural control on interseismic coupling in the North Chile subduction zone[END_REF]), Maule 2010[START_REF] Fulton | Low coseismic friction on the Tohoku-Oki fault determined from temperature measurements[END_REF]; Iquique 2014 and its aftershock (Ruiz et al., 2014b), Illapel 2015 [START_REF] Tilmann | The 2015 Illapel earthquake, central Chile: A type case for a characteristic earthquake?[END_REF]. Orange stars: gCMT solutions, blue stars: CSN catalog. High frequencies radiations (yellow, [START_REF] Wang | Frequency-dependent energy radiation and fault coupling for the 2010 Mw8. 8 Maule, Chile[END_REF][START_REF] Meng | Dual megathrust slip behaviors of the 2014 Iquique earthquake sequence[END_REF][START_REF] Meng | Double pincer movement: Encircling rupture splitting during the 2015 Mw 8.3 Illapel earthquake[END_REF], identified seamounts (green, [START_REF] Geersen | Subducting seamounts control interplate coupling and seismic rupture in the 2014 Iquique earthquake area[END_REF][START_REF] Maksymowicz | Density-depth model of the continental wedge at the maximum slip segment of the Maule Mw8. 8 megathrust earthquake[END_REF], slab depths (dashed black, [START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF] and slab-Moho intercept (orange, [START_REF] Tassara | Anatomy of the Andean subduction zone: three-dimensional density model upgraded and compared against global-scale models[END_REF]. (d) Interpreted patches with distributed plate interface compared to coupling model [START_REF] Metois | Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38-18 S)[END_REF] with known active faults [START_REF] Santibáñez | Crustal faults in the Chilean Andes: geological constraints and seismic potential[END_REF]. (green, Santibáñez et al., 2018;[START_REF] Tréhu | Geologic controls on up-dip and along-strike propagation of slip during subduction zone earthquakes from a high-resolution seismic reflection survey across the northern limit of slip during the 2010 Mw 8.8 Maule earthquake, offshore Chile[END_REF]. (c) Quaternary uplift rates of marine terraces (black bars, largest average rate since terrace abandonment, [START_REF] Saillard | From the seismic cycle to long-term deformation: linking seismic coupling and Quaternary coastal geomorphology along the Andean megathrust[END_REF], coastal uplift rate obtained from a landscape evolution model (black curve, [START_REF] Melnick | Rise of the central Andean coast by earthquakes straddling the Moho[END_REF] and extent of known active faults (NF: normal fault, TF: thrust fault, [START_REF] Santibáñez | Crustal faults in the Chilean Andes: geological constraints and seismic potential[END_REF]. (d) Alongstrike extent of historical megathrust ruptures [START_REF] Saillard | From the seismic cycle to long-term deformation: linking seismic coupling and Quaternary coastal geomorphology along the Andean megathrust[END_REF]. (e) Comparison between the cumulative length of segments at critical state (thin black), with underplating (thick blue), with both underplating and basal erosion (light blue) and the average coupling value along profiles from [START_REF] Metois | Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38-18 S)[END_REF] corrected by [START_REF] Klein | Bridging the gap between North and Central Chile: insight from new GPS data on coupling complexities and the Andean sliver motion[END_REF] between 24.5 o S and 28 o S (plain red curve). Red bands show segments of decreasing coupling. Correlation factor between coupling and underplating: r = -0.36, between coupling and underplating + basal erosion: r = -0.39. 20 Figure 6. Areas prone to distributed deformation along the plate interface compared to various published co-seismic slip models [START_REF] Duputel | The Iquique earthquake sequence of April 2014: Bayesian modeling accounting for prediction uncertainty[END_REF][START_REF] Schurr | The 2007 M7. 7 Tocopilla northern Chile earthquake sequence: Implications for along-strike and downdip rupture segmentation and megathrust frictional behavior[END_REF][START_REF] Pritchard | An aseismic slip pulse in northern Chile and along-strike variations in seismogenic behavior[END_REF][START_REF] Tilmann | The 2015 Illapel earthquake, central Chile: A type case for a characteristic earthquake?[END_REF][START_REF] Vigny | The 2010 Mw 8.8 Maule megathrust earthquake of central Chile, monitored by GPS[END_REF]Ruiz et al., 2014b;[START_REF] Chlieh | Crustal deformation and fault slip during the seismic cycle in the North Chile subduction zone, from GPS and InSAR observations[END_REF][START_REF] Moreno | Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone[END_REF][START_REF] Hayes | Continuing megathrust earthquake potential in Chile after the 2014 Iquique earthquake[END_REF][START_REF] Béjar-Pizarro | Andean structural control on interseismic coupling in the North Chile subduction zone[END_REF][START_REF] Fulton | Low coseismic friction on the Tohoku-Oki fault determined from temperature measurements[END_REF]. Orange stars: gCMT solution, blue stars: CSN solution.

Figure 1 .

 1 Figure 1. (a) Sediment thickness at trench (Santibáñez et al., 2018; Tréhu et al., 2019). (b) Difference between internal and basal effective

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Critical taper theory with different mechanical states. The grey rectangle represents the standard range of topographic slope α and slab dip β of subduction zones. With extremely low effective friction along the megathrust (red curve), the wedge systematically falls in the stable domain. With larger effective friction, and large difference of effective friction between the megathrust and the wedge (green curve),the wedge can reach the accretionary critical state characterized by internal faulting. With a reduced difference (blue curves), a wedge can still reach the critical state but the dip of internal faults within the wedge will decrease and become parallel to the plate interface. High pore fluid pressure facilitates interface deformation for low slab dips (dashed blue curve), whereas lower pore fluid pressure induces interface deformation at larger slab dip (plain blue curve).

Figure 5 .

 5 Figure 5. (a) Along-dip cumulative length of segments at accretionary (green), underplating (dark blue) or basal erosion (orange) critical state (length smoothed by a rectangular window over 8 profiles). The light blue curve shows the cumulative length of segments with interface deformation (basal erosion + underplating). Regions with limited underplating (and interface deformation) length are represented by blue bands. (b) Thickness of sediment income(green, Santibáñez et al., 2018;[START_REF] Tréhu | Geologic controls on up-dip and along-strike propagation of slip during subduction zone earthquakes from a high-resolution seismic reflection survey across the northern limit of slip during the 2010 Mw 8.8 Maule earthquake, offshore Chile[END_REF]. (c) Quaternary uplift rates of marine terraces

TABLE No

 No

	internal def.

Wesnousky, S. G.: Predicting the endpoints of earthquake ruptures,Nature, 444, 358-360, 2006.

Acknowledgements. We thank the two anonymous reviewers for their helpful comments. This project has been funded by French National Research Agency (ANR) grant SEAFRONTTERA ANR-13-PDOC-0013-01.

275

Code and data availability. Data available at to be completed and code available on demand.

Author contributions. NC designed the study; NC and RT ran the simulations; NC, RT, PA analyzed the results; NC and PA wrote the paper.

Competing interests. Authors declare no competing interest.


L Daukiya

M N Nair

S Hajjar-Garreau

F Vonau

D Aubel

J L Bubendorff

M Cranney

E Denys

A Florentin

G Reiter

L Simon

Highly n-Doped Graphene Generated Through Intercalated Terbium Atoms

Keywords: numbers: 68.65.-k, 81.07.-b, 81.05.ue, 81.16.Rf, 73.22.Pr

We obtained highly n-type doped graphene by intercalating terbium atoms between graphene and SiC(0001) through appropriate annealing in UHV. After terbium intercalation angle-resolvedphotoelectron spectroscopy (ARPES) showed a drastic change in the band structure around the K-points of the Brillouin zone: the well-known conical dispersion band of a graphene monolayer was superposed by a second conical dispersion band of a graphene monolayer with an electron density reaching 10 15 cm -2 . In addition, we demonstrate that atom intercalation proceeds either below the buffer layer or between the buffer layer and monolayer graphene. The intercalation of terbium below a pure buffer layer led to the formation of a highly n-doped graphene monolayer decoupled from the SiC substrate, as evidenced by ARPES and XPS measurements. The band structure of this highly n-doped monolayer graphene showed a kink (a deviation from the linear dispersion of the Dirac-cone) which has been associated to electron-phonon coupling constant one order of magnitude larger than those usually obtained for graphene with intercalated alkali metals.

Based on the possibility of intercalation of various atoms, strongly inspired by the results observed in graphite intercalation compounds (GICs) and motivated by possible superconductivity observed for such intercalated systems, graphene became the investigation playground for a large variety of alkali metals such as Li, K, Ca, Cs and Rb [3][4][5][6][7]. For graphene, particularly for epitaxial graphene on silicon carbide (SiC), the band dispersion after intercalation can be measured by ARPES. For most of the above mentioned elements a complex band-structure was found, showing a combination of a new linear highly n-doped dispersion and the dispersion of pristine monolayer graphene. The new dispersion curve was either attributed to a decoupling of the buffer layer or to the doping of monolayer graphene on top when no intercalation was expected [8]. However, no systematic studies were performed to prove experimentally these assumptions and to try to distinguish the contribution of each graphene layers and to isolate the highly n-doped graphene. Indeed, it is not really clear how the type of element modifies the Fermi level of the highly n-doped graphene, if the observed two bands correspond to the mixing of the band structures of two separated graphene layers or if it could be associated to a new type of stacking of graphene layers as proposed by Kim et al. [8]. In the case of the intercalation of alkali metals a kink (a deviation of the linear dispersion band near the Fermi level) was observed and discussed in term of electron-phonon coupling in the framework of BCS-theory of superconductivity.

For these elements, the strength of electron-phonon coupling was systematically studied and found to depend on the density of transferred electrons [10]. In our study we focus on much more rarely studied lanthanide elements. Previous studies have measured the band structure of graphene with intercalated Ytterbium [11] and more recently Ce [12]. In both cases a complex band structure is observed and in the case of Ytterbium a faint highly n-doped Dirac cone was also observed but a possible electron-phonon coupling was neither observed nor discussed.

In this work we examined intercalation of terbium atoms (Tb) atoms using epitaxial graphene on SiC(0001). Using X-ray photoelectron spectroscopy (XPS) and angleresolved photoelectron spectroscopy (ARPES), we showed the appearance of a highly ndoped graphene after intercalation of terbium atoms. A monolayer graphene (ML) on SiC(0001) consists in two graphene layers. One is partially covalently bonded to the upper silicon atoms of the SiC substrate. This layer is called the buffer layer (BL), which is covered with a true ML of graphene on top. After intercalation, ARPES spectra exhibited a complex band structure with contribution from a highly n-doped graphene layer. XPS measurements revealed that the intercalated Tb remained metallic, and suggested that it was protected against oxidation. We also observed no Tb-C bond formation. Intercalation occurred also below the buffer layer generating the highly n-doped graphene layer. In contrast to previous studies, we were able to clearly distinguish the origin of the two dispersion bands and showed that it is possible to form a single highly n-doped ML graphene. The observed strength of the electron-phonon coupling was one order of magnitude higher than expected for intercalated alkali metal with our doping level [10].

EXPERIMENTAL METHODS

Monolayer graphene samples were prepared by annealing of 6H-SiC(0001) at 1100 o C at a base pressure of 1×10 -10 mbar. By varying the annealing time and temperature, we were able to prepare samples with different ratio of graphene to buffer layer. The terbium atoms were deposited by electron-beam evaporation in ultra-high Vacuum (UHV) at a base pressure of 8×10 -10 mbar. The deposition rate was monitored by quartz-crystal microbalance.

We deposited a terbium layer of approximately 0.3nm. After deposition of terbium, samples were first annealed at 600 o C for 10 minutes. Subsequent annealing for 1 hour resulted in intercalation of terbium atoms. In-house physical characterizations were performed in UHV by low-energy electron diffraction (LEED), by scanning tunneling microscopy (STM) and by photoelectron spectroscopy techniques (XPS, ARPES).

LEED experiments were performed in situ with a Omicron SpectaLEED. STM experiments were performed in situ with a LT-STM from Omicron at 77 K at a base pressure in the 10 -11 mbar range. Images were acquired in constant current mode with bias voltage applied to the sample and employing chemically etched W tips. Photoelectron spectroscopy (XPS, ARPES) measurements were performed in situ using a VG Scienta R3000 spectrometer equipped with a hemispherical analyzer. An energy resolution of 16meV is reached with a pass energy of 20eV. For XPS we used a monochromatic X-ray source (AlKα, 1486.6 eV) and for ARPES a UV source HeII (40.8eV). For more clarity, periodicity expressed with respect to SiC(0001) and graphene will be referred in the following to -SiC and to -G respectively.

Once deposited, terbium atoms which are highly reactive get progressively oxidized, even in UHV as shown by the Tb3d core level peaks in figure 1. Just after deposition the Tb3d core level spectra displayed a metallic character, as shown Fig. 1 a). Two components of the spectra, located around 1241.8 eV and 1275.5 eV, were attributed to Tb3d5/2 and Tb3d3/2 with a spin orbit coupling of 33.7eV, respectively The metallic character was evidenced by multiplet-splitting in the photoemission process due to the unfilled shell at the 4f core level [13]. These multiplets were well defined at 1238.8 eV and 1249.4 eV as indicated by arrows on figure 1. The O1s component in the XPS spectra increased over time, multiplet splitting disappeared and an oxidized Tb3d peak emerged at 1250.9 eV as shown in figure 1 d). This is attributed to the effect of hydridization between the final state 4f and O2p core level [14]. After annealing at 600 o C, the amount of oxide was reduced and the Tb3d core level peak showed again a more pronounced metallic character (Fig. 1 b)). However, this Tb layer was not stable in time and oxidized again even in UHV. After annealing at 800 o C (figure 1 c)), contribution from oxides were completely removed. The resulting structure was stable in UHV. This suggests that the Tb atoms were intercalated below the graphene layer and so were protected against oxidation. Figure 2 we deconvoluted all components in the spectra resulting from strong multiplet-splitting and C1s, respecting the atomic ratio. After Tb intercalation it was not possible to identify the components S1-sp3 and S2-sp2 of the BL. However, a new component labelled "NewBL" emerged at 285.16 eV suggesting the formation of a new layer resulting from the decoupling of the BL from SiC after the intercalation of Tb atoms.

Figure 3 displays the STM and LEED measurements done on pristine epitaxial ML graphene (fig. 3 A)) and after intercalation of terbium atoms (after deposition and annealing at 800 o C, Figs. 3 B) to F)). The pristine graphene consists in large areas of ML graphene as schematized in figure 2 C) with some part of buffer layer graphene (typically less than few percent), which appears darker in the figure (labeled BL). The LEED pattern in the inset shows the typical 1x1-G diagram surrounded by the so-called 6x6-SiC reconstruction which is due to the BL reconstruction under the ML graphene. After intercalation of terbium, terraces are no more well defined but seems to be indented-like with small recessed areas observed everywhere (less than 10% of the surface, encircled in fig. 3 B)), at the former positions of BL areas and even in the middle of terraces. Figures 3 C) andD) are zoom on top of a terrace and on such recessed area respectively. Corresponding insets (Figs. 3 E) and F)), showing all 6 C atoms, confirm the presence of a graphene layer on top of intercalated structures, a layer which is more corrugated in the case of Fig. 3 D). In both cases the BL 6x6-SiC reconstruction has disappeared and intercalated terbium atoms (below ML and BL) are not ordered. The absence of underlying reconstruction confirms that the BL is completely decoupled from the substrate. This is proven as well by the LEED pattern in Fig. 3 B) which shows the 1x1-G pattern but the 6x6-SiC spots are no more visible.

ARPES spectra for pristine ML graphene (shown in Fig. 4 B)) taken along the ΓK direction as shown in the schematic in Figure 4 A), showed typical dispersion of π band with the Dirac point at 0.45 eV below the Fermi level due to n-doping from substrate [1] (ndoping = 1 × 10 13 cm -2 ). After deposition of terbium atoms, no dispersion was observed by ARPES around the K points, which demonstrates that terbium atoms completely covered the graphene layer. After annealing at 800 o C a modification in the band structure was observed. The ARPES spectra showed that the contribution from pristine graphene to the dispersion curve was then accompanied by a new band with a large Fermi surface, as shown in Figure 4 B). The n-doping and the Fermi surface of this new band increased when the sample was additionally annealed at 800 o C for 1 hour (see Fig. 4 C)). The resulting complex band structure can be interpreted as the superposition of two dispersion curves originated from two distinct graphene layers, as previously discussed for intercalation of ytterbium and aluminum [2,11] or for the upper graphene layer on top of bilayer graphene after doping [8]. The Fermi velocity (v F ) of pristine graphene was measured to be around 0.83 ± 0.07 × 10 6 m.s -1 . v F and E D2 of the new graphene layer observed after 10 min of annealing at 800 o C were 0.80 ± 0.09 × 10 6 m.s -1 and -1.57 eV, respectively. After annealing at 800 o C for 1 hour v F and E D2 changed to 0.65 ± 0.07 × 10 6 m.s -1 and -1.59 eV respectively. The density of electrons at the Fermi level can be calculated by the following relation between E D and v F as

The schemes in

n e = E 2 D πh 2 v 2 F ( 1
)
Using Equation (1) the electron density calculated for the pristine graphene is 2.6x10 13 cm -2 . The electron density for the new graphene is around 4.3 × 10 14 cm -2 .

An important question in this study concerns the origin of this new graphene layer. In earlier intercalation studies on graphene, highly n-doped spectra were always attributed to the decoupling of the buffer layer from the underlying SiC substrate. However, to the best of our knowledge, no clear experimental evidence was provided. The contributions of pristine ML graphene and highly n-doped graphene to the band structure have never been separated. In order to identify the exact nature of the component labelled "NewBL" we have generated a pure pristine buffer layer for which we studied the decoupling of the buffer layer resulting from the intercalation of terbium atoms. Samples of pristine buffer layer were prepared by annealing of SiC substrates at 900 o C-1000 o C in UHV at a base pressure of 1 x 10 -10 mbar. The corresponding C1s core level spectrum is given figure 5 A). We observed only three components, the SiC peak and components S1-sp3 and S2-sp2 corresponding to the BL. After the terbium intercalation, the XPS spectra revealed a clear modification of the buffer layer confirming the deconvolution done previously on the C1s core level peaks for the ML graphene after intercalation of terbium atoms. Figure 5 B) shows the C1s core level spectrum after intercalation of terbium atoms under buffer layer graphene. The peaks S1-sp3 and S2-sp2 corresponding to the BL have again disappeared and the peak called "NewBL" was clearly observable at 285 eV. A slight shift of around 0.36 eV toward lower energy was observed for the C1s peak of SiC. The full width at half maximum (FWHM) of this peak was 0.97eV, which was larger than one of the C-Sp2-G peak of the graphene layer (0.65eV). We attribute this increase in FWHM to inhomogeneities of the doping across the whole sample.

In order to clearly confirm that the new structure of the buffer layer corresponds to a highly n-doped graphene layer as observed in the previous complex band structure found by ARPES, we have systematically varied the graphene to buffer layer ratio in the pristine samples before terbium intercalation by adjusting the annealing conditions. We have prepared samples with ratios of graphene to buffer layer of 80/20, 50/50, 30/70 and 0/100 and performed intercalation of terbium atoms under the same conditions for every sample.

The graphene-to-buffer layer ratio can be verified by the C1s core level spectra, acquired on the samples before intercalation experiments, given in the figures 6 A) to D). The resulting ARPES curves are shown in the figures 6 E) to H).

The intensity of the new n-doped graphene dispersion band increased with increasing percentage of BL content on the sample while the contribution of the pristine graphene layer disappeared. Finally, for the sample with 100% buffer layer a single dispersion curve was observed, which corresponds to a highly n-doped graphene layer. This confirms that the origin of this highly n-doped band has to be attributed to the decoupling of the buffer layer caused by intercalation of Tb atoms. We show also that it is possible to synthesize an isolated highly n-doped graphene layer with a level of doping which can reach 10 15 cm -2 , a value similar to values obtained by ionic liquid doping.

Then a second question arises about the origin of the second band observed in Figs. 4 C) and D) and schematized by a black curve in Figs. 4 F) and G). Indeed, we may assume that it is due to the presence of areas of pristine monolayer graphene remaining after the intercalation of terbium atoms, as its ARPES spectrum looks like the one of pristine monolayer graphene seen in Fig. 4 B). But XPS and LEED results have proven that BL is decoupled everywhere by Tb atoms. In this case, the band structure in ARPES should be more parabolic (more massive fermions) and STM images at low bias should resolved only one per two C atoms of the top ML, as in the case of bilayer graphene in A-B stacking.

ARPES spectra shown in Figs. 4 C) and D) having linear dispersion and STM images, where every 6 C atom is imaged (see Figs. 3 B), C) and E), representing more than 90% of the surface), rule out this possibility. Moreover, disordered intercalated structures are clearly visible under ML. Our observations strongly support the conclusions of K. S. Kim et al. [8] and recent calculations done in the case of Yb intercalation by M. Kim et al. [START_REF] Kim | Manipulation of Dirac cones in intercalated epitaxial graphene[END_REF], where it is supposed that the intercalation between the bottom and top graphene layer perturbs the stacking leading to a A-A stacking. We must here conclude that the intercalation happens everywhere on the sample surface, under the buffer layer and also between the top and the bottom graphene layer.

In addition, a kink near the Fermi level was observed in the ARPES results. As shown in figure 3 C), it became more pronounced when both band structures were observed. Such a kink may reveal collective excitations and possible electron-phonon coupling. The electron phonon coupling strength (λ) is a key parameter in the determination of the critical temperature of superconductivity in the framework of the BCS theory, particularly for the graphite/graphene intercalation compounds [15]. The strength of electron-phonon coupling depends on the density of states at the Fermi level, on the phonon frequency mode and on how much the intercalated atoms are able to bend the in-plane graphene vibrations into out of plane graphene phonons mode (deformation potentials). In a systematic study of λ as a function of the intercalated alkali metals for graphene compounds, Fedorov et al. [10] found experimental evidence for a dependence of λ on the density of electron (D) at the Fermi level following the relation λ=10 -9 √ D with D in cm -2 . They found for the highest electron density 5.10 14 cm -2 in their system a value of λ=0.16. Considering the highly n-doped graphene layer, and in the hypothesis that the kink would be associated to e-phonon coupling and not other e-e correlation effect, we have extracted the electron-phonon coupling strength from ARPES measurements following the method used by Watchrinyanon et al. [10] and found λ=1.2. This value is one order of magnitude higher than the one expected in the case of alkali metals. Indeed for our level of electron doping reaching 10 15 cm -2 we should expect λ=0.25 [10]. As discussed in [10] such a high value of λ indicates either a possible ordering of the intercalated species, the involvement of more than one layer of graphene (this is the case for GiCs), or a reduced distance between the intercalated atoms and the graphene layer.

CONCLUSIONS

In conclusion, we performed the intercalation of terbium atoms under the epitaxial graphene on SiC(0001). By using ARPES, XPS and STM measurements, we have shown that terbium atoms can be intercalated under the buffer graphene layer which is transformed in a new graphene monolayer with a doping level reaching 10 15 cm -2 . This doping level is of the order of magnitude of the highest doping level achieved with ionic liquid gate field.

A strong electron-phonon coupling strength λ is deduced from APRES measurements.

Compared to other systems, and particularly to the intercalated alkali metals which have been systematically studied and for which a superconductivity has been observed, our value is one order of magnitude higher than expected for our high doping level and also higher than for graphite intercalation compounds. Of course the link with the superconductivity is still not done here but these value exemplarily demonstrate the interest of terbium as intercalated elements. Further work is required to improve the homogeneity of the sample and then perform transport measurements. The role of terbium coverage, the intercalation temperature and how the n-doping happens have to be studied systematically.

  scheme for ML graphene with the BL on SiC(0001) is shown in figure 2 (c). The typical C1s core level peak shows four components (figure 2 a)). One, labeled SiC, corresponds to carbon atoms in the SiC substrate at 283.76 eV and the C-sp2-G component at 284.68 eV corresponds to the top graphene monolayer. The BL is a graphene layer partially covalently bonded to the underlying SiC substrate and contributes to two components: S1-sp3 at 285.02 eV is attributed to covalently bonded carbon atoms, S2-sp2 at 285.62 eV to C-C bonds in the BL. After intercalation as shown in figure 2 (b), the C-sp2-G component did not change. However the contribution of the SiC component was shifted by 0.36 eV towards lower energies, attributed to the modification of the work function due to the presence of the intercalated terbium atoms on top of the SiC substrate. Such observations were also made for the intercalation of other elements [2, 11]. The contribution of Tb4p3/2 is at the same position as the one of C1s. Following the procedure proposed by Padalia et al. [13],

Figures 4 E

 4 ) to G) provide a simplistic interpretation of the ARPES spectra where the dispersion curve resulting from the initial graphene layer is shown in black. The Dirac point (E D1 ) of the initial graphene layer remained located at 0.45 eV below the Fermi level. The dispersion curve corresponding to the new graphene layer is represented by the red line. This new graphene layer had a large Fermi surface increasing with annealing time and a Dirac point (E D2 ) located 1.57 eV below Fermi level, as shown in Figures 4 F) to G).

FiguresFIG. 1 :

 1 Figures

FIG. 4 :FIG. 5 :

 45 FIG. 4: (A) Schemes of the Graphene reciprocal lattice and the K point where the ARPES spectra are acquired. ARPES spectra of (B) Pristine graphene and after deposition of Tb and annealing at 800 o C for 5 min (C) and for (D) 30 min. In (C) and (D), line scans along different energies (near E f and near -2.4 eV) are represented, in order to show clearly the existence of a new band structure, in particular in the lowest part of the ARPES spectra. The inset in D shows a short window represented by the dashed line rectangle in order to increase the contrast locally. (E) to (G) are the corresponding scheme showing an increase of the Fermi surface of the new band with annealing time, as evidenced by the line scan near E f in (D) where the green arrows represent the positions of the new band in the line scan near E f in (C).

FIG. 6 :

 6 FIG. 6: C1s core level spectra of pristine graphene samples with varying Graphene to Buffer layer ratio (A) 80/20 (B) 50/50 (C) 30/70 (D) 0/100. E)-H) resulting ARPES spectra after terbium intercalation on these samples. This clearly shows that the new highly n-doped band structure results from the decoupling of BL by intercalated terbium atoms.

This work was supported by the Région Alsace, the ANR ChimigraphN and the Université Franco-Allemande).