Aster Team

Georges Aumaître

Didier L Bourlès

Karim Keddadouche

Transient hillslope erosion in slow evolution landscapes

Transient evolution and adjustment to changing tectonic and climatic boundary conditions is an essential attribute of landscapes and characterizing transient behavior is a key to understand their dynamics and history. Developing new approaches to detect such transience has been explored by various methods, in particular to identify landscape response to Late Cenozoic and Quaternary climatic changes. Such studies have often focused on regions of high relief and/or active tectonic activity where interferences between tectonic and climatic signals might complicate the interpretation of the observations. We investigated the case of the hillslopes of the Serra do Cipó quartzitic range in SE Brazil in order to detect and quantify transience in a tectonically-quiescent landscape over 100 ka timescales. We determined hilltop curvature from a high resolution Digital Surface Model derived from Pléiades imagery and measured cosmogenic nuclide concentrations ( 10 Be and 26 Al) at these hilltop 1 sites. We compare both observations with predictions of hillslope diffusion theory, and observe a distinctive signature of an acceleration of denudation. We performed a joint inversion of topographic and isotopic data to retrieve an evolution of hillslope sediment transport coefficient through time. The timing of the increase in denudation cannot be unequivocally associated with a single climatic event but is consistent with climatically-modulated important fluctuations in precipitation and erosion in this area during Middle and Late Pleistocene.

Introduction

The Earth surface acts as a major interface between internal and external geodynamics, and is the focus of many investigations concerning the complex network of processes connecting surface processes, tectonics and climate [e.g. [START_REF] Champagnac | Tectonics, climate, and mountain topography[END_REF]. Notably, a key element of this system is the degree of coupling between climatic variations and denudation, which is fundamental to assess the strength of many postulated feedbacks between the atmosphere and the dynamics of weathering, sediment production and transport across continental surfaces. Over short time scales, many observations support the key role of temperature and precipitation conditions in directly controlling the modes and rates of weathering of materials transiting through regolith profiles, as well as their export across and beyond hillslopes domain [e.g. [START_REF] Gaillardet | Global silicate weathering and CO2 consumption rates deduced from the chemistry of large rivers[END_REF][START_REF] Hales | Climate-controlled variations in scree production, Southern Alps, New Zealand[END_REF][START_REF] Gabet | Modern erosion rates in the High Himalayas of Nepal[END_REF][START_REF] West | Aspect-dependent variations in regolith creep revealed by meteoric 10Be[END_REF]. However, when integrating these processes over larger spatial and temporal frames, highly contrasting propositions have been made concerning the actual relationship between surface processes acting on continental relief and climatic evolution, over the Late Cenozoic. In particular, an important question still actively debated concerns whether or not Late Cenozoic cooling was associated with a global increase in denudation rates [e.g. [START_REF] Herman | Plio-Pleistocene increase of erosion rates in mountain belts in response to climate change[END_REF][START_REF] Willenbring | The null hypothesis: globally steady rates of erosion, weathering fluxes and shelf sediment accumulation during Late Cenozoic mountain uplift and glaciation[END_REF][START_REF] Norton | Lack of a weathering signal with increased Cenozoic erosion?[END_REF], suggesting that additional investigations using novel methodologies are needed to better resolve denudation variations and transient landscape dynamics through time.

Regions of high relief or tectonic activity have often been the primary target of studies focusing on climate-denudation links, using the sedimentary archives resulting from the erosional fluxes in such areas [e.g. [START_REF] Charreau | Paleoerosion rates in Central Asia since 9Ma: A transient increase at the onset of Quaternary glaciations?[END_REF][START_REF] Lenard | Steady erosion rates in the Himalayas through late Cenozoic climatic changes[END_REF]. However, important limitations can be encountered in this type of settings, such as high denudation rates reaching the analytical limits of many measurement methods, or gaps in the sedimentary records and stochastic events introducing a high degree of variability in the denudation signal. Additionally, temporal and spatial variations of the tectonic boundary condition can be difficult to deconvolve from the record and complicate interpretations in terms of climatic influence on surface processes.

Conversely, areas of low tectonic activity may provide propitious conditions allowing the investigation of the coupling between climate variations and surface processes. However, long-duration continental sedimentary archives, spanning several 100s ka or Ma, are often scarce in such settings so alternative approaches are required to develop longer-term records of landscape transience and assess its link with climatic evolution [START_REF] Struck | Soil production and transport on postorogenic desert hillslopes quantified with 10Be and 26Al[END_REF]. When measuring active denudation, as opposed to paleo-denudation records in sedimentary archives, slow landscape evolution also implies long integration timescales for most geochronological techniques, such as cosmogenic nuclides, which might provide insights on the dynamics of surface processes over 100s ka to several Ma. This time frame is at least one order of magnitude longer to that usually encountered in active orogenic settings.

Terrestrial Cosmogenic Nuclides (TCN) are one of the most widely used tools to constrain the distribution of denudation rates at the Earth surface over 1 ka-1 Ma timescales [START_REF] Von Blanckenburg | The control mechanisms of erosion and weathering at basin scale from cosmogenic nuclides in river sediment[END_REF]. Attempts to recover information on past denudation conditions have mostly focused on TCN time-series in the sedimentary record. However, combinations of several cosmogenic nuclides, with different decay constants, offer the possibility to identify transient dynamics resulting from changes in denudation rates through time [START_REF] Mudd | Detection of transience in eroding landscapes[END_REF].

For example, the promising 10 Be-14 C pair has been used to decipher complex dynamics of surface processes over Holocene time scales [START_REF] Hippe | Constraining processes of landscape change with combined in situ cosmogenic 14 C-10 Be analysis[END_REF][START_REF] Skov | Detecting landscape transience with in situ cosmogenic 14C and 10Be[END_REF]. The more common 26 Al-10 Be system, which is classically used for burial dating studies, is much less sensitive to fluctuations of denudation on this time frame [START_REF] Mudd | Detection of transience in eroding landscapes[END_REF], but has been shown to offer important insights over longer time scales [Knudsen and Egholm, 2018;[START_REF] Struck | Soil production and transport on postorogenic desert hillslopes quantified with 10Be and 26Al[END_REF], and its actual potential to decipher transient denudation evolution in slow denudation landscapes needs to be explored more systematically. Independent of chronological approaches, another way to characterize transience in landscapes, resulting from changes in denudation or incision processes, is the study of the imprint of such perturbations on their morphology. Such approaches has been intensively used in fluvial systems , with the analysis of river long-profiles [e.g. [START_REF] Roberts | Estimating uplift rate histories from river profiles using African examples[END_REF][START_REF] Fox | Abrupt changes in the rate of Andean Plateau uplift from reversible jump Markov Chain Monte Carlo inversion of river profiles[END_REF], and the increasing availability of high resolution topographic data has paved the way for similar developments in the hillslope domain [START_REF] Roering | Functional relationships between denudation and hillslope form and relief[END_REF][START_REF] Hurst | Using hilltop curvature to derive the spatial distribution of erosion rates[END_REF][START_REF] Hurst | Hillslopes record the growth and decay of landscapes[END_REF]Grieve et al., 2016b;[START_REF] Mudd | Detection of transience in eroding landscapes[END_REF]. Notably, the combination of TCN data with high resolution hillslope morphological analysis has proven to be a fruitful source of information on landscape dynamics and their response to past perturbations [START_REF] Godard | Transient landscape dynamics across the Southeastern Australian Escarpment[END_REF].

In this study, we investigated the slowly eroding quartzitic landscapes of the Serra do Cipó Range in Minas Gerais, Brazil, by performing both high resolution morphometric analysis and cosmogenic nuclides concentration measurements ( 10 Be and 26 Al) at hilltop sites. The observed relationships between these parameters imply an important acceleration of denudation over 100s ka timescales and illustrate a new approach to probe into the transient dynamics of landscapes as well as expanding our knowledge of Quaternary paleoclimate of the central-eastern portion of South America.

Theoretical background

We present here the key concepts and ideas supporting our analysis of transient hillslope evolution in the Serra do Cipó. As reviewed by [START_REF] Mudd | Detection of transience in eroding landscapes[END_REF], transient dynamics is a fundamental attribute of landscapes involving changes in their morphology, sediment fluxes or rates of denudation, and diagnostic observations can be made on these properties to identify such transient state in terms of deviation from morphological or erosional steady states. In the case of hillslopes, insights can be obtained from high resolution morphometric analysis of their shape or from denudation rates calculated from Terrestrial Cosmogenic Nuclides (TCN) concentration measurements in bedrock or stream sediments, which are two types of observations we combine in our study. [START_REF] Roering | Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology[END_REF] proposed an hillslope transport model, where sediment flux q s is non-linearly dependant on topographic gradient dz/dx,

Hillslope erosion formulation

q s = -K(dz/dx) 1 -[(dz/dx)/S c ] 2 , ( 1
)
where K is a transport efficiency, or diffusion, coefficient, z is the elevation of the topographic surface, x the horizontal coordinate and S c a threshold gradient. Numerous field investigations have supported the validity of this formulation in high relief landscapes [e.g. [START_REF] Ouimet | Beyond threshold hillslopes: Channel adjustment to base-level fall in tectonically active mountain ranges[END_REF][START_REF] Dibiase | Landscape form and millennial erosion rates in the San Gabriel Mountains, CA[END_REF]. In part of the landscape where topographic gradients are small compared to the threshold S c this relationship can simplify as,

q s = -K dz dx , (2)
Equations 1 and 2 are Geomorphic Transport Laws [START_REF] Dietrich | Geomorphic transport laws for predicting landscape form and dynamics[END_REF], which can be combined with a continuity statement,

dq s dx = βE, (3)
with β the rock-to-regolith density ratio and E the long-term denudation rate, which is equal to rock uplift rate when steady state is achieved. Combining equations 2 and 3 yields,

E = KC HT β , (4)
where C HT is hilltop curvature, the second-derivative of the topographic surface. Equation 4connects together two fundamental and distinct types of observations, which can be made about landscapes: measurements of their spatial topographic properties, in this case C HT , and measurements of their temporal rate of evolution, in this case erosion rate E.

For a given landscape, assessing whether these two types of observations can be reconciled within the framework of equation 4 can be used as a test of the underlying steady-state assumption. It should be noted that equation 4 was derived for the particular case of a ridge-like hillslope system. In the following, we will rather consider isolated hills and a radial symmetry for the hillslopes, which implies a slight modification of equation 4 [START_REF] Anderson | Modeling the tor-dotted crests, bedrock edges, and parabolic profiles of high alpine surfaces of the Wind River Range, Wyoming[END_REF],

E = 2 KC HT β .
(5)

Detection of transience from the confrontation of morphometric and TCN data

Terrestrial Cosmogenic Nuclides (TCN) are now routinely used for the determination of denudation rates in a variety of settings, as their steady-state concentration C in various materials, such as stream sediment or outcropping bedrock, can be related to denudation rate E with,

C = i P i Eρ Λ i + λ , (6)
with P i the TCN production rate at the surface, ρ the material density, Λ i the cosmic ray attenuation length and λ the decay constant for the considered nuclide. Subscripts i refers to the various production pathways, such as spallation by neutrons or muons capture.

The integration time scale of these denudation rate measurements is usually calculated as the time required to erode a thickness equivalent to the neutrons attenuation lengthscale, or ∼60 cm in most situations [START_REF] Von Blanckenburg | The control mechanisms of erosion and weathering at basin scale from cosmogenic nuclides in river sediment[END_REF]. It results in a wide range of integration windows from Ma to ka timescales for denudation rates ranging from 1 to 100 m/Ma, respectively. This dependency of integration time on denudation rates implies that sites with different denudation rates will respond more or less rapidly, in term of TCN-derived denudation rates, to perturbations of similar amplitude, and can be used to identify transient adjustment in landscapes as we illustrate below.

If a local cosmogenic nuclides steady state is achieved at the hilltop scale, the relationship between denudation rate E and hilltop curvature C HT is linear and set by equation 5, assuming a spatially constant value for the sediment transport coefficient K (Fig. 1).

Transience resulting from a change in the sediment transport coefficient through time can be detected by the non-alignment of the E vs C HT or [ 10 Be] vs C HT relationships to theoretical steady-state predictions. For the case of a slow denudation landscape, we assume that hillslope morphology does not change significantly over a 100 ka time-scale, and do not take into account eventual feedback between denudation and hilltop curvature. Be-derived denudation rate E, as a response to a 5-fold increase in sediment transport coefficient K. Thick grey lines are the starting and final steady-state situations. Thin grey lines are steady state relations for intermediate sediment transport coefficients (equation 5). Red curves display the transient evolution at 100 ka intervals. (B) Same as previous panel but for 10 Be concentration against C HT (log axes). Note that we consider hillslope morphology, and thus hilltop curvature, to be stable over a time scale of a few 100s ka. This approximation is justified by the hillslope response time of the order of 1 Ma implied by the low sediment transport coefficient considered here.

Detection of transience from TCN pairs

As discussed by [START_REF] Mudd | Detection of transience in eroding landscapes[END_REF], independently of morphological considerations, transience in denudation can also be detected by the use of pairs of cosmogenic nuclides with different decay constants. 10 Be-26 Al is the most widely used TCN pair, notably in burial dating applications [START_REF] Granger | Dating sediment burial with in situ-produced cosmogenic nuclides: theory, techniques, and limitations[END_REF]. However, their decay constant are of the same order of magnitude, which greatly limit the sensitivity of this system for the identification of transient adjustment conditions [START_REF] Mudd | Detection of transience in eroding landscapes[END_REF]. Instead, this has oriented research efforts toward the 10 Be-14 C pair as the most promising system to resolve recent changes in denudation [START_REF] Hippe | Constraining processes of landscape change with combined in situ cosmogenic 14 C-10 Be analysis[END_REF][START_REF] Skov | Detecting landscape transience with in situ cosmogenic 14C and 10Be[END_REF][START_REF] Hippe | Cosmogenic in situ 14 C-10 Be reveals abrupt Late Holocene soil loss in the Andean Altiplano[END_REF]. Nevertheless, it should be noted that in the particular case of slowly eroding landscapes, where the 10 Be and 26 Al half-lives are comparable with the TCN denudation integration timescale, the amplitude of the deviation of isotopic ratio away from steady-state denudation predictions can become analytically resolvable (Fig. 2). This deviation is particularly noticeable for an increase in denudation rate, as it significantly lowers the ratio below the steady-state denudation curve. The time-scale of the perturbation being also dependent on the denudation rate, the transient response can thus be observed over a few 100s ka in the case of slowly evolving landscapes (1 m/Ma).

Settings

Our study area is located in the Serra do Cipó national park, in the highlands of the Brazilian Southeast. The Serra do Cipó mountain is a narrow (30-km wide) quartzitic range located in the southern portion of the Espinhaço Range (Fig. 3), which is the most extensive and continuous Precambrian orogenic belt of the South American Platform and originated 600 Ma ago during the Brasilian Orogeny [START_REF] Junior | Southern Serra do Espinhaço: The Impressive Plateau of Quartzite Ridges[END_REF]. The Espinhaço Range is an important hydrographic and landscape divider in South America. Al surface concentrations and ratio to 2, 5 and 10fold increases or decreases in denudation rates. Left panels display the evolution through time of the 26 Al/ 10 Be ratio, since the onset of the perturbation. Right panels display the corresponding evolution in a two-nuclides plot. White and colored circles indicate the starting points and 0.5 Ma mark, respectivelly. Solid and dashed brown curves correspond to expectations for steady-state erosion and a constantly exposed surfaces, respectively. Production rates are calculated for the average conditions at the sampled sites in the Serra do Cipó (Fig. 3), using a time-independent scaling scheme [START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF].

To the west is located the large hydrographic basin of the São Francisco River, which drains a craton characterized by smooth relief, semi-humid climate, and savanna vegetation. To the east, drainage is organized around smaller watersheds, such as the Doce and Jequitinhona hydrographic basins, that drain a half-orange hills plateau-type landscape, under humid climate and tropical forests [START_REF] Salgado | The physical geography of Brazil[END_REF].

After its uplift during the Brasilian Orogeny, the Serra do Cipó has not been subjected to new orogenic phases and the subsequent erosion of the range resulted in a relief characterized by a wide quartzitic plateau (average elevation 1200 m), comprising a large sequence of ridges, reaching 2000 m in altitude, depressions and dissected fluvial valleys. There are no evidence of a direct impact of the opening of the South Atlantic on our study area.

Although landscape properties are globally homogeneous across the Serra do Cipó, the portion of the range that drains towards the Doce River is more dissected than that draining into the São Francisco River.

Nowadays, the climate of Serra do Cipó Range is mesothermal, with Mean Annual Temperature between 18 o C-19 o C and Mean Annual Precipitation ranging from 1250 to 1550 mm, as well as dry winters and wet summers. However, there are lines of evidence, from pedological and geomorphological records, that the region underwent many climatic changes during the Holocene and late Pleistocene, with occurrences of periods wetter and dryer than the present-day conditions [START_REF] Chueng | Landscape paleodynamics in siliciclastic domains with the use of phytoliths, sponge spicules and carbon isotopes: The case of southern Espinhaço Mountain Range, Minas Gerais, Brazil[END_REF][START_REF] Barros | Late quaternary landscape evolution in the Atlantic Plateau (Brazilian highlands): Tectonic and climatic implications of fluvial archives[END_REF][START_REF] Silva | Altitude, vegetation, paleoclimate, and radiocarbon age of the basal layer of peatlands of the Serra do Espinhaço Meridional, Brazil[END_REF].

Due to the great resistance of the quartzite bedrock to erosion, denudation rates are <6 m/Ma [START_REF] Barreto | Denudation rates of the Southern Espinhaço Range, Minas Gerais, Brazil, determined by in situ-produced cosmogenic beryllium-10[END_REF][START_REF] Peifer | Growing topography due to contrasting rock types in a tectonically dead landscape[END_REF], soils are thin and poor, hosting the typical vegetation of Campos Rupestres, a high altitude rocky field biome [START_REF] Junior | Southern Serra do Espinhaço: The Impressive Plateau of Quartzite Ridges[END_REF]. This combination of slow denudation rates and significant relief is well suited for the coupling of high-resolution morphometry and cosmogenic nuclides analyses due to long-integration time scales for denudation, a limited contribution of stochastic events to surface processes, a limited storage of erosion products and long relief response timescales.

Methods

General sampling strategy

Our strategy was to select prominent sites on top of soil-mantled vegetation-free hillslopes (Fig. 4), where the observed denudation rates E and morphological properties such as hilltop curvature C HT could be linked into the simple physical framework of equation 5 [START_REF] Anderson | Modeling the tor-dotted crests, bedrock edges, and parabolic profiles of high alpine surfaces of the Wind River Range, Wyoming[END_REF][START_REF] Hurst | Using hilltop curvature to derive the spatial distribution of erosion rates[END_REF]. The slope of this relationship between E and C HT is set by the sediment transport coefficient K and the rock-to-regolith ratio β.

We primarily selected regolith-mantled hilltop sites, with regular morphologies and no bedrock outcrops (Fig. 4A-B). We sampled quartz fragments resulting from the dismantlement of quartz veins embedded in the fine-grained quartzite bedrock and sandy regolith (Fig. 4C). At some sites, the quartz veins were directly outcropping, but not significantly protruding over the surface, and the quartz fragments were detached directly from the veins (Fig. 4D). At other sites, the veins were not directly outcropping, but covered by 10 to 20 cm of regolith. We observed that the quartz veins remain coherent and resistant during their exhumation toward the surface, whereas the fine grained quartzite bedrock is converted into sandy regolith. Site O is an exception, where we directly sampled the outcropping quartzite bedrock. O). Samples were crushed and sieved to extract the 250-1000 µm fraction, which was submitted to sequential magnetic separations. The remaining fraction was leached with 37%

Cosmogenic nuclides analysis

HCl to remove eventual carbonate fragments. The samples were then repetitively leached with H 2 SiF 6 and submitted to vigorous mechanical shaking until pure quartz was obtained.

Decontamination from atmospheric 10 Be was achieved by a series of three successive leachings in concentrated HF, each removing 10% of the remaining sample mass [START_REF] Brown | Examination of surface exposure ages of Antarctic moraines using in situ produced 10Be and 26Al[END_REF]. After addition of Be and Al carriers, the samples were digested in concentrated HF.

Be and Al were isolated for measurements using ion-exchange chromatography. 10 Be/ 9 Be and 26 Al/ 27 Al measurements were performed at the French AMS National Facility, located at CEREGE in Aix-en-Provence. All results and technical characteristics of the measurements and calculations are provided in tables 1 and 2.

For reference, and before investigating time-variations in denudation, we compute steadystate denudation rates (tables 1 and 2), which were calculated with the online calculator described in [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10Be and 26Al measurements[END_REF] and the nuclide specific LSD scaling scheme of Lifton et al.

[2014], using the CRONUS-Earth calibration dataset [START_REF] Borchers | Geological calibration of spallation production rates in the CRONUS-Earth project[END_REF] for the calculation of spallation production rates and muon production rates according to [START_REF] Balco | Production rate calculations for cosmic-ray-muon-produced 10Be and 26Al benchmarked against geological calibration data[END_REF].

We used 160 g/cm 2 for the effective neutron attenuation length in rock, and a density of 2.65 g/cm 3 . No shielding correction was considered for the hilltop sites we sampled.

Morphometric analysis

We acquired a triplet of Pléiades images covering our sampling area (Fig. 3B), and used photogrametric techniques to construct a 1 m resolution Digital Surface Model (DSM) and perform hilltop curvature calculations. Pléiades satellites provide panchromatic images with a 0.7 m resolution resampled at 0.5 m pixel size, and have the capacity to acquire multi-stereo views of a scene during the same orbit [START_REF] Bernard | 3D capabilities of PLEIADES satellite[END_REF]. In low vegetation settings, they provide an efficient and cost-effective alternative to airborne LiDAR to produce high-resolution Digital Elevation Models over large areas [e.g. [START_REF] Ansberque | The Longriqu fault zone, eastern Tibetan Plateau: Segmentation and Holocene behavior[END_REF]. While this type of imagery, and derived topographic data, has been extensively used in tectonic geomorphology investigations since the launch of the satellites in 2011, its use in quantitative process-oriented geomorphological studies is less common [START_REF] Godard | Transient landscape dynamics across the Southeastern Australian Escarpment[END_REF]. A triplet of Pléiades images was acquired by instrument PHR1A on the 5 th of July 2019 over a ∼250 km 2 area, covering the hilltop sites studied here.

The image triplet was processed with NASA AMES Stereo pipeline [START_REF] Broxton | The Ames Stereo Pipeline: Automated 3D Surface Reconstruction from Orbital Imagery[END_REF][START_REF] Moratto | Ames Stereo Pipeline, NASA's Open Source Automated Stereogrammetry Software[END_REF] in order to produce a Digital Surface Model (DSM). We used the Rational Polynom Coefficients (RPC) provided with the images for the internal calibration of the pushbroom sensor. No Ground Control Points were used during the processing. The images were first roughly aligned by orthorectification on WGS84 datum using a SRTM Digital Elevation Model. Low-resolution disparity maps for couples of images were computed to determine coarse correspondence between pixels. Sub-pixel refinement was then performed using an adaptive window correlator [START_REF] Nefian | A bayesian formulation for sub-pixel refinement in stereo orbital imagery[END_REF]. The disparity map was filtered by removing pixels with a value differing from the mean of their neighborhood by more than 3. The point cloud was triangulated and gridded to produce a 1-m resolution Digital Surface Model (Figs. 3 and5).

In order to compute hilltop curvature C HT , we fitted a quadratic surface to the DSM elevation over a circular window centered on the sampling point [START_REF] Hurst | Using hilltop curvature to derive the spatial distribution of erosion rates[END_REF]. The equation of the surface is,

z = ax 2 + by 2 + cxy + dx + ey + f. (7)
Topographic curvature C HT is computed as the two-dimensional Laplacian of the quadratic surface as,

C HT = 2a + 2b (8)
The radius of the analysis window was chosen to be large enough to filter short-wavelength surface roughness, and small enough to avoid perturbation from the ridge-and-valley topographic signal [START_REF] Lashermes | Channel network extraction from high resolution topography using wavelets[END_REF][START_REF] Roering | Evidence for biotic controls on topography and soil production[END_REF][START_REF] Hurst | Using hilltop curvature to derive the spatial distribution of erosion rates[END_REF][START_REF] Godard | Spatially averaged long-term erosion rates measured from in situ-produced cosmogenic nuclides in alluvial sediment[END_REF]Grieve et al., 2016a]. In our case, this radius was set at 10 m (Fig. 6). In order to account for uncertainties on the position of the sampling sites and the influence of small wavelength topographic fluctuations, we computed C HT at each pixel within a 10 m radius of the sampling point, and then calculated the mean and standard deviation of these values.

Numerical modeling

We predicted the concentrations in 10 Be and 26 Al at the surface by considering the progressive accumulation of these cosmogenic nuclides during exhumation as a response to imposed time-varying denudation, using a Lagrangian point of view [START_REF] Knudsen | Time-integrating cosmogenic nuclide inventories under the influence of variable erosion, exposure, and sediment mixing[END_REF].

We used the CRONUS-Earth calibration dataset [START_REF] Borchers | Geological calibration of spallation production rates in the CRONUS-Earth project[END_REF] and nuclides specific time-dependent LSD scaling scheme [START_REF] Lifton | Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes[END_REF], for the calculation of spallation production rates. Muons production was computed following Model 1A and associated parameters from [START_REF] Balco | Production rate calculations for cosmic-ray-muon-produced 10Be and 26Al benchmarked against geological calibration data[END_REF]. We used 160 g/cm 2 for the effective neutron attenuation length in rock, and a density of 2.65 g/cm 3 .

The originality of our approach is that the denudation rate E at our sites of interest is

• • • • • • • • • • • • • • • • • • • yαx -1.4 R 2 =0.94 yαx -1.1 R 2 =0.99
Figure 6: Evolution of the standard deviation of curvature as a function of the radius of the analysis window used to fit the quadratic surface of equation 7. 80×10 3 measurement points were sampled from the Pléiades DSM. A slight scaling break occurs at ∼10 m radius for the analysis window and this size is used in the curvature computation. We note that the scaling break in less pronounced than in the results of [START_REF] Roering | Evidence for biotic controls on topography and soil production[END_REF], due to the gentler variations in elevation of the Serra do Cipó hillslopes and the absence of trees, as opposed to the rugged and densely forested landscape of the Oregon Coast Range.

computed from their measured hilltop curvature C HT according to equation 5. This allows us to impose time-variations in denudation through changes in K, affecting all sites synchronously, rather than considering the evolution of denudation rates at each site individually. In other words, this allowed us to perform a joint inversion of a single K history to assess the changes in denudation through time, rather than dealing with the 11 independent denudation histories at each site separately.

We parameterized the evolution of the sediment transport coefficient K through time as a two-steps history, which is thus defined by three parameters : (1) the sediment transport coefficient K 1 during an initial slow denudation phase, (2) the sediment transport coefficient K 2 during a recent faster denudation phase and (3) the transition time T between the two phases. Similar scenarios of hillslope evolution based on a two-steps history for the sediment transport coefficient were explored by [START_REF] Roering | Constraining climatic controls on hillslope dynamics using a coupled model for the transport of soil and tracers: Application to loess-mantled hillslopes, South Island, New Zealand[END_REF] over shorter timescales, and allowed to capture the impact of a climate changing event while limiting the number of free parameters. A significant amount of scatter is expected in our dataset, which is usual in the type of complex natural system we are dealing with. Our forward model is based on a simplified definition of the problem, which does not account for numerous processes complicating hillslope dynamics and are not captured by idealized relationships such as equation 5. Additionally, the two-steps scenario we used is a coarse, but necessary, simplification of the actual denudation history. Such epistemic uncertainty is not represented by the measurement errors on either the 10 Be concentrations or hilltop curvatures. For that reason, we do not expect the modelled concentrations C mod to be able to fit the observed concentrations C obs within these measurement errors. We introduce a residual error σ e as a 4 th parameter which is sampled along the ones described above [START_REF] Laloy | Bayesian inversion of a CRN depth profile to infer Quaternary erosion of the northwestern Campine Plateau (NE Belgium)[END_REF], with a likelihood defined as,

L(K 1 , K 2 , T ) = 1 √ 2πσ e exp -1 2σ 2 e n i=1 (C obs,i -C mod,i (K 1 , K 2 , T )) 2 , ( 9
)
where n is the total number of observed TCN concentrations used in the inversion. We used a standard Monte Carlo-Markov Chain (MCMC) approach to move through the parameter space and estimate the posterior distributions [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF].

Data

At hilltops, the measured 10 Be concentrations range from 0.8±0.02 to 6.9±0.1 10 6 at/g (Fig. 7A and table 1), which corresponds to steady denudation rates of 6±0.5 to 0.35±0.05 m/Ma, respectively. The 26 Al/ 10 Be ratios range from 4.4±0.1 to 6.4±0.2, with 6 samples plotting within the erosion island, slightly above the steady state erosion curve (Fig. 7A and table 2). These samples are evenly distributed over the range of observed 10 Be concentrations. One sample (R) presents an anomalously high ratio, which is impossible to obtain by TCN accumulation and decay. It is probably affected by an analytical issue on the determination of the 27 Al concentration, which is significantly higher than the rest of the dataset (except sample O). Only the 10 Be concentration will be used for this sample in the following. Finally, 4 samples display ratios clearly below steady-state denudation predictions, including the 3 samples with the lowest 10 Be concentrations.

At our sampling sites, measured hilltop curvatures span an order of magnitude from 0.005 to 0.039 m -1 (Fig. 7B andtable 3). Denudation rates, calculated under a steady state assumption, increase with hilltop curvature. We note that sample O is slightly off the general trend and corresponds to the only site where we sampled the outcropping quartzite bedrock, instead of the quartz veins. 6 Discussion

Interpretation of observed TCN concentrations

The denudation rates calculated under a steady-state assumption span a large range, from 0.5 to 6 m/Ma. [START_REF] Barreto | Denudation rates of the Southern Espinhaço Range, Minas Gerais, Brazil, determined by in situ-produced cosmogenic beryllium-10[END_REF] reported denudation rates of 1.8 to 2.4 m/Ma and 4.0 to 6.3 m/Ma for catchments located in the São Francisco and Doce river basins, respectively.

The global range of denudation rates reported for these catchment wide denudation rates is largely overlapping with our local measurements near the drainage divide. We do not observe systematic gradient in denudation rates at our sites between the two sides of the divide, as in the results of [START_REF] Barreto | Denudation rates of the Southern Espinhaço Range, Minas Gerais, Brazil, determined by in situ-produced cosmogenic beryllium-10[END_REF], and the diversity of denudation rates in our dataset probably reflects an underlying long-standing diversity in hillslope and hilltop morphologies across the relict landscape of the Serra do Cipó uplands. Some of our samples present 26 Al/ 10 Be ratios that are significantly lower than steady-state denudation predictions (sites Q, M, U and V). Such low ratios are usually interpreted as reflecting temporary burial of the material [START_REF] Granger | Dating sediment burial with in situ-produced cosmogenic nuclides: theory, techniques, and limitations[END_REF] where, in the absence of TCN production, the differential rate of decay between the two isotopes is progressively lowering the ratio. However, such explanation is highly unlikely for our samples collected at the top of continuously eroding prominent hilltop sites, where there is no evidence for any significant potential for long-term storage and burial. Regolith mixing in slowly eroding soil profiles can also induce a lowering of 26 Al/ 10 Be ratios below the steady-state denudation line [e.g. [START_REF] Makhubela | Effects of long soil surface residence times on apparent cosmogenic nuclide denudation rates and burial ages in the Cradle of Humankind, South Africa[END_REF]. While this process could explain the measured ratio at site Q, it is difficult to reconcile with field observations at the other sites where we directly sampled outcropping veins. Furthermore, these three sites (M, U and V) display the highest denudation rates, and are eroding almost an order of magnitude faster than the slowest sites of the dataset. The last possible process to explain the observed low 26 Al/ 10 Be ratios at sites M, U and V is an increase of denudation rate [START_REF] Struck | Soil production and transport on postorogenic desert hillslopes quantified with 10Be and 26Al[END_REF], which would dominantly impact the fastest sites with short integration time scales (100-200 ka, table 1) and much less the slowest sites, which average denudation over a longer time-frame (Fig. 2). Independent of any consideration on the morphology of the studied sites, our TCN data present evidence of transience and acceleration of denudation rates in the past. The amplitude and timing of such increase is however difficult to constrain from

26 Al-10 Be data alone.

Relationship between TCN data and hilltop morphology

The relationship between denudation rate E and hilltop curvature C HT (equation 5) is not linear (Fig. 7B), as would be expected from hillslope diffusion theory [START_REF] Anderson | Modeling the tor-dotted crests, bedrock edges, and parabolic profiles of high alpine surfaces of the Wind River Range, Wyoming[END_REF][START_REF] Hurst | Using hilltop curvature to derive the spatial distribution of erosion rates[END_REF]. No single sediment transport coefficient can explain the data within uncertainty, and we observe a systematic tendency for the sites with higher denudation rates (lower integration time scale) to plot on a higher K line, when compared to the slower ones. In line with our previous observations on the 26 Al/ 10 Be ratios, this constitutes an additional evidence for an acceleration of denudation over our area of interest (Fig. 1). We note that a recent large-scale compilation by [START_REF] Gabet | Hilltop Curvature Increases with the Square Root of Erosion Rate[END_REF] have also identified nonlinearity in the relationship between C HT and TCN-derived denudation rates but further investigations are needed to assess its signification in terms of hillslope dynamics.

In order to obtain more quantitative insights on this postulated acceleration of denudation, we tested a 2-stage denudation history defined by 3 parameters : the transition time between the slow and fast periods of denudation and their respective sediment transport coefficients (Fig. 8). We predicted 10 Be concentrations at the sites of interest for any combinations of these parameters, using the measured hilltop curvature (C HT ) to compute denudation rates evolution through time using equation 5, and perform a MCMC driven inversion. The sampling of the parameter space is driven by the Metropolis-Hasting algorithm, with an acceptance rate of 32%. We ran 16 independent chains, each 5 × 10 5 , in length with a 10 4 length burn-in phase. The multivariate Potential Scale Reduction Factor is 1.001 suggesting convergence was achieved [START_REF] Gelman | Inference from Iterative Simulation Using Multiple Sequences[END_REF].

The inversion results support a low initial K value (< 8 × 10 -5 m 2 /a), followed by a very marked increase (5 × 10 -5 to 5 × 10 -4 m 2 /a). A pronounced trade-off exists between the onset time of the second phase and the value of its sediment transport coefficient, which makes it difficult to directly confront these results with climatological and morphological records. Indeed, both an intense and recent phase with a transition as young as 10 ka or a much older up to 1 Ma event with lower sediment transport coefficient could provide reasonable fits to the data (Fig. 8B). We note that solving for TCN accumulation using an simplified Eulerian approach [START_REF] Knudsen | Time-integrating cosmogenic nuclide inventories under the influence of variable erosion, exposure, and sediment mixing[END_REF] with exponential attenuation profile for muonic production and a time-independent scaling [START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF], instead of the Lagrangian scheme we used in our reference model, does not substantially modify the outcome of the inversion (Fig. 9B).

Additional independent constraints can be obtained from the low 26 Al/ 10 Be ratios of the fastest denudation sites Assuming an initial K 1 = 4 × 10 -5 m 2 /a, for each value of K 2 we computed the best transition time T allowing to explain both 10 Be and 26 Al concentrations within 1σ uncertainty, for each sample independently (Fig. 8A, red lines). For each of these samples (M, U, V), we observed the same tradeoff between K 2 and T as for the C HT -10 Be inversion of the whole dataset. The lines defining the best K 2 -T combinations for these 26 Al-10 Be data appears to intersect the region of maximum posterior density for the C HT -10 Be inversion at a transition time T ∼ 100 ka. Finally, a joint inversion of C HT , 10 Be and 26 Al data for the whole dataset narrows down the range of probable transition time and is also consistent with a value close to 100 ka (Fig. 9C). As discussed below, this time scale for a postulated shift in surface dynamics is in agreement with reported environmental changes in the area, notably precipitation [e.g. Cruz et al., 2009b].

We note that our results could be affected by the processes occurring during the exhumation of the quartz fragments sampled at the surface. Notably, a mixing of these fragments in the regolith layer, rather than the monotonous exhumation toward the surface assumed in our models, would significantly impact their accumulation of cosmogenic nuclides through time [Granger et al., 1996;[START_REF] Makhubela | Effects of long soil surface residence times on apparent cosmogenic nuclide denudation rates and burial ages in the Cradle of Humankind, South Africa[END_REF]. As discussed above, we have numerous lines of evidence that those fragments result from the dismantlement of massive quartz veins, which remained mostly cohesive during their exhumation history embedded in the regolith, implying little mixing. The 26 Al/ 10 Be ratios of most of the high-concentration samples are not compatible with a significant mixing of the fragments in the regolith (Fig 7A ), with the 8). Result of a MCMC inversion based on a 2-stages history set by 3 parameters: the initial (K 1 ) and final (K 2 ) sediment transport coefficients and the transition time T between the two. Blue shaded contours are highest posterior density region for T against K 2 . Solid lines are marginal probability density for K 2 and T , dashed line corresponds to K 1 . (B) Same as reference model but forward problem is solved with an Eulerian approach [START_REF] Knudsen | Time-integrating cosmogenic nuclide inventories under the influence of variable erosion, exposure, and sediment mixing[END_REF], using the time-independent [START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF] scaling and an exponential model for the muonic production at depth, with parameters from [START_REF] Braucher | Production of cosmogenic radionuclides at great depth: A multi element approach[END_REF]. This simplified and computationally much faster Eulerian scheme provides essentially the same results as the reference model. (C) Same as reference model but using both 10 Be and 26 Al concentrations (only 10 Be data is used in the reference model). (D) Same as reference model, but the density decreases to 1.7 during the last 1 m of the exhumation history, to test the influence of density change associated with a regolith layer. Response time τ (a)

E = 1 0 m /M a E = 1 m / M a E= 0. 1 m /M a L = 1 5 0 m L = 2 0 0 m L = 1 0 0 m l i n e a r
Figure 10: Adjustment time-scales of hillslopes based on linear (dashed lines) and nonlinear (solid lines) evolution models [START_REF] Roering | Hillslope evolution by nonlinear, slope-dependent transport: Steady state morphology and equilibrium adjustment timescales[END_REF][START_REF] Hurst | Using hilltop curvature to derive the spatial distribution of erosion rates[END_REF], for different values of the erosion rate E and hillslope length L.

exception of sample Q. We also note that our reference model is designed with a constant density during the exhumation history. In order to test the influence of a change in density when moving through the regolith for the final stage of exhumation [START_REF] Struck | Soil production and transport on postorogenic desert hillslopes quantified with 10Be and 26Al[END_REF],

we tested the effect of a 1-m thick layer of density 1.7, and we observed that such modification does not substantially alter the general pattern retrieved from the inversion for the various parameters (Fig. 9D). It is important to note that our approach relies on the hypothesis that hillslope morphology does not change substantially over the investigated time period, so that there is no feedback between denudation and hilltop curvature, which is considered constant here. Calculated hillslope response times for the range of very low denudation rates and sediment transport coefficient considered here are at least several Ma, supporting the hypothesis of a stable landscape over the time-frame of the climatic perturbations scenarios we investigate (Fig. 10).

Comparison with regional climatic and geomorphic records

The two types of observations we made (C HT -10 Be relationship and 26 Al/ 10 Be ratios) suggest a significant increase of denudation occurring over the last few 100s ka. However, it is difficult to pinpoint an actual event from each approach alone. Indeed, onset of rapid denudation as young as a few 10s ka, close to the Last Glacial Maximum, or as old as 800 ka, near the Mid-Pleistocene Transition, provide comparable fits to the C HT -10 Be data (Fig. 8). We note that we also used a very simplified 3-parameters definition of the evolution of K, which implies a single step-change and does not allow the possibility of multiple events or of a gradual change. We tested this type of scenario by mapping the evolution of the sediment transport coefficient on the δ 18 O stack of [START_REF] Lisiecki | A Pliocene-Pleistocene stack of 57 globally distributed benthic d18O records[END_REF], with a prescribed range of variations. Such evolution of K yields a relationship between C HT and 10 Be concentrations that does not fit the data as well as the abrupt evolution associated with the 3-parameters scenarios. However, it still reproduces the general trend and in particular the fact that it does not correspond to a steady state situation (Fig. 8C, green line). When considered alone, the 26 Al/ 10 Be ratios of the high-denudation samples also display a similar tradeoff between T and K 2 , but the combination with C HT -10 Be inversion results significantly narrows down the range of likely age for the onset of the increase in denudation rates to 100-300 ka (Figs. 8 andS4).

There are lines of evidence supporting important Quaternary fluctuations of precipitation, temperature and vegetation in SE Brazil [START_REF] Salgado-Labouriau | A dry climatic event during the late Quaternary of tropical Brazil[END_REF][START_REF] Behling | South and southeast Brazilian grasslands during Late Quaternary times: A synthesis[END_REF][START_REF] Ledru | Paleoclimate changes during the last 100,000 yr from a record in the Brazilian Atlantic rainforest region and interhemispheric comparison[END_REF], and in particular in the Espinhaço Range [START_REF] Chueng | Landscape paleodynamics in siliciclastic domains with the use of phytoliths, sponge spicules and carbon isotopes: The case of southern Espinhaço Mountain Range, Minas Gerais, Brazil[END_REF][START_REF] Silva | Altitude, vegetation, paleoclimate, and radiocarbon age of the basal layer of peatlands of the Serra do Espinhaço Meridional, Brazil[END_REF]. For example, precipitation has been shown to be strongly modulated by orbital fluctuations [START_REF] Wang | Millennial-scale precipitation changes in southern Brazil over the past 90,000 years[END_REF]Cruz et al., 2009a], with some speleothem records suggesting a wetter climate starting at 70 ka [Cruz et al., 2009b], a timing comparable with the transition time inferred from our inversion results. Early and Middle Quaternary geomorphological markers, such as river terraces, are often poorly preserved in the area, but display a strong climatic control on cycles of erosion and deposition [START_REF] Barros | Late quaternary landscape evolution in the Atlantic Plateau (Brazilian highlands): Tectonic and climatic implications of fluvial archives[END_REF]. The hypothesis that climatic fluctuations is the main driver of the increase in denudation rates is reinforced by the fact that there was a reduction in in-traplate tectonic activity beginning in the Middle Quaternary or earlier in the southeastern highlands of Brazil [START_REF] Salgado | Relief evolution of the Continental Rift of Southeast Brazil revealed by in situ-produced 10Be concentrations in river-borne sediments[END_REF]. This idea of a dominant climatic control is further reinforced by the non-linear E-C HT relationship observed in figure 7, which can be considered as symptomatic of a top-down climatic forcing, where climatically-triggered changes in precipitation or vegetation modulate hillslope sediment transport coefficient. Conversely, a bottom-up perturbation associated with tectonically-driven base level changes would preserve the co-linearity of the samples in an E-C HT plot. While we cannot unambiguously associate a particular event to our observations of transience in the landscape of the Serra do Cipó, they probably result from the integration of climatic fluctuations and an increased denudation trend over the Late Pleistocene.

Implications for the detection of transient evolution in landscapes

Combining the analysis of detailed high resolution hillslope morphometry with TCN data is a promising way to investigate landscape dynamics but has been only used in a limited number of situations [START_REF] Dibiase | Hillslope response to tectonic forcing in threshold landscapes[END_REF][START_REF] Hurst | Using hilltop curvature to derive the spatial distribution of erosion rates[END_REF][START_REF] Godard | Transient landscape dynamics across the Southeastern Australian Escarpment[END_REF][START_REF] Godard | Hillslope denudation and morphologic response to a rock uplift gradient[END_REF].

Such an approach allows to combine the two types of observation into a physical, processbased framework, such as equation 5, in order to test the underlying assumption and evaluate its limitations. Our study illustrates the potential of such an approach to identify transient dynamics in slowly eroding landscapes, and assess its links with climaticallydriven top-down forcings. We take advantage of the differences in integration timescales of the two methods: a few 100 ka for TCN concentrations and of the order of 1 Ma for hilltop curvature measurements (Fig. 10), which allows to unravel time variations in the intensity of surface processes. Hilltops of the Serra do Cipó range display several evidences of Late Quaternary acceleration of erosion, as the relationship between 10 Be concentrations and hilltop morphology, as well as the 26 Al/ 10 Be ratios, both displaying significant deviation from steady state predictions. There are however limits to the resolving power of such combination with respect to the details of denudation history that can be retrieved, which are manifest in the observed trade-offs between parameters (Fig 8A and9). Using additional TCN, such 36 Cl or in situ 14 C might allow to extract more details about the denudation history. However, in situ 14 C has not yet been used to investigate such low denudation settings and modeling studies anticipate important limitations to its resolving power in such settings [START_REF] Skov | Detecting landscape transience with in situ cosmogenic 14C and 10Be[END_REF]. 36 Cl would require target minerals such has carbonates or feldspars which are not present in the setting we investigated, but might be feasible in granitic landscapes. While it does not provide a time-series for the evolution of denudation, our study illustrates the possibilities and potential of investigating paleo-denudation of slowly evolving landscapes in addition to the often explored tectonicaly active mountain ranges.

Conclusion

We studied the denudation of hilltop sites in the Serra do Cipó range, by collecting high resolution morphometric data and measuring 10 Be and 26 Al concentrations in surface materials. We compared our observations with prediction of simple hillslope diffusion theory for the relationship between morphological attributes and denudation rates derived from these TCN concentrations. We identified substantial deviation from steady state configuration, suggesting an important change in surface denudation over recent time periods, which results in the observed transient state. The originality of our approach is to conduct a joint inversion of both isotopic and morphometric data to derive the temporal evolution of the hillslope sediment transport coefficient. Our inversions suggest an important increase in the sediment transport coefficient and denudation rate over 100 ka timescale, which most likely results from a modulation of soil transport efficiency of climatic origin. However, a trade-off between the timing of the acceleration and denudation intensity limits our ability to precisely pin point a specific climatic event. Our study illustrates the strong potential of landscape analysis that link geochronological and topographic data into a single physically-based framework, as for example hillslope diffusion theory. In our case, it allows to identify transient dynamics at the scale of hillslope, and the persistence of instabilities, which is key for longer-wavelength landscape-scale analysis of the influence of climatic and tectonic forcings. Be in-house carrier mass, ∼150 µl at 3.025 x 10 -3 g/g [START_REF] Merchel | Towards more precise 10Be and 36Cl data from measurements at the 10-14 level: Influence of sample preparation[END_REF].

c 10

Be/ 9 Be ratios were calibrated against the STD-11 standard by using an assigned value of 1.191±0.013×10 -11 [START_REF] Braucher | Preparation of ASTER in-house 10Be/9Be standard solutions[END_REF].

d
Uncertainties are reported at the 1σ level.

e Uncertainties on isotopic ratios are calculated according to the standard error propagation method using the quadratic sum of the relative errors and include a conservative 0.5% external machine uncertainty [START_REF] Arnold | The French accelerator mass spectrometry facility ASTER: Improved performance and developments[END_REF], the uncertainty on the certified standard ratio, a 1σ uncertainty associated to the mean of the standard ratio measurements during the measurement cycles, a 1σ statistical error on counted events and the uncertainty associated with the chemical and analytical blank corrections.

f One process blank was treated and measured with our samples, yielding a 10

Be/ 9

Be ratio of 3.53±0.34×10 -15

. It corresponds to an upper 1σ bound of 120 ×10 3 10 Be atoms for the background level in this blank, which is at least 120 times lower than the number of 10 Be atoms in the dissolved sample masses (600 times lower in average over our dataset).

g Denudation rates were calculated with the online calculator described in [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10Be and 26Al measurements[END_REF] and the nuclide specific LSD scaling scheme [START_REF] Lifton | Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes[END_REF], using the CRONUS-Earth calibration dataset [START_REF] Borchers | Geological calibration of spallation production rates in the CRONUS-Earth project[END_REF] for the calculation of spallation production rates and muon production rates according to [START_REF] Balco | Production rate calculations for cosmic-ray-muon-produced 10Be and 26Al benchmarked against geological calibration data[END_REF]. We use 160 g/cm 2

for the effective neutron attenuation length in rock, and a density of 2.65 g/cm 3

. No shielding correction was considered for the hilltop sites we sampled.

h Numbers in parenthesis are the denudation rates computed with Lal/Stone time-independent scaling.

i Integration time scales for 10 Be denudation rates [START_REF] Von Blanckenburg | The control mechanisms of erosion and weathering at basin scale from cosmogenic nuclides in river sediment[END_REF]. Al/ 27 Al ratios were normalized to the in-house standard SM-Al-11 whose ratio of 7.401±0.064×10 -12 [START_REF] Arnold | The French accelerator mass spectrometry facility ASTER: Improved performance and developments[END_REF] has been cross-calibrated against primary standards from a round-robin exercise [START_REF] Merchel | First international 26Al interlaboratory comparison -Part I[END_REF].

d
Uncertainties are reported at the 1σ level.

Figure 1 :

 1 Figure1: Illustration of the manifestation of transience in landscapes from the confrontation of cosmogenic nuclide data with high resolution topographic metrics. (A) Evolution of the relationship between hilltop curvature (C HT ) and 10 Be-derived denudation rate E, as a response to a 5-fold increase in sediment transport coefficient K. Thick grey lines are the starting and final steady-state situations. Thin grey lines are steady state relations for intermediate sediment transport coefficients (equation 5). Red curves display the transient evolution at 100 ka intervals. (B) Same as previous panel but for10 Be concentration against C HT (log axes). Note that we consider hillslope morphology, and thus hilltop curvature, to be stable over a time scale of a few 100s ka. This approximation is justified by the hillslope response time of the order of 1 Ma implied by the low sediment transport coefficient considered here.

Figure 2 :

 2 Figure2: Responses of 10 Be and 26 Al surface concentrations and ratio to 2, 5 and 10fold increases or decreases in denudation rates. Left panels display the evolution through time of the 26 Al/ 10 Be ratio, since the onset of the perturbation. Right panels display the corresponding evolution in a two-nuclides plot. White and colored circles indicate the starting points and 0.5 Ma mark, respectivelly. Solid and dashed brown curves correspond to expectations for steady-state erosion and a constantly exposed surfaces, respectively. Production rates are calculated for the average conditions at the sampled sites in the Serra do Cipó (Fig.3), using a time-independent scaling scheme[START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF].

Samples for 10

 10 Figure 3: (A) Situation of the Southern Espinhaço Range in Minas Gerais State. (B) Location of sampled hilltop sites in the Serra do Cipó Range. Color scale indicates calculated steady-state denudation rates from 10 Be concentrations (m/Ma). Dashed rectangle is the extent of the Pléiades images triplet used to compute a high resolution Digital Surface Model (DSM). Pink line is the divide between the São Francisco and Doce river basins.

Figure 4 :

 4 Figure 4: Sampling sites pictures. (A) View on Serra do Cipó landscape, toward the East from site M. (B) View on Serra do Cipó landscape, toward the East from site T. (C) Quartz vein at hilltop sampling site N. (D) Surface of hilltop sampling site W, with quartz vein embedded in the regolith and small quartz fragments. (E) Hilltop sampling site M. (F) View toward site W.

Figure 5 :

 5 Figure 5: Shaded Pléiades Digital Surface Model (DSM) for 4 selected hilltop sites. Spacing of elevation contours is 1 m.

  Figure7: (A) Two-nuclides plot for the samples collected at hilltops (1σ confidence ellipses). See figure3for sites locations. Grey-filled area indicates the region between the steady-state surface denudation and constant exposure scenarios. Grey and black lines correspond to various regolith mixing-depth and denudation rates, respectively, according to the model of[START_REF] Foster | Hillslope lowering rates and mobile-regolith residence times from in situ and meteoric 10Be analysis[END_REF]. Scaling according to[START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF]. (B) Denudation rate calculated from measured 10 Be concentrations as a function of hilltop curvature computed from the Pléiades DSM. Grey lines correspond to equation 5 with various values of the sediment transport coefficient.

Figure 8 :

 8 Figure8: (A) Result of the MCMC inversion based on a 2-stages history set by 3 parameters: the initial (K 1 ) and final (K 2 ) sediment transport coefficients and the transition time T between the two. Blue shaded contours are highest posterior density region for T against K 2 , and the filled blue circle indicates the Maximum Aposteriori Probability (MAP). Solid lines are marginal probability densities for K 2 and T , dashed line corresponds to K 1 . Red lines correspond to parameters combination allowing to predict the 10 Be and 26 Al concentrations within 1σ of the measured values (K 1 is set at 4 × 10 -5 m 2 /a). Red filled circle is a manually set parameters combination. (B) Marine δ 18 O stack from[START_REF] Lisiecki | A Pliocene-Pleistocene stack of 57 globally distributed benthic d18O records[END_REF] (log time scale). (C) Comparison of hilltop curvature (C HT ) and 10 Be concentrations. Grey lines are predictions for different time-constant sediment transport coefficients scenarios. Blue and red curves correspond to the predictions associated with sets of parameters from the filled circles on panel A. Dark green curve corresponds to the predictions with a sediment transport coefficient evolution mapped on the δ 18 O curve, ranging from 10 -5 to 3 × 10 -4 m 2 /a.

Figure 9 :

 9 Figure9: Tests of the influence of various aspects of the modelling approach on the inverted distributions for K 1 , K 2 and T . (A) Reference model (same as Fig.8). Result of a MCMC inversion based on a 2-stages history set by 3 parameters: the initial (K 1 ) and final (K 2 ) sediment transport coefficients and the transition time T between the two. Blue shaded contours are highest posterior density region for T against K 2 . Solid lines are marginal probability density for K 2 and T , dashed line corresponds to K 1 . (B) Same as reference model but forward problem is solved with an Eulerian approach[START_REF] Knudsen | Time-integrating cosmogenic nuclide inventories under the influence of variable erosion, exposure, and sediment mixing[END_REF], using the time-independent[START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF] scaling and an exponential model for the muonic production at depth, with parameters from[START_REF] Braucher | Production of cosmogenic radionuclides at great depth: A multi element approach[END_REF]. This simplified and computationally much faster Eulerian scheme provides essentially the same results as the reference model. (C) Same as reference model but using both 10 Be and 26 Al concentrations (only 10 Be data is used in the reference model). (D) Same as reference model, but the density decreases to 1.7 during the last 1 m of the exhumation history, to test the influence of density change associated with a regolith layer.

eA

  process blank yielded a ratio of 2.01±0.45×10 -14 . It corresponds to an upper 1σ bound of 730 ×10 3 26 Al atoms for the background level in this blank, which is at least 130 times lower than the number of 26 Al atoms in the dissolved sample masses (610 times lower in average over our dataset) f Procedures for the calculation of denudation rates are identical to the ones described in table 1.g Integration time scales for 26 Al denudation rates[START_REF] Von Blanckenburg | The control mechanisms of erosion and weathering at basin scale from cosmogenic nuclides in river sediment[END_REF].

Table 1 :

 1 Cosmogenic nuclides data : 10

	Be results

a Dissolved pure quartz mass. b

Table 2 :

 2 Cosmogenic nuclides data : 26 Al results and comparison with 10 Naturally occurring Al measured by ICP-OES. Note the high concentration for sample O, which is the only one from the quartzite bedrock instead of the quartz veins.

	Be results from table 1

a b Al Carrier mass, 993±4 ppm concentration. c The measured 26

Table 3 :

 3 Site locations and hilltop curvature measurements

	Site Latitude Longitude Altitude Hilltop curvature
		( o )	( o )	(m)	(m -1 )
	M -19.2637 -43.5558	1313	0.028±0.010
	N	-19.2614 -43.5508	1354	0.020±0.008
	O -19.2612 -43.5508	1353	0.039±0.010
	P	-19.2361 -43.5058	1355	0.006±0.002
	Q -19.2301 -43.5051	1365	0.005±0.001
	R	-19.2119 -43.5027	1335	0.008±0.001
	S	-19.1594 -43.5156	1380	0.008±0.003
	T	-19.1582 -43.5174	1379	0.007±0.003
	U	-19.2206 -43.4989	1337	0.038±0.011
	V	-19.2218 -43.4968	1301	0.024±0.007
	W -19.2337 -43.5050	1357	0.014±0.003

Acknowledgments

This work was supported by CAPES-COFECUB (Research Project 869/15) and an IUF fellowship to VG. The 10 Be and 26 Al measurements were performed at the ASTER AMS national facility (CEREGE, Aix en Provence) which is supported by the INSU/CNRS, the ANR through the "Projets thématiques d'excellence" program for the "Equipements d'excellence" ASTER-CEREGE action and IRD. The Pléiades images used in this study were obtained through the DINAMIS framework and GEOSUD EQUIPEX. We thank Rodolphe Cattin for insightful discussions. Constructive comments by Simon Mudd, an anonymous reviewer and the associate editor greatly helped to improve this manuscript.

This article is dedicated to the memory of our friend and colleague Didier L. Bourlès who sadly passed away during its preparation.


Ying Wang

Laurent Pedesseau

Claudine Katan

Jacky Even

Oleg V Prezhdo

Sergei Tretiak

Dibyajyoti Ghosh

Amanda J Neukirch
email: ajneukirch@lanl.gov

Nonadiabatic molecular dynamics analysis of hybrid Dion-Jacobson 2D lead Iodide perovskites

The past 6 years have witnessed the rapid growth of interest in Dion-Jacobson (DJ) phase two dimensional (2D) hybrid halide perovskites as optoelectronic materials with considerable intrinsic stability. The precise relationships between structural variations and the resulting charge carrier dynamics at finite temperature in these materials are key to practical applications and are not yet completely understood. Here, we study 3-(aminomethyl) piperidinium (3AMP) and 4-(aminomethyl) piperidinium (4AMP) spacer cation-based lead iodide DJ-phase systems and find these spacer cations to have a profound impact on the structural dynamics. Particularly, large conformational dynamics of the 3AMP based perovskite compared to that of the 4AMP at room temperature leads to pronounced state energy fluctuation near band edges and further results in a shorter quantum coherence. The faster quantum decoherence of the 3AMP spacer-based perovskite underpins a longer nonradiative lifetime, offering insight into its superior performance as an optoelectronic material. This work sheds light on the relationship between structural fluctuations and charge carrier dynamics that can help in designing 2D perovskites with superior photophysical properties.

Here, we perform an in-depth comparison of structural and electronic properties that dominate the photophysical properties of two different DJ phase perovskites, (3AMP) (MA) 3 Pb 4 I 13 and (4AMP)(MA) 3 Pb 4 I 13 . To consider dynamical processes influencing charge carriers at ambient conditions, we employ ab initio real-time time-dependent density functional theory (DFT) and nonadiabatic molecular dynamics (NAMD) modeling. We find that different spacer cations significantly impact the thermal structural variations in these halide perovskites. (3AMP) (MA) 3 Pb 4 I 13 has larger structural fluctuations causing faster quantum dephasing compared to (4AMP) (MA) 3 Pb 4 I 13 . Comparable nonadiabatic coupling with faster quantum dephasing leads to a long carrier lifetime in 3-AMP based perovskites compared to 4-AMP based perovskites. In addition, electron-phonon coupling in DJ phases may have a decisive impact on the reduction of extrinsic non-radiative recombination and carrier trapping as shown previously for RP phases [START_REF] Gong | Electron-phonon interaction in efficient perovskite blue emitters[END_REF] .

We start our simulations with previously reported structures [START_REF] Mao | Hybrid Dion-Jacobson 2D lead iodide perovskites[END_REF] as the initial geometries for (3AMP)(MA) 3 Pb 4 I 13 and (4AMP)(MA) 3 Pb 4 I 13 . In these geometries, every slab with four layers of [PbI 6 ] 4-octahedra remain separated by a layer of 3AMP or 4AMP dications. We further refer to these materials as 3AMP and 4AMP perovskites, respectively. Due to the presence of light atoms and thermally activated molecular dynamics, X-ray diffraction-based structural determination techniques cannot discern the exact conformation of the asymmetric organic cations in between the two inorganic layers with absolute certainty, as well as structural distortions related to inversion symmetry breaking. As a result, recent studies have reported two different interlayer cation alignments in the spacer region [START_REF] Park | Ferroelectricity and Rashba effect in a two-dimensional Dion-Jacobson hybrid organic-inorganic perovskite[END_REF][START_REF] Kang | Dynamic Disorder and Potential Fluctuation in Two-Dimensional Perovskite[END_REF] . Thus, we initially consider both spacer-cation conformations (Fig. 1 (c) -(f) and Supplemental Fig. S1) because these two conformations may coexist. In these two conformations, the nearest-neighbor spacer cations (that is 3AMP or 4AMP molecules) along the b-axis remain aligned (Fig. 1(c), (d)) or anti-aligned (Fig. 1(e), (f)) to each other in the ab-plane. Our structure optimizations (see SI Computational Methods) with fixed cell parameters and without any constraints on atomic positions reveal that for both perovskites, the aligned and anti-aligned structures are energetically very similar as tabulated in Supplemental Table S1. The experimental bandgap values [START_REF] Mao | Hybrid Dion-Jacobson 2D lead iodide perovskites[END_REF] of these two perovskites are 1.87 eV and 1.89 eV respectively. Simulated band gap values without considering spin-orbit coupling (SOC) of aligned structures and anti-aligned structures are shown in Table 1. The band gap values without SOC correction are 10's of meV away from the experimentally measured data for the 4AMP system. This seemingly excellent agreement is actually due to a well-documented mutual cancelation of errors that stems from leaving out spin-orbit coupling [START_REF] Even | Importance of spin-orbit coupling in hybrid organic/inorganic perovskites for photovoltaic applications[END_REF][START_REF] Pedesseau | Electronic properties of 2D and 3D hybrid organic/inorganic perovskites for optoelectronic and photovoltaic applications[END_REF][START_REF] Even | DFT and k• p modelling of the phase transitions of lead and tin halide perovskites for photovoltaic cells[END_REF][START_REF] Meggiolaro | First-principles modeling of defects in lead halide perovskites: best practices and open issues[END_REF][START_REF] Ghosh | Mixed A-cation perovskites for solar cells: atomic-scale insights into structural distortion, hydrogen bonding, and electronic properties[END_REF][START_REF] Ghosh | Putting the squeeze on lead iodide perovskites: pressure-induced effects to tune their structural and optoelectronic behavior[END_REF] (band gap narrowing) and electronic many-body interactions [START_REF] Onida | Electronic excitations: density-functional versus many-body Green'sfunction approaches[END_REF] (band gap widening). Unfortunately, with 260 atoms per simulation cell, the inclusion of SOC is prohibitive for many of the calculations in this work (see SI Computational Methods). Nonetheless, with or without the inclusion of SOC correction, we notice two important trends from our simulations. First, the calculated band gap is largely insensitive to the cation orientation in the case of the 4AMP perovskite, however, the band gap of 3AMP significantly changes with the alignment of its cations. It is known experimentally that the shape of the 3AMP cation by comparison to its 4AMP counterpart is at the origin of strain accumulation in the corresponding DJ compounds or intermediate n values [START_REF] Kepenekian | Concept of lattice mismatch and emergence of surface states in two-dimensional hybrid perovskite quantum wells[END_REF][START_REF] Blancon | Semiconductor physics of organicinorganic 2D halide perovskites[END_REF] . Second, computational findings for the anti-aligned conformation of 3AMP matches closely with the available experimental structural data. To analyze structural features, we rely on previous extensive work relating the backbone structure to the electronic properties of 2D perovskites [START_REF] Pedesseau | Advances and Promises of Layered Halide Hybrid Perovskite Semiconductors[END_REF][START_REF] Katan | Quantum and dielectric confinement effects in lower-dimensional hybrid perovskite semiconductors[END_REF][START_REF] Blancon | Semiconductor physics of organicinorganic 2D halide perovskites[END_REF][START_REF] Fu | Cation Engineering in Two-Dimensional Ruddlesden-Popper Lead Iodide Perovskites with Mixed Large A-Site Cations in the Cages[END_REF] . First, we note that the average Pb-I-Pb bond angle and Pb-I bond length (Supplemental Table S2) for both 3AMP and 4AMP perovskites do not change significantly between the aligned and anti-aligned structures. In order to quantify how much of the band gap opening observed for 3AMP (Table 1) is cation orientation vs perovskite backbone, we place the anti-aligned cations in the perovskite backbone that was optimized for the aligned 3AMP and optimized only the interlayer organic atoms while keeping the perovskite atoms frozen. The band gap from this test is 1.84 eV, indicating that the cation orientation is almost entirely responsible for the band gap opening that matches experiment. Finally, inspection of non-covalent hydrogen-bonding (H-bonding) shows that the average H-bond length between inorganic Pb4I13 layers and spacer cations is shorter in the anti-aligned conformations (Supplemental Table S3) compared to the aligned ones. This further indicates that the anti-aligned structures have stronger H-bonding interactions that stabilize this conformation. detail. First, we examine the dipole moments for the two molecules (see SI Computational Methods). As shown in Supplemental Table S4, when the interlayer cations are aligned, the dipole moments of two nearest-neighbor cations for 3AMP and 4AMP are 9.09 and 6.06 Debye, respectively. In these aligned perovskites, the dominant part of the dipole moments of interlayer cations lies along the out-of-plane direction (c axis), pointing to the aminomethyl. This aminomethyl is more positively charged than (NH 2 ) + within the piperidinium based on the dipole moment direction. When inspecting the partial charge density of the band edge states for the aligned systems we do observe noticeable asymmetries for 3AMP and 4AMP in Supplemental Figure S2. 3AMP are very asymmetric within the perovskite layers with the electron and hole wavefunctions being separated. This separation rationalizes the erroneously calculated low band gap, owing to a known problem of DFT when dealing with charge-separated states [START_REF] Perdew | Physical content of the exact Kohn-Sham orbital energies: band gaps and derivative discontinuities[END_REF][START_REF] Perdew | Density functional theory and the band gap problem[END_REF] .

When the interlayer cations are anti-aligned, the dipole moment of two nearest-neighbor cations for 3AMP and 4AMP is 1.01 and 4.15 Debye, respectively. For anti-aligned structures, dipoles are directed within the ab-plane, which leads to a more symmetric partial charge distribution (Supplemental Figure S3). This large difference in the dipole moment between the aligned and anti-aligned two nearest-neighbor 3AMP rationalizes the band gap mismatch between the aligned conformation and the experimentally measured band gap in this system. These results demonstrate that to accurately model these asymmetric spacer cation-based systems, it is very important to check their possible mutual orientations, corresponding energetics, and electronic properties. As such for subsequent calculations, we only consider the anti-aligned conformations of both DJ phase perovskites. While we show that for the 3AMP system the cations need to be in an antialigned configuration due to limitations of DFT, the actual orientation of the cations in the experiment remains an open question.

Experimentally, 3AMP and 4AMP bandgaps (around 1.9 eV) are in the visible light range. Calculated band structures are shown in Supplemental Figure S4. I-5p orbitals dominate the valence band maximum (VBM) and Pb-6p orbitals dominate the conduction band minimum (CBM) for both perovskites as illustrated in Supplemental Figures S3 andS5. There is still a slight orbital asymmetry for the anti-aligned structures in the charge density distributions that is most likely arises from the residual dipole moments of these two systems. This is likely due to variations in the octahedral tilting angles along the stacking axis and preventing perfect inversion symmetry [START_REF] Mao | Seven-layered 2D hybrid lead iodide perovskites[END_REF] .

To get some insight into the charge carrier transport, we calculated the effective masses from the band-edge structures computed including SOC. As shown in Supplemental Table S5, the effective masses of electrons and holes along with the stacking direction range from 0.22 m0 to 0.63 m0, which is greater than the effective masses of in-plane carriers (0.17-0.39 m0). In Supplemental Figure S2 (a, b), bands near CBM and VBM look almost flat due to the short distance between Γ and B in the reciprocal space. However, if we zoom in to the part of Γ→B (see Supplemental Figure S6), we can observe clear curvatures, suggesting that these DJ-phase perovskites have better charge transport along the stacking direction compared to other 2D perovskites [START_REF] Cheng | Layer-edge device of two-dimensional hybrid perovskites[END_REF] , which typically have very large effective masses and almost flat dispersion curves compared to their 3D counterparts along the stacking axis [START_REF] Katan | Quantum and dielectric confinement effects in lower-dimensional hybrid perovskite semiconductors[END_REF] . When we include SOC corrections, two bands near the VBM of 3AMP remain nearly doubly degenerate, Figure S6 (a). The substantially lighter effective hole of one of the bands near the VBM of this perovskite can improve the charge mobility along the stacking direction significantly compared with that of 4AMP.

Once 3AMP and 4AMP are heated up, such that the following ab initio molecular dynamics This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0066087

(AIMD) trajectory remains stable at 300±5 K for 5 ps , we observe several important structural features. From the T=0K optimized geometries, we find that the averaged [PbI 6 ] 4-octahedra of 4AMP become more distorted compared to that of 3AMP. Further, 3AMP has smaller timeaveraged Pb-I-Pb bond angle deviations from 180° over the simulation cell (Table S6). With a temperature increase, there is less octahedral distortion for both materials, however, 4AMP remains the more distorted as seen from the Pb-I-Pb bond angle, Supplemental Table S6. The asymmetry in the VBM and CBM charge densities that were observed at equilibrium balances out at RT, and charge densities become evenly distributed across the first and third sublayers as shown in Figure 2. Another point to note is the modification of the Jahn-Teller distortions [START_REF] Park | Excited-state vibrational dynamics toward the polaron in methylammonium lead iodide perovskite[END_REF][START_REF] Neukirch | Geometry distortion and small polaron binding energy changes with ionic substitution in halide perovskites[END_REF] in 3AMP.

At equilibrium, the average axial Pb-I bond length is shorter than the equatorial Pb-I bond length for both 3AMP and 4AMP materials. At RT, 4AMP remains the same, however, in 3AMP the average equatorial Pb-I bond length becomes shorter than the average axial bond lengths, Supplemental Table S7. This is one indication that while the Pb-I octahedra are more disordered in 4AMP, 3AMP undergoes more obvious structural changes at ambient conditions. In Supplemental Table S2, 3AMP at equilibrium has a shorter average N-H•••I bond length than 4AMP, which might be the origin of dynamically enhanced organic-inorganic sublattice coupling.

When calculating the root-mean-square deviation of distances between atoms within the simulation cells, there is much more variation along the MD trajectory for 3AMP, Supplemental We next perform NAMD simulations to investigate the nonradiative recombination of charge carriers near the band edges of DJ-phase perovskites. Here, we only consider the intrinsic nonradiative processes (extrinsic processes, such as doping/defects induced carrier trapping, are not included). Since our simulations are performed around RT, we include all energy states within 0.025 eV from the band edges in our excited state dynamics algorithm. This includes states ranging from VBM-1 to CBM+3. We first estimate the radiative lifetime by applying the Einstein coefficient for spontaneous emission, 𝐴 21 (see SI Computational Methods for details). At this level of theory, we find that both materials have a radiative lifetime on the order of μs (in agreement with experiments [START_REF] Dequilettes | Photoluminescence lifetimes exceeding 8 μs and quantum yields exceeding 30% in hybrid perovskite thin films by ligand passivation[END_REF][START_REF] Yang | Superior Carrier Lifetimes Exceeding 6 µs in Polycrystalline Halide Perovskites[END_REF] ), and nonradiative lifetimes on the order 10 ns as summarized in Table 2. Thus, we mostly focus on the nonradiative recombination processes in these materials. Note that, multiexciton-driven recombinations have not been considered here. Excited-state dynamics necessarily involve the coupling of electronic and lattice degrees of freedom. Both elastic and inelastic carrier scattering matter during the electron-phonon process at ambient conditions. The elastic electron-phonon scattering assists the quantum decoherence, resulting in the collapse of superposition to a singular state. The significance and importance of decoherence in condensed-phase nonradiative recombination have been emphasized extensively 53- 56 . The interaction between the system and environment dampens the nonadiabatic effects and, consequently, controls the transition rates between electronic states [START_REF] Jaeger | Decoherence-induced surface hopping[END_REF] . The decoherence time between electronic states i and j can be obtained as the pure-dephasing time in optical response theory using the second-order cumulant approximation [START_REF] Jaeger | Decoherence-induced surface hopping[END_REF] ,

𝐷 𝑖𝑗 (𝑡) = 𝑒𝑥𝑝 (- 1 ħ 2 ∫ 𝑑𝑡′ 𝑡 0 ∫ 𝑑𝑡 ′′ 𝐶 𝑖𝑗 (𝑡 ′′ ) 𝑡′ 0 ) (1) 𝐶 𝑖𝑗 (𝑡) = ⟨𝛿𝐸 𝑖𝑗 (𝑡 ′ ) 𝛿𝐸 𝑖𝑗 (𝑡 -𝑡′)⟩ 𝑡′ (2) (a) (b)
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0066087

where 𝐶 𝑖𝑗 (𝑡) is the unnormalized autocorrelation function of thermal bandgap fluctuation 𝛿𝐸 𝑖𝑗 between two states i and j along the trajectory. The pure dephasing functions are shown in Figure 4 (a). To calculate the pure dephasing times, we fit these functions with a Gaussian function exp(-0.5( 𝑡 𝜏

) [START_REF] Laboratory | Best Research-Cell Efficiency Chart[END_REF] ). As shown in Table 2, 3AMP has a pure dephasing time of 2.7 fs that is substantially smaller than 4.0 fs calculated for 4AMP, being consistent with the fact that the dephasing times are inversely proportional to the magnitude of the electronic energy gap fluctuations [START_REF] Akimov | Persistent electronic coherence despite rapid loss of electron-nuclear correlation[END_REF] . Usually, the fast decoherence process slows the electronic relaxation. This phenomenon can be well understood from the quantum Zeno effect, where a quantum transition halts if the decoherence time is infinitesimal [START_REF] Itano | Quantum zeno effect[END_REF][START_REF] Prezhdo | Relationship between quantum decoherence times and solvation dynamics in condensed phase chemical systems[END_REF][START_REF] Prezhdo | Quantum anti-zeno acceleration of a chemical reaction[END_REF] . Thus, a shorter dephasing time of 3AMP points to a slower intrinsic nonradiative recombination.

Given the importance of electron-phonon scattering processes, the numerical study of the phonon modes that are coupled to the electronic degrees of freedom provides an additional description of excited-state dynamics. The rapid decay of coherence analyzed above suggests that many phonon modes are involved in this process. The spectral density, which is the Fourier transform of the normalized autocorrelation function of the band gaps, identifies the active phonon modes that are coupled to the associated electronic transition. Figure 4 (b) shows that for both 3AMP and 4AMP, the interband transition couples strongly to many low-frequency phonon modes that appear under 500 𝑐𝑚 -1 . This is because the band edge states consist of orbitals of heavy atoms, Pb and I that give a strong contribution. 3AMP generally exhibits a larger magnitude in spectral densities compared to 4AMP, which indicates a stronger interaction between electron subsystem and phonons in 3AMP. The major peak (53 𝑐𝑚 -1 ) for the electron-phonon coupling in 3AMP can be assigned to bending of Pb-I-Pb bonds [START_REF] Park | Excited-state vibrational dynamics toward the polaron in methylammonium lead iodide perovskite[END_REF][START_REF] Pérez-Osorio | Vibrational properties of the organic-inorganic halide perovskite CH3NH3PbI3 from theory and experiment: factor group analysis, first-principles calculations, and low-temperature infrared spectra[END_REF] . The major peak for 4AMP (13 𝑐𝑚 -1 ) can be assigned to the stretching of Pb-I bonds [START_REF] Park | Excited-state vibrational dynamics toward the polaron in methylammonium lead iodide perovskite[END_REF] . We now focus on the inelastic electron-phonon scattering that reveals the extent of the exchange of energy between the electronic and vibrational subsystems and consequently affects the recombination processes in materials. The definition of non-adiabatic coupling (NAC) [START_REF] Tully | Molecular dynamics with electronic transitions[END_REF] (𝑑 𝑖𝑗 • 𝑅 ̇) between i and j states is as follows,

𝑑 𝑖𝑗 • 𝑅 ̇= 〈𝜑 𝑖 | 𝜕 𝜕𝑅 |𝜑 𝑗 〉 𝑅 ̇= 〈𝜑 𝑖 |∇ 𝑅 𝐻|𝜑 𝑗 〉 𝜀 𝑗 -𝜀 𝑖 𝑅 ̇ (3)
where 𝜑 𝑖/𝑗 are wavefunction of i th /j th state of the system, 𝑅 ̇ are phonon velocities, and 𝜀 𝑗/𝑖 are energies of i th /j th states. Since the intra-band relaxations typically happen on hundreds of femtosecond timescales, it is reasonable to assume that all excited carriers relax to the band edges before recombining across the band gap in these two halide perovskites. Generally, the stronger the NAC, the shorter the nonradiative carrier recombination lifetime is. Thus, we focus specifically on the NAC between the states at the band edge. The root means square (RMS) NAC between the VBM and CBM of 3AMP is 0.21 meV, which is slightly stronger than 0.19 meV calculated for 4AMP. 4AMP has a larger time-averaged bandgap than 3AMP (Table 2) and smaller nuclei position fluctuation (Figure S7). According to equation (3), this is a natural explanation for the smaller NAC value of 4AMP. being greater than 0.1 meV. The area under the curve above 0.1 meV is nearly identical for the two materials in question. The NAC values large enough to induce the transition among states have almost the same occurrence in these two materials, as also shown in Supplemental Figure 4 (c). Thus, NAC does not play a decisive role in the determination of the relative rates of nonradiative processes in these two materials.

Based on the analysis above, we find that 3AMP has a noticeably longer intrinsic nonradiative lifetime when compared with 4AMP, see Figure 4(d). A faster ground state population growth means indeed faster nonradiative recombination [START_REF] Ghosh | Charge carrier dynamics in two-dimensional hybrid perovskites: Dion-Jacobson vs. Ruddlesden-Popper phases[END_REF][START_REF] Ghosh | Lattice Expansion in Hybrid Perovskites: Effect on Optoelectronic Properties and Charge Carrier Dynamics[END_REF] . This difference originates almost exclusively from the difference in the dephasing processes in these halide perovskites. Further, a reduction of an intrinsic nonradiative lifetime may contribute to an overall larger carrier lifetime and Voc of 3AMP as found in previous experimental work [START_REF] Mao | Hybrid Dion-Jacobson 2D lead iodide perovskites[END_REF] .

In summary, our time-dependent ab initio NAMD modeling work shines some light on how intrinsic electron-lattice dynamics may affect solar cell device performance for 3AMP and 4AMP perovskite devices. In solar cells, the efficiency is strongly relevant to the dynamics of photoexcited charge carriers. It is important to get an estimation of the timescale of electron-hole recombination and pathways using theoretical and simulation aspects. The current NAMD method we adopted in this work provides a good option, which is also applicable to other 2D hybrid organic-inorganic perovskites. The larger structural dynamics of 3AMP lead to more obvious state energy fluctuation near the band edges and further result in a faster loss of coherence. Combined with nearly identical band gaps and nonadiabatic coupling constants, the faster decoherence of 3AMP leads to a longer intrinsic nonradiative lifetime, which might be a partial explanation for the experimentally higher PCE and performance in solar cells of 3AMP. The energy transfer between electronic and vibrational subsystems is substantially weakened and inhibits the energy loss in these two systems, which leads to longer nonradiative lifetimes. Both materials have an extremely long (microseconds) radiative lifetime at the chosen level of theory. The intense structural fluctuation of 2D DJ perovskites accelerates the decoherence process, leading to slower intrinsic nonradiative recombination, which is beneficial for optoelectronic performances. Our simulations further relate the band-gap sensitivity to a presence of a specific organic cation. But other practical aspects are related to the choice of the organic cation including the material growth: 4AMP-based DJ pure phase crystals have been synthesized up to n=7 like BA-based RP phases [START_REF] Mao | Seven-layered 2D hybrid lead iodide perovskites[END_REF] . The easy growth of 4AMP DJ phases with large n by comparison to 3AMP is attributed to a reduced accumulation of strain within the 4AMP-DJ lattices at room temperature for intermediate n values [START_REF] Kepenekian | Concept of lattice mismatch and emergence of surface states in two-dimensional hybrid perovskite quantum wells[END_REF][START_REF] Blancon | Semiconductor physics of organicinorganic 2D halide perovskites[END_REF] , related to a very small lattice mismatch between the n=1 compound and the 3D perovskite. On the other hand, longer carrier lifetimes may be a natural explanation for the superior photovoltaic performances of 3AMP based solar cells, which exhibit a larger open-circuit voltage (Voc) than their 4AMP based counterparts. In addition, the larger short-circuit current (Jsc) for 3AMP was attributed to an enhanced carrier transport due to shorter interlayer I-I distances and smaller effective masses [START_REF] Mao | Hybrid Dion-Jacobson 2D lead iodide perovskites[END_REF] . This leads to insight into ways to engineer even better-performing devices. Possible interlayer cation engineering or design could be: finding interlayer cations which lead to shorter X•••X (X=halide) distance [START_REF] Cohen | Dion-Jacobson two-dimensional perovskite solar cells based on benzene dimethanammonium cation[END_REF] . Or shortening the I-I distance by replacing the singleband ring of 3AMP/4AMP with a benzene ring [START_REF] Li | Two-dimensional Dion-Jacobson hybrid lead iodide perovskites with aromatic diammonium cations[END_REF] . Thus, observed atomic-scale insights on the relationship between the dynamic structure and charge carrier recombination processes may help the community to screen a variety of mono-and di-cationic spacer molecules to optimize the optoelectronic properties of 2D-halide perovskites, in addition to constraints related to material This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0066087

growth and heterostructure optimization for device applications. Overall this work provides valuable guiding principles for experimentalists when choosing 2D perovskite materials for optoelectronic applications.

SUPPLEMENTARY MATERIAL

Detailed computational method, structure, partial (band decomposed) charge densities, density of states, band structures of aligned and anti-aligned 3AMP and 4AMP, data for total energies, Pb-I-Pb angles, N-H•••I bond lengths, effective masses for different systems.

Figure 1 .

 1 Figure 1. Molecular structure of (a) 3AMP and (b) 4AMP spacer cations (here, in order to present the structures of these two cations clearly, we show the 3AMP and 4AMP cations without coordination directions). Optimized structures of aligned (c) (3AMP)(MA) 3 Pb 4 I 13 (3AMP) and (d) (4AMP)(MA) 3 Pb 4 I 13 (4AMP) perovskites. (e) and (f) show their anti-aligned conformations, respectively. To emphasize the spacer cations, Pb and I atoms are shown with smaller radii. Yellow, blue, canyon, purple, and pink represent C, N, Pb, I, and H atoms, respectively. The green octahedra represent the Pb-I octahedra.

Figure S7 .

 S7 The larger structural dynamical change leads to larger orbital energy fluctuations along the MD trajectory of 3AMP, Figure3. This difference in energy fluctuations is relevant to the dephasing time, which we discuss below.

Figure 2 .Figure 3 .

 23 Figure 2. Average partial (band decomposed) charge distribution (without SOC) at the Γ point of (a) VBM and (b) CBM of 3AMP at room temperature. (c) and (d) show the same but for 4AMP.

Figure 4 .

 4 Figure 4. (a) Dephasing functions and (b) the spectral density of the transition between VBM and CBM. (c) The zoomed-in histogram of the full range of NAC values. Plotted NACs have been calculated between VBM and CBM of 3AMP and 4AMP. (d) Ground state populations in NAMD relaxation of 3AMP and 4AMP.

Figure 4 (

 4 c) shows a plotted histogram of the frequency of occurrence NAC values of various strengths, where a high NAC value is empirically defined as This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0066087

Table 1 .

 1 Calculated band gap values (without SOC) of aligned/anti-aligned 3AMP and 4AMP perovskites. BG refers to Bandgap, which is the same for the following tables.

	species BG (eV)	aligned-3AMP aligned-4AMP anti-aligned-3AMP 1.42 1.86 1.83	anti-aligned-4AMP 1.89
	Since the cation alignment orientation is of interest to the community, and somewhat surprisingly important to obtain consistent basic electronic structure for 3AMP, we explore this factor in more

Table 2 .

 2 The

	species 3AMP 4AMP	BG (eV, 300K) NAC (meV, rms) Dph time (fs) 1.40 0.21 2.69 1.52 0.19 4.03	RLT (μs) 1.03 1.14	NRLT (ns) 8.20 6.51

calculated time-averaged band gaps along the AIMD trajectories, root means square of non-adiabatic coupling constants (NAC) over simulation time, dephasing time (Dph), radiative lifetime (RLT), and nonradiative lifetime (NRLT) for anti-aligned 3AMP and 4AMP.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0066087

AUTHOR INFORMATION

Corresponding Authors

*E-mail: dibyajnc@gmail.com andajneukirch@lanl.gov

Author Contributions

All authors have approved the final version of the manuscript.