Alessandro Genoni
email: alessandro.genoni@univ-lorraine.fr

On the termination of the X-ray constrained wavefunction procedure: reformulation of the method for a univocal determination of 𝝀 Authors

Keywords: X-ray constrained/restrained wavefunction, Lagrange multipliers, constraints, restraints, quantum crystallography

The long-standing problem of determining the correct point at which stopping the calculations in the X-ray constrained wavefunction (XCW) approach is formally solved on solid mathematical grounds and a new and more physically meaningful variant of the XCW method is proposed.

Introduction and theoretical background

The X-ray constrained wavefunction (XCW) approach (Jayatilaka, 1998;Jayatilaka & Grimwood, 2001;Grimwood & Jayatilaka, 2001;Bytheway, Grimwood & Jayatilaka, 2002;Bytheway, Grimwood, Figgis et al., 2002;[START_REF] Grimwood | [END_REF] is undoubtedly one of the cornerstone methods of modern quantum crystallography (Genoni et al., 2018;Grabowsky, Genoni & Bürgi, 2017;Massa & Matta, 2017;Genoni & Macchi, 2020;[START_REF] Grabowsky | 21st Century Challenges in Chemical Crystallography II -Structural Correlations and Data Interpretation, Structure and Bonding[END_REF]Macchi, 2020). It was originally proposed by Jayatilaka in 1998 and it can be considered as a computational technique having the goal of determining electronic wavefunctions compatible with sets of experimental X-ray diffraction data. To accomplish this task, the XCW strategy aims at looking for wavefunctions that minimise the energy of the investigated system and that, at the same time, reproduce sets of experimental structure factor amplitudes within a desired level of agreement. This is equivalent to minimising the energy of the system under the constraint that the squared goodness-of-fit (GoF 2 ) is equal to a desired value Δ. 1 Namely, the imposed constraint is

GoF 2 = Δ (1)
with GoF 2 given by the following expression:

GoF 2 = 1 𝑁 " -𝑁 # . /𝜂 1𝐹 𝒉 %&'% 1 -1𝐹 𝒉 ()# 13 * 𝜎 𝒉 * 𝒉 (2)
where 𝑁 " is the number of considered experimental structure factor amplitudes, 𝑁 # the number of adjustable parameters, h the triad of Miller indices labelling the reflection, 1𝐹 𝒉 %&'% 1 and 1𝐹 𝒉 ()# 1 the calculated and experimental structure factor amplitudes, respectively, 𝜂 an external h-independent scale factor that is used to put the calculated structure factor amplitudes on the same scale of the experimental ones, and 𝜎 𝒉 the experimental uncertainty associated with 1𝐹 𝒉 ()# 1 (provided that the statistics obeys the normal (Gaussian) law and as long as the estimated error bars are close enough to the true standard deviation of the law).

According to the original formulation of the technique, it is suggested to set Δ in equation ( 1) equal to 1.0 (Jayatilaka, 1998;Jayatilaka & Grimwood, 2001). This means that the goal of the technique is to obtain "constrained wavefunctions" that provide structure factor amplitudes that are, on average, within one standard deviation of the experimental values. Therefore, to minimise the total energy of the system in presence of constraint (1), the following functional was introduced:

𝐽 = 𝐸 +, + 𝜆 9GoF 2 -Δ: (3)
1 Note that throughout the text it was decided to adopt the new terminology recently introduced by Davidson and coworkers (Davidson, Grabowsky & Jayatilaka, 2022a), who correctly recommended to use GoF 2 instead of 𝜒 ! .

with 𝐸 +, as the quantum mechanical energy of the investigated system and 𝜆 as a Lagrange multiplier.

By inspecting equation (3), it is also worth noting that 𝜆 is not dimensionless, but it has the physical dimensions of an energy.

Although multi-determinant variants of the XCW method have been also proposed in recent years (Genoni, 2017;Casati et al., 2017;Genoni, Franchini et al., 2018;Genoni, Macetti, Franchini et al., 2019), the original and most used versions of the Jayatilaka approach assume that the wavefunction to be determined is a single Slater determinant (Jayatilaka, 1998;Jayatilaka & Grimwood, 2001;[START_REF] Grimwood | [END_REF]Hudák et al., 2010;[START_REF] Jayatilaka | Modern Charge-Density Analysis[END_REF]Genoni, 2013a;Genoni, 2013b;Dos Santos et al., 2014;Genoni & Meyer, 2016;Bučinský et al., 2016). This will also be the choice that will be considered throughout this work. Under this hypothesis, obtaining a wavefunction that minimises the energy of the system under exam and that simultaneously reproduces experimental Xray diffraction data within a desired limit is equivalent to making functional 𝐽 (see equation (3) above)

stationary with respect to the molecular orbitals of the single Slater determinant and with respect to the Lagrange multiplier 𝜆.

In the case of a 2N-electron closed-shell system, making 𝐽 stationary with respect to the molecular orbitals (under the additional condition that the molecular orbitals are also orthonormal) leads to the following modified Hartree-Fock equations:

<- - * ∇ * + 𝑣(𝒓) + 𝑣 . (𝒓) + 𝐾 A + 𝜆 𝑣 /01 (𝒓)B 𝜙 2 (𝒓) = 𝜖 2 𝜙 2 (𝒓) (4)
with 𝑣(𝒓) = ∑ being the electron density of the examined system), and with 𝐾 A as the exchange integral operator defined as follows:

𝐾 A 𝜙 2 (𝒓) = - 1 2 J 𝑑𝒓 < 𝛾(𝒓, 𝒓 < ) |𝒓 -𝒓 < | 𝜙 2 (𝒓 < ) (5)
where

𝛾(𝒓, 𝒓 < ) = 2 . 𝜙 = * ? =@- (𝒓 < ) 𝜙 = (𝒓) (6)
Furthermore, again in equation (4), 𝑣 /01 (𝒓) can be expressed as

𝑣 /01 (𝒓) = . 𝐾 𝒉 <Re9𝐹 𝒉 %&'% : cos[2𝜋 (𝑸 A 𝒓 + 𝒒 A ) • (𝑩𝒉)] 𝒉 + Im9𝐹 𝒉 %&'% : sin[2𝜋 (𝑸 A 𝒓 + 𝒒 A ) • (𝑩𝒉)]B (7)
which derives from the introduction of the so-called one-electron scattering operator:

𝐼 b 𝒉 = 𝐼 b 𝒉,C + 𝑖 𝐼 b 𝒉,0 = . exp[𝑖2𝜋 (𝑸 A 𝒓 + 𝒒 A ) • (𝑩𝒉)] (8)
? #

A@with and 𝐼 b 𝒉,C and 𝐼 b 𝒉,0 as its real and imaginary parts (which are hermitian, while 𝐼 b 𝒉 is not), 𝑁 D as the number of symmetry-unique positions in the unit-cell, {𝑸 A , 𝒒 A } as the roto-translations associated with the unit-cell symmetry operations (being 𝑸 A and 𝒒 A a rotation matrix and a translation vector, respectively), and 𝑩 as the reciprocal lattice matrix. Finally, 𝐾 𝒉 in equation ( 7) is given by the following expression:

𝐾 𝒉 = 2𝜂 𝑁 " -𝑁 # 𝜂 1𝐹 𝒉 %&'% 1 -1𝐹 𝒉 ()# 1 𝜎 𝒉 * 1𝐹 𝒉 %&'% 1 (9).
For the sake of completeness, it is also worthwhile to point out that, while 𝜆 has the physical dimensions of an energy (see above), the potential 𝑣 /01 (𝒓) is dimensionless.

In the current implementation of the XCW method, equation ( 4) is solved self-consistently for different values of 𝜆 in a sort of "trial & error" procedure. In other words, the molecular orbitals resulting from the calculation with 𝜆 (2) at the i-th step are used as guess for the computation at the following step with

𝜆 (2E-) = 𝜆 (2) + Δ𝜆.
The procedure is usually iterated, but without a clear and definitive criterion to stop.

Different suggestions have been proposed over the years: from the simple and arbitrary idea of halting at the 𝜆 value for which the SCF process stops converging to a very recent and sophisticated approach introduced by Davidson and coworkers (Davidson, Grabowsky & Jayatilaka, 2022b) who adapted an asymptotic extrapolation technique originally devised by Tozer, Ingamells and Handy [START_REF] Tozer | [END_REF] in methods aiming at obtaining Kohn-Sham-like wavefunctions from highlevel ab initio electron densities given as sets of grid points in real space (Zhao & Parr, 1993;Zhao, Morrison & Parr, 1994;[START_REF] Tozer | [END_REF]. For a more comprehensive overview on the different attempts of establishing a halting criterion, the reader can refer to a recent review on the XCW strategy by Davidson, Grabowsky and Jayatilaka (Davidson, Grabowsky & Jayatilaka, 2022a).

Nevertheless, despite all the remarkable efforts, none of the proposed approaches seems the definitive solution to the halting problem of the XCW method. In fact, although in principle the procedure should be stopped when GoF 2 = 1, this goal is seldom achieved (i.e., GoF 2 > 1) or it may also occur that the computations continue beyond the desired limit, thus leading to final values of GoF 2 that are lower than 1.0. Therefore, in the current XCW implementations, the starting constraint given by equation ( 1) is never really considered and it does not play any role in the working equations. Due to this reason and since the calculated structure factor amplitudes would never exactly match the experimental ones even in the case in which we fully accounted for constraint (1), all the variants of the Jayatilaka approach developed so far should be more correctly renamed as X-ray restrained wavefunction (XRW) methods [START_REF] Jayatilaka | Modern Charge-Density Analysis[END_REF]Grabowsky, Genoni & Bürgi, 2017;Ernst, Genoni & Macchi, 2020;[START_REF] Macetti | [END_REF] or X-ray regularisation procedures (Davidson, Grabowsky & Jayatilaka, 2022a;Davidson, Grabowsky & Jayatilaka, 2022b).

Therefore, based on the previous considerations, it is as the following functional was optimised in the current implementation of the Jayatilaka technique (not the functional given in equation ( 3)):

𝐽 = 𝐸 +, + 𝜆 GoF 2 (10)
In other words, in the current XCW calculations, one tries to determine the wavefunction that minimises the energy of the system and that (on average, not exactly) reproduces as much as possible the given experimental X-ray structure factor amplitudes. Consequently, 𝜆 is not a Lagrange multiplier, but only an external parameter that is manually adjusted during the XCW procedure to modulate the weight of the experimental data in the computations.

The above-described situation stems from the fact that, in the current formulation of the XCW approach, the stationary condition of the original functional 𝐽 (i.e., equation ( 3)) with respect to 𝜆 is completely neglected. In the next section, we will show that if we account for this stationary condition, 𝜆 regains its original meaning of Lagrange multiplier and a clear and univocal way to determine its value simply arises from the derived equations. This will also lead to a slight reformulation of the XCW algorithm.

Determination of 𝝀 and reformulation of the XCW approach

We start from equations (4) with the idea of solving the constrained problem set by relation (3). First of all, by inspecting equations (4), we can observe that they are like usual Hartree-Fock equations except for the additional term 𝜆𝑣 /01 (𝒓), which can be considered as a "constraint potential". Furthermore, 𝑣 /01 (𝒓) is known explicitly (see equation ( 7)) and, for a given set of molecular orbitals, it is predefined.

On the contrary, 𝜆 is known only implicitly and its correct value is the one that should make the wavefunction (and, in our particular case, the molecular orbitals) satisfy the condition expressed by means of equation ( 1). It is also worth noting that each value of 𝜆 is uniquely associated with a set of molecular orbitals. When this set is used to calculate the wavefunction and the corresponding electron density, then functional 𝐽 becomes a function of 𝜆 only. Now, following a reasoning analogous to the one adopted in constrained Density Functional Theory [START_REF] Wu | [END_REF]Wu & Van Voorhis, 2006;Kaduk, Kowalczyk & Van Voorhis, 2012), it is possible to show that 𝐽(𝜆) is a strictly concave function of 𝜆, which has an important repercussion on the way in which 𝜆 can be determined.

Let us consider the first derivative of 𝐽 with respect to 𝜆:

2 𝑑𝐽 𝑑𝜆 = 2 . l 𝛿𝐽 𝛿𝜙 2 * 𝜕𝜙 2 * 𝜕𝜆 + c.c. o ? 2@- + 𝜕𝐽 𝜕𝜆 (11)
However, for a given set of molecular orbitals that are solution of the modified Hartree-Fock equations (4), 𝛿𝐽/𝛿𝜙 2 * = 0 and, consequently, equation ( 11) simply becomes:

𝑑𝐽 𝑑𝜆 = 𝜕𝐽 𝜕𝜆 = GoF 2 -Δ (12)
The meaning of equation ( 12) is that the stationary point of 𝐽(𝜆) (for which 𝑑𝐽 𝑑𝜆 ⁄ = 0) provides the constraint expressed by equation (1). To determine the nature (minimum or maximum) of the stationary point, one needs to evaluate the second derivative of 𝐽(𝜆):

𝑑 * 𝐽 𝑑𝜆 * = 𝜕 𝜕𝜆 /GoF 2 -Δ 3 = 𝜕GoF 2 𝜕𝜆 = . 2𝜂 𝑁 " -𝑁 # 𝜂 1𝐹 𝒉 %&'% 1 -1𝐹 𝒉 ()# 1 𝜎 𝒉 * 𝜕1𝐹 𝒉 %&'% 1 𝜕𝜆 𝒉 (13) Considering that 1𝐹 𝒉 %&'% 1 = r𝐹 𝒉 %&'% /𝐹 𝒉 %&'% 3 * s -/*
, and exploiting the definitions of one-electron scattering operator (see equation ( 8)), we obtain:

𝑑 * 𝐽 𝑑𝜆 * = . 𝐾 𝒉 Re t𝐹 𝒉 %&'% 𝜕/𝐹 𝒉 %&'% 3 * 𝜕𝜆 u 𝒉 (14),
with 𝐾 𝒉 defined by equation ( 9). Afterwards, exploiting again the definition of one-electron scattering

operator and the definition of potential 𝑣 /01 (𝒓) (see equation ( 7)), it is possible to show that 2 Hereinafter, we will extensively use the concept of functional derivative (e.g., 𝛿𝐽/𝛿𝜙 " * in equation ( 11)). The functional derivative can be easily seen as a mathematical quantity that connects the variation in a functional (which is a function having another function as variable) to the variation in the function on which the functional depends.

More mathematically, given a functional 𝐹[𝑓] (namely, a function 𝐹 having another function 𝑓(𝑥) as variable), its differential 𝛿𝐹 is that part of the difference 𝐹[𝑓 + 𝛿𝑓] -𝐹[𝑓] that depends linearly on the variation 𝛿𝑓 of function 𝑓(𝑥); since, in principle, each 𝛿𝑓(𝑥) may contribute to the previous difference, when 𝛿𝑓 is small the differential 𝛿𝐹 can be expressed as follows:

𝛿𝐹 = 1 𝛿𝐹 𝛿𝑓(𝑥) 𝛿𝑓(𝑥) 𝑑𝑥
where 𝛿𝐹 𝛿𝑓(𝑥) ⁄ is the so-called functional derivative of 𝐹[𝑓] with respect to function 𝑓(𝑥) at point 𝑥. For further details, the reader may refer to the appendix dedicated to the concepts of functional and functional derivative in the monography on density functional theory written by Parr and Yang [START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF].

𝑑 * 𝐽 𝑑𝜆 * = 4 . w 𝑑𝒓 𝜙 2 * (𝒓) 𝑣 /01 (𝒓) 𝜕𝜙 2 (𝒓) 𝜕𝜆 ? 2@- (15)
For interested readers, the full details of the derivation of equation ( 15) from ( 13

Consequently, exploiting equation ( 19) in ( 15), the expression of the second derivative 𝑑 * 𝐽 𝑑𝜆 * ⁄ becomes:

3 It is worth noting that equation ( 16) corresponds to the special case in which the variable of functional 𝐹 (i.e., the function 𝑓) depends on a parameter 𝜆, namely 𝐹=𝐹[𝑓(𝑥, 𝜆)]. In this situation, one can write:

𝜕𝐹 𝜕𝜆 = 1 𝛿𝐹 𝛿𝑓(𝑥, 𝜆) 𝜕𝑓(𝑥, 𝜆) 𝜕𝜆 𝑑𝑥
More details can be found again in the monography about density functional theory by Parr and Yang [START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF].

𝑑 * 𝐽 𝑑𝜆 * = 4 . . |∫ 𝑑𝒓 𝜙 2 * (𝒓) 𝑣 /01 (𝒓) 𝜙 & (𝒓)| * 𝜖 2 -𝜖 & &G2 ? 2@- (20),
The sum over 𝑎 can be divided into two contributions (Szabo & Ostlud, 1996): i) one involving a summation over all the virtual orbitals, and ii) another one with a summation over all the occupied orbitals (except 𝜙 2 (𝒓)):

𝑑 * 𝐽 𝑑𝜆 * = 4 . . |∫ 𝑑𝒓 𝜙 2 * (𝒓) 𝑣 /01 (𝒓) 𝜙 & (𝒓)| * 𝜖 2 -𝜖 & H2"I & + 4 . . |∫ 𝑑𝒓 𝜙 2 * (𝒓) 𝑣 /01 (𝒓) 𝜙 = (𝒓)| * 𝜖 2 -𝜖 = ? =@- =G2 ? 2@- ? 2@- (21)
If we consider only the second term and we interchange the indices 𝑖 and 𝑗, we get [START_REF] Szabo | Modern Quantum Chemistry. Introduction to Advances Electronic Structure Theory[END_REF]:

𝑌 = 4 . . |∫ 𝑑𝒓 𝜙 2 * (𝒓) 𝑣 /01 (𝒓) 𝜙 = (𝒓)| * 𝜖 2 -𝜖 = ? =@- =G2 ? 2@- = 4 . . |∫ 𝑑𝒓 𝜙 2 * (𝒓) 𝑣 /01 (𝒓) 𝜙 = (𝒓)| * 𝜖 = -𝜖 2 ? 2@- 2G= ? =@- = -𝑌 (22),
namely the second term on the right-hand side of equation ( 21) is equal to zero. Therefore, the second derivatives 𝑑 * 𝐽 𝑑𝜆 * ⁄ can be simply written as follows:

𝑑 * 𝐽 𝑑𝜆 * = 4 . . |∫ 𝑑𝒓 𝜙 2 * (𝒓) 𝑣 /01 (𝒓) 𝜙 & (𝒓)| * 𝜖 2 -𝜖 & H2"I & ? 2@- ( 23
)
where it is worth stressing again that the sum over i runs over the occupied molecular orbitals, while the sum over a runs over the virtual molecular orbitals. If we assume that, as usual, the occupied molecular orbitals are the lowest eigenfunctions (i.e., 𝜖 2 < 𝜖 & ∀ 𝑖, 𝑎), the second derivative 𝑑 * 𝐽 𝑑𝜆 * ⁄

given by equation ( 23) is always negative. This means that, with respect to 𝜆, we only have a maximum as stationary point.

By optimising 𝐽 with respect to 𝜆, it is thus possible to find the correct value of 𝜆 that provides the electronic wavefunction of the investigated system in presence of the experimental constraint introduced through equation ( 1). This can be done through rapidly converging numerical techniques (e.g., Newton's method) that make use of the easy-to-calculate first and second derivatives considered above. Therefore, following what is currently done in constrained DFT [START_REF] Wu | [END_REF]Wu & Van Voorhis, 2006;Kaduk, Kowalczyk & Van Voorhis, 2012), one could envisage to re-implement the XCW method as an optimisation procedure consisting in an algorithm with two nested loops: i) an

outer loop that will be very similar to the current XCW calculations, with the SCF iterations that are carried out to optimise the molecular orbitals (see equation ( 4)); ii) an inner loop of microiterations (that will be nested within each iteration of the outer loop) with the goal of determining the value of the Lagrange multiplier 𝜆 that makes the wavefunction/electron density satisfy the constraint given by equation ( 1). When convergence is achieved, the process should provide both the desired X-ray constrained wavefunction and the potential that is necessary to fulfil the constraint.

Conclusions and perspectives

In this short communication, we have shown that, by introducing the stationary condition of functional 𝐽 with respect to 𝜆, the long-standing problem of determining the point at which halting the XCW procedure can be elegantly solved. In fact, it has been shown analytically that the right 𝜆-value is a maximum stationary point of 𝐽 with respect to the auxiliary variable 𝜆. In principle, this value could be determined by transforming the current XCW procedure into an algorithm essentially made of two nested loops, an outer one for the optimisation of the molecular orbitals and an inner one for the estimation of 𝜆. Operating in this way, at convergence the desired constraint would be necessarily fulfilled and, unlike the current implementations of the Jayatilaka technique, 𝜆 would regain its original meaning of Lagrange multiplier and we would have a computational strategy more similar to a real Xray constrained wavefunction approach.

Concerning the last point, it is again important to observe that, by introducing a constraint in the form of equation ( 1), even the new version of the method would not allow an exact reproduction of all the considered experimental structure factor amplitudes and, for this reason, it should be also not defined as "X-ray constrained wavefunction technique". However, to distinguish the new variant of the strategy from the old one (which should be better called as "X-ray restrained wavefunction approach"; see above), we propose to introduce the term "weakly X-ray constrained wavefunction (wXCW) method".

In fact, through the new version, at least the final GoF 2 would be indeed constrained to its desired value Δ. If one wanted to reproduce the values of all the structure amplitudes, it would be necessary to introduce one Lagrange multiplier for each experimental observation. This would not modify the mathematical derivation outlined above (actually, as a generalisation, the values of the different Lagrange multipliers would still define a maximum stationary point of functional 𝐽) and would give rise to a fully (or strongly) "X-ray constrained wavefunction technique". It is also very interesting to note that a strongly XCW method would be practically identical to a strategy originally proposed in the 1960s by Weiss, who envisaged a technique based on perturbation theory to obtain a wavefunction corrected by experimental X-ray diffraction intensities [START_REF] Weiss | X-ray determination of electron distributions[END_REF]. In this context, an interesting reevaluation of Weiss' approach and a thorough discussion of its differences as compared to the Jayatilaka method can be found in a pamphlet on quantum crystallography recently published by Macchi [START_REF] Macchi | Quantum Crystallography: Expectations vs Reality[END_REF].

Finally, we want to point out that, in the new variant of the XCW technique, the incertitude shifts from 𝜆 to Δ (i.e., the desired level of agreement), but with a significant advantage in terms of physical meaning for the method. In fact, instead of varying 𝜆, in the new version of the strategy one could imagine considering different Δ values (with lower bound equal to 1.0 in case of experimental data) and performing a wXCW calculation for each of them. The lowest value of Δ for which convergence was achieved would then represent the highest degree of confidence with which the experimental X-ray diffraction data could be reproduced on average, given the wavefunction ansatz and the basis set chosen for the calculations.

Appendix A. Further details of the mathematical derivation

As mentioned in the main text, here we report the full derivation of equation ( 15) from equation ( 13).

Therefore, the starting point is the initial expression of the second derivative 𝑑 * 𝐽 𝑑𝜆 * ⁄ :

𝑑 * 𝐽 𝑑𝜆 * = 𝜕 𝜕𝜆 /GoF 2 -Δ 3 = 𝜕GoF 2 𝜕𝜆 = . 2𝜂 𝑁 " -𝑁 # 𝜂 1𝐹 𝒉 %&'% 1 -1𝐹 𝒉 ()# 1 𝜎 𝒉 * 𝜕1𝐹 𝒉 %&'% 1 𝜕𝜆 𝒉 (A1) Considering that 1𝐹 𝒉 %&'% 1 = r𝐹 𝒉 %&'% /𝐹 𝒉 %&'% 3 * s -/*
, we have:

𝜕1𝐹 𝒉 %&'% 1 𝜕𝜆 = 1 2 1𝐹 𝒉 %&'% 1 t 𝜕𝐹 𝒉 %&'% 𝜕𝜆 /𝐹 𝒉 %&'% 3 * + 𝐹 𝒉 %&'% 𝜕/𝐹 𝒉 %&'% 3 * 𝜕𝜆 u (A2)
Plugging equation (S2) into equation (S1), we obtain:

𝑑 * 𝐽 𝑑𝜆 * = . 𝜂 𝑁 " -𝑁 # 𝜂 1𝐹 𝒉 %&'% 1 -1𝐹 𝒉 ()# 1 𝜎 𝒉 * 1𝐹 𝒉 %&'% 1 t 𝜕𝐹 𝒉 %&'% 𝜕𝜆 /𝐹 𝒉 %&'% 3 * + 𝐹 𝒉 %&'% 𝜕/𝐹 𝒉 %&'% 3 * 𝜕𝜆 u (A3)
𝒉 Now, let us consider 𝐹 𝒉 %&'% and /𝐹 𝒉 %&'% 3 *

. By exploiting the definition of one-electron scattering operator (see equation ( 8) in the main text), we have:

𝐹 𝒉 %&'% = 2 .~𝜙 2 1𝐼 b 𝒉 1𝜙 2 • = 2 .~𝜙 2 1𝐼 b 𝒉,C 1𝜙 2 • + 𝑖 ~𝜙2 1𝐼 b 𝒉,0 1𝜙 2 • ? 2@- ? 2@- (A4) /𝐹 𝒉 %&'% 3 * = 2 .~𝜙 2 1𝐼 b 𝒉 1𝜙 2 • * ? 2@- = 2 .~𝐼 b 𝒉 𝜙 2 1𝜙 2 • ? 2@- = 2 .~𝜙 2 1𝐼 b 𝒉 J 1𝜙 2 • = ? 2@- 2 .~𝜙 2 1𝐼 b 𝒉,C 1𝜙 2 • -𝑖 ~𝜙2 1𝐼 b 𝒉,0 1𝜙 2 • (A5) ? 2@-
Starting from (A4) and (A5) we can determine

KL 𝒉 %&'% KM and KNL 𝒉 %&'% O * KM : 𝜕𝐹 𝒉 %&'% 𝜕𝜆 = 2 . € 𝜕𝜙 2 𝜕𝜆 •𝐼 b 𝒉 •𝜙 2 ' + €𝜙 2 •𝐼 b 𝒉 • 𝜕𝜙 2 𝜕𝜆 ' ? 2@- = 2 . € 𝜕𝜙 2 𝜕𝜆 •𝐼 b 𝒉 •𝜙 2 ' + € 𝜕𝜙 2 𝜕𝜆 •𝐼 b 𝒉 J •𝜙 2 ' * ? 2@- = 2 . € 𝜕𝜙 2 𝜕𝜆 •𝐼 b 𝒉,C •𝜙 2 ' + 𝑖 € 𝜕𝜙 2 𝜕𝜆 •𝐼 b 𝒉,0 •𝜙 2 ' + € 𝜕𝜙 2 𝜕𝜆 •𝐼 b 𝒉,C •𝜙 2 ' + 𝑖 € 𝜕𝜙 2 𝜕𝜆 •𝐼 b 𝒉,0 •𝜙 2 ' ? 2@- = 4 . € 𝜕𝜙 2 𝜕𝜆 •𝐼 b 𝒉 •𝜙 2 ' ? 2@- (A6) 𝜕/𝐹 𝒉 %&'% 3 * 𝜕𝜆 = 2 . € 𝜕𝜙 2 𝜕𝜆 •𝐼 b 𝒉 J •𝜙 2 ' + €𝜙 2 •𝐼 b 𝒉 J • 𝜕𝜙 2 𝜕𝜆 ' ? 2@- = 2 . €𝐼 b 𝒉 J 𝜙 2 • 𝜕𝜙 2 𝜕𝜆 ' * + €𝜙 2 •𝐼 b 𝒉 J • 𝜕𝜙 2 𝜕𝜆 ' = ? 2@- 2 . €𝜙 2 •𝐼 b 𝒉 • 𝜕𝜙 2 𝜕𝜆 ' * + €𝜙 2 •𝐼 b 𝒉 J • 𝜕𝜙 2 𝜕𝜆 ' ? 2@- = 2 . €𝜙 2 •𝐼 b

  potential, 𝑣 . (𝒓) = ∫ 𝑑𝒓′ 9(𝒓 " ) |𝒓6𝒓 " | as the Hartree potential (𝜌(𝒓)

  ) are given in Appendix A. However, at this point, it is more important to observe that 𝜙 2 (𝒓) depends on the actual constraint potential 𝜆𝑣 /01 (𝒓). Therefore, we can write:

	𝜕𝜙 2 (𝒓) 𝜕𝜆	= w 𝑑𝒓′	𝛿𝜙 2 (𝒓) 𝛿[ 𝜆𝑣 /01 (𝒓′)]	𝜕/ 𝜆𝑣 /01 (𝒓′)3 𝜕𝜆	= w 𝑑𝒓′	𝛿𝜙 2 (𝒓) 𝛿[ 𝜆𝑣 /01 (𝒓′)]	𝑣 /01 (𝒓′)	(16)
	In the previous equation, 3 𝛿𝜙 2 (𝒓) 𝛿[ 𝜆𝑣 /01 (𝒓′)] ⁄	is functional derivative of the i-th occupied
	molecular orbital with respect to the constraint potential 𝜆𝑣 /01 (𝒓) and it can be essentially seen as the
	first-order response of the orbital (𝛿𝜙 2 (𝒓) ) caused by a small change in the constraint potential
	Therefore, using relation (17), we have
			𝛿𝜙 2 (𝒓) 𝛿[ 𝜆𝑣 /01 (𝒓′)]	= . 𝜙 & (𝒓) &G2	𝜙 & * (𝒓′) 𝜙 2 (𝒓′) 𝜖 2 -𝜖 &	(18),
	and substituting into (16) we get:
			𝜕𝜙 2 (𝒓) 𝜕𝜆	= . 𝜙 & (𝒓) w 𝑑𝒓′ &G2	𝜙 & * (𝒓′) 𝜙 2 (𝒓′) 𝜖 2 -𝜖 &	𝑣 /01 (𝒓′)

(𝛿[ 𝜆𝑣 /01 (𝒓′)]). Therefore, first-order perturbation theory can be exploited to evaluate this functional derivative, with 𝛿𝜙 2 (𝒓) expanded in terms of the eigenstates of the modified Fock operator of the XCW method:

𝛿𝜙 2 (𝒓) = . 𝜙 & (𝒓) ∫ 𝑑𝒓′′ 𝜙 & * (𝒓′′) 𝛿[ 𝜆𝑣 /01 (𝒓′′)] 𝜙 2 (𝒓′′) 𝜖 2 -𝜖 & &G2 (17),

where 𝜖 2 and 𝜖 & are the orbital energies associated with molecular orbitals 𝜙 2 (𝒓) and 𝜙 & (𝒓) , respectively, obtained through the resolution of equations (

4

).

  which is equation (14) in the main text and where 𝐾 𝒉 is defined by equation (9).Re t𝐹𝒉 %&'% 𝜕/𝐹 𝒉 %&'% 3 * 𝜕𝜆 u = 4 Re t𝐹 𝒉 %&'% . €𝜙 2 •𝐼 b €𝜙 2 • ‡ 𝐾 𝒉 rRe9𝐹 𝒉 %&'% : 𝐼 b 𝒉,C + Im9𝐹 𝒉 %&'% : 𝐼 b 𝒉,0 sHowever, due to the definition of 𝑣 /01 (𝒓) (see equations (7) and (8) in the main text), we have:𝑣 ˆ/01 = . 𝐾 𝒉 rRe9𝐹 𝒉 %&'% : 𝐼 b 𝒉,C + Im9𝐹 𝒉 %&'% : 𝐼 b 𝒉,0 s

						?
							𝒉 J •	𝜕𝜙 2 𝜕𝜆	'	u
						2@-
											?	?
		= 4 Re trRe9𝐹 𝒉 %&'% : + 𝑖 Im9𝐹 𝒉 %&'% :s …. €𝜙 2 •𝐼 b 𝒉,C •	𝜕𝜙 2 𝜕𝜆 '	-𝑖 . €𝜙 2 •𝐼 b 𝒉,0 •	𝜕𝜙 2 𝜕𝜆 ' †u
											2@-	2@-
					?					?
		= 4 Re9𝐹 𝒉 %&'% : . €𝜙 2 •𝐼 b 𝒉,C •	𝜕𝜙 2 𝜕𝜆	'	+ 4 Im9𝐹 𝒉 %&'% : . €𝜙 2 •𝐼 b 𝒉,0 •	𝜕𝜙 2 𝜕𝜆	'
					2@-					2@-
			?
		= 4 . €𝜙 2 •Re9𝐹 𝒉 %&'% : 𝐼 b 𝒉,C + Im9𝐹 𝒉 %&'% : 𝐼 b 𝒉,0 •	𝜕𝜙 2 𝜕𝜆	'	(A10)
			2@-
	By substituting (A10) into (A9), we obtain:
	? 2@-𝑑𝜆 * = 4 . 𝑑 * 𝐽		𝒉							•	𝜕𝜙 2 𝜕𝜆	'	(A11)
		? 2@-	? €𝜙 2 •𝑣 ˆ/01 •	𝜕𝜙 2 𝜕𝜆 ' = 4 . w 𝑑𝒓 𝜙 2 * (𝒓) 𝑣 /01 (𝒓) ? 2@-	𝜕𝜙 2 (𝒓) 𝜕𝜆	(A13)
	𝒉,C • which is exactly equation (15) in the main text. 𝜕𝜙 2 𝜕𝜆 ' -𝑖 €𝜙 2 •𝐼 b 𝒉,0 •	𝜕𝜙 2 𝜕𝜆 ' + €𝜙 2 •𝐼 b 𝒉,C •	𝜕𝜙 2 𝜕𝜆 ' -𝑖 €𝜙 2 •𝐼 b 𝒉,0 •	𝜕𝜙 2 𝜕𝜆 '
			2@-
			?
		= 4 . €𝜙 2 •𝐼 b 𝒉 J •	𝜕𝜙 2 𝜕𝜆 '					(A7)
			2@-
	And it is easy to show that:
	𝜕/𝐹 𝒉 %&'% 3 * 𝜕𝜆	= 4 . €𝜙 2 •𝐼 b 𝒉 J • ? 2@-	𝜕𝜙 2 𝜕𝜆 '	= 4 . €𝐼 b 𝒉 J 𝜕𝜙 2 𝜕𝜆 ? 2@-	•𝜙 2 ' *	= 4 . € 𝜕𝜙 2 𝜕𝜆 ? 2@-	•𝐼 b 𝒉 •𝜙 2 ' *	= l 𝜕𝐹 𝒉 %&'% 𝜕𝜆	* o	(A8)
	Therefore, exploiting equations (A4-A8), equation (A3) becomes
				𝑑 * 𝐽 𝑑𝜆 * = . 𝐾 𝒉 Re t𝐹 𝒉 %&'% 𝜕/𝐹 𝒉 %&'% 3 * 𝒉 𝜕𝜆	u	(A9),

Now, let us analyse the term Re ƒ𝐹 𝒉 %&'% KNL 𝒉 %&'% O * KM " in equation (A9): 𝒉 (A12) and, therefore, relation (A11) becomes: 𝑑 * 𝐽 𝑑𝜆 * = 4 .

Acknowledgements

The author is indebted to some colleagues for useful and enlightening discussions on the concepts of Lagrange multiplier in the context of Jayatilaka's X-ray constrained/restrained wavefunction approach: Marco Giovanni Genoni (Department of Physics "Aldo Pontremoli", University of Milan, Italy), Lorenzo Ugo Ancarani & Arnaud Leclerc (Laboratory of Theoretical Physics and Chemistry, University of Lorraine, France). All of them are gratefully thanked.

data -to preserve data integrity, please ensure where possible that the IUCr Word tools (available from http://journals.iucr.org/services/docxtemplate/) are installed when editing this document.


Mirela Vasile
email: mirela.vasile@geo.unibuc.ro

Alfred Vespremeanu-Stroe

Daniela Pascal

Regis Braucher

Alin Pleșoianu

Răzvan Popescu

Bernd Etzelmüller

Rock walls distribution and Holocene evolution in a mid-latitude mountain range (the Romanian Carpathians

Keywords: rock wall morphometry, lithology, rock slope failures, Romanian Carpathians rock wall morphometry, lithology, rock slope failures, Romanian Carpathians

Rock walls in high mountain areas are the expression of long-term slopes response (10 3 -10 5 years) to tectonics, weathering and denudation and a major source of sediment and hazard. Mountain rock walls (RW) characteristics and evolution at mountain-range scale is rarely discussed in the literature. Using a database of 791 RW mapped in the Romanian Carpathians, we present their distribution and morphometry in respect to lithological class, structural features and topography and relate them to post-Younger Dryas (Holocene) rock slope failure chronology. Morphometric data indicate that metamorphic and igneous RW (linked to a great extent to glacial valleys and cirques headwalls) are usually restricted to the highest sectors of the mountain slopes, are characterized by reduced relative heights and have an asymmetrical distribution, being common on the North-exposed slopes but extremely rare on the South. Statistical analysis results show the high significance of structural and tectonic control on RW distribution in sedimentary units which imposes the predominance of West and North orientations and RW dimensions up to a degree higher than in other lithologies. Based on 38 10 Be surface exposure ages obtained on metric boulders from the Southern and Eastern Carpathians, we hypothesise that metamorphic and igneous RW in the formerly glaciated Carpathian valleys were significantly shaped during Early Holocene (before 9 ka) by rock slope failures events that followed the deglaciation of the highest cirques and the intense RW permafrost degradation. We associate the long-term imprints of frost weathering to the significant North/South RW and rock glaciers distribution asymmetry, also identified in other mid-latitude mountain sites with similar topographic constraints.

Powered by Editorial Manager® and ProduXion Manager® from Aries Systems Corporation

Rock walls in high mountain areas are the expression of long-term slopes response (10 3 -10 5 years) to tectonics, weathering and denudation and a major source of sediment and hazard.

Mountain rock walls (RW) characteristics and evolution at mountain-range scale is rarely discussed in the literature. Using a database of 791 RW mapped in the Romanian Carpathians, we present their distribution and morphometry in respect to lithological class, structural features and topography and relate them to post-Younger Dryas (Holocene) rock slope failure chronology.

Morphometric data indicate that metamorphic and igneous RW (linked to a great extent to glacial valleys and cirques headwalls) are usually restricted to the highest sectors of the mountain slopes, are characterized by reduced relative heights and have an asymmetrical distribution, being common on the North-exposed slopes but extremely rare on the South. Statistical analysis results

show the high significance of structural and tectonic control on RW distribution in sedimentary units which imposes the predominance of West and North orientations and RW dimensions up to a degree higher than in other lithologies. Based on 38 10 Be surface exposure ages obtained on metric boulders from the Southern and Eastern Carpathians, we hypothesise that metamorphic and igneous RW in the formerly glaciated Carpathian valleys were significantly shaped during Early Holocene (before 9 ka) by rock slope failures events that followed the deglaciation of the highest cirques and the intense RW permafrost degradation. We associate the long-term imprints of frost weathering to the significant North/South RW and rock glaciers distribution asymmetry, also identified in other mid-latitude mountain sites with similar topographic constraints.

INTRODUCTION

Mountain RW are landforms highly sensitive for mechanical weathering and erosional processes [START_REF] Hales | Climatic controls on frost cracking and implications for the evolution of bedrock landscapes[END_REF][START_REF] Matsuoka | Frost weathering and rockwall erosion in the southeastern Swiss Alps: Long-term (1994-2006) observations[END_REF][START_REF] Allen | Extremely warm temperatures as a potential cause of recent high mountain rockfall[END_REF][START_REF] Phillips | Rock slope failure in a recently deglaciated permafrost rock wall at Piz Kesch (Eastern Swiss Alps)[END_REF], the rates of which are dictated by the interplay of lithology, climate and the local uplift regime [START_REF] Willett | Orogeny and orography: The effects of erosion on the structure of mountain belts[END_REF][START_REF] Seong | Rates of basin-wide rockwall retreat in the K2 region of the Central Karakoram defined by terrestrial cosmogenic nuclide 10Be[END_REF][START_REF] Bartosch | Topographic evolution of the Eastern Alps: The influence of strike-slip faulting activity[END_REF]. Tectonics and structure significantly influence the extent and morphometry of the exposed rock surfaces in high mountain areas [START_REF] Lifton | Influence of rock strength on the valley morphometry of Big Creek, central Idaho, USA[END_REF][START_REF] Ellis | A global perspective on the topographic response to fault growth[END_REF][START_REF] Sauchyn | Structural control of the morphometry of open rock basins, Kananaskis region, Canadian Rocky Mountains[END_REF]) which determine differential resistance to weathering and erosion. In formerly glaciated European mountain ranges, weathering and denudation rates are reported quantitatively after the Last Glacial period [START_REF] Curry | Lateglacial and Holocene talus slope development and rockwall retreat on Mynydd Du, UK[END_REF][START_REF] Hughes | Geological controls on Pleistocene glaciation and cirque form in Greece[END_REF][START_REF] Messenzehl | Regional-scale controls on the spatial activity of rockfalls (Turtmann Valley, Swiss Alps) -A multivariate modeling approach[END_REF][START_REF] Matthews | Small rock-slope failures conditioned by Holocene permafrost degradation: a new approach and conceptual model based on Schmidt-hammer exposure-age dating, Jotunheimen, southern Norway[END_REF]. The occurrence of numerous rock slope failures (RSF) in response to local deglaciation debuttressing has been documented by absolute age dating [START_REF] Soldati | Landslides and climate change in the Italian Dolomites since the Late glacial[END_REF][START_REF] Prager | Age distribution of fossil landslides in the Tyrol (Austria) and its surrounding areas[END_REF][START_REF] Ballantyne | Enhanced rock-slope failure following ice-sheet deglaciation: Timing and causes[END_REF] with responses varying from immediate to millennial time lags [START_REF] Ballantyne | Enhanced rock-slope failure following ice-sheet deglaciation: Timing and causes[END_REF]. The reconstruction of RSF chronology in Tatra Mountains [START_REF] Pánek | Cosmogenic age constraints on post-LGM catastrophic rock slope failures in the Tatra Mountains (Western Carpathians)[END_REF] shows that lower magnitude events within steep topography occurred in the highest sectors of slopes hundreds of years after glacier retreat and are likely triggered by ice mass disappearance, whereas complex RSF producing at millennial time-scale in lower topography are associated with climate changing to warmer and more humid conditions during the onset of the Holocene and the Sub-Boreal period [START_REF] Soldati | Landslides and climate change in the Italian Dolomites since the Late glacial[END_REF][START_REF] Ivy-Ochs | Surface exposure dating of the Flims landslide, Graubünden, Switzerland[END_REF][START_REF] Hermanns | Rapid rock-slope failures[END_REF].

High-mountain rock slopes are a continuous source area for geomorphic processes that trigger natural hazards like debris flows, rockfalls or rock avalanches [START_REF] Loye | Identification of potential rockfall source areas at a regional scale using a DEMbased geomorphometric analysis[END_REF][START_REF] Corona | Parameterization of rockfall source areas and magnitudes with ecological recorders: When disturbances in trees serve the calibration and validation of simulation runs[END_REF][START_REF] Kromer | Identifying and Monitoring Rockfall Precursors Using Terrestrial Laser Scanning for Improved Rockfall Hazard Management[END_REF], especially when affected by permafrost degradation. RW stability is responsive to climate variables such as changes of permafrost conditions [START_REF] Krautblatter | Why permafrost rocks become unstable: A rock-ice-mechanical model in time and space[END_REF][START_REF] Girard | Environmental controls of frost cracking revealed through in situ acoustic emission measurements in steep bedrock[END_REF]) and global climate change influencing periglacial processes [START_REF] Gruber | Permafrost thaw and destabilization of Alpine rock walls in the hot summer of 2003[END_REF][START_REF] Messenzehl | Regional-scale controls on the spatial activity of rockfalls (Turtmann Valley, Swiss Alps) -A multivariate modeling approach[END_REF][START_REF] Phillips | Rock slope failure in a recently deglaciated permafrost rock wall at Piz Kesch (Eastern Swiss Alps)[END_REF], which has also been documented in warming conditions during the Holocene in the European and Scandinavian Alps [START_REF] Hormes | 10Be exposure ages of a rock avalanche and a late glacial moraine in Alta Valtellina, Italian Alps[END_REF][START_REF] Nagelisen | Post-glacial rock avalanches in the Obersee Valley, Glarner Alps, Switzerland[END_REF][START_REF] Hilger | Permafrost as a first order control on long-term rock-slope deformation in (Sub-)Arctic Norway[END_REF]. This raises the question of RW evolution and subsequent debris production induced by the post-Younger Dryas permafrost retreat, especially in mid-latitude where permafrost in northerly slopes exists at lower altitudes than in southerly slopes [START_REF] Magnin | Determination of warm, sensitive permafrost areas in near-vertical rockwalls and evaluation of distributed models by electrical resistivity tomography[END_REF].

Range-scale morphometric studies in the Romanian Carpathians have been documenting the distribution of glacial cirques by high-accuracy mapping [START_REF] Mîndrescu | Cirque form and development in Romania: Allometry and the buzzsaw hypothesis[END_REF]Mîndrescu et al. 2010) and object oriented image analysis [START_REF] Ardelean | Clasificarea semi-automată a unor forme de relief pentru cartarea geomorfologică[END_REF]. Geochronology studies based on absolute ages resume to deglaciation history, pointing to a Younger Dryas glacial advance at 12.9-12.1 ka only in the highest massifs (Popescu et al., 2017a). During the last decade, the intensive monitoring of the thermal regime of permafrost susceptible sites has shown the restrictive conditions for permafrost preservation in RW and rock glaciers [START_REF] Vespremeanu-Stroe | The wind regime of Romania -Characteristics, trends and North Atlantic oscillation influences[END_REF][START_REF] Ardelean | Quantifying postglacial sediment storage and denudation rates in a small alpine catchment of the Făgăraș Mountains (Romania)[END_REF]Onaca et al., 2017;Popescu et al., 2017b). Recent studies on RW present state, in terms of stability or thermal regime, emphasized the rockfall hazard imposed by seasonal thawing in steep North-exposed RW (Vasile et al., 2014;[START_REF] Onaca | Near surface thermal characteristics of alpine steep rockwalls in the Retezat Mountains[END_REF]Vasile and Vespremeanu-Stroe, 2017).

In this paper we i) present the distribution of RW in relation to lithology, structure and topography for assessing differential RW retreat control factors and ii) provide a first insight on Holocene RSF occurrence as the last major shaping agent of RW in the Romanian Carpathians.

The objectives are achieved by RW mapping and statistical analysis of distribution controlling parameters, and by in-situ 10 Be surface exposure data analysis in glacial cirques and valleys from five different massifs.

STUDY AREA

The Romanian Carpathians stand as a geographical subdivision of the Carpathian Mountain Arc that stretches in Central and Eastern Europe (44° 30 ' -47° 45' N and 21° 30' -27° 10' E). They expand to a length of 900 km and reach the maximum altitude of 2544 m above sea level (a.s.l.). The three main subdivisions (i.e., the Eastern and Southern Carpathiansabbreviated EC and SC further onand the Apuseni Mountains, Fig. 1a,b) show lithological and topographic differences that reflect the complexity of the geological evolution, structural characteristics, and influence of the Pleistocene glaciations.

The Carpathians are part of the Alpine Orogeny and include tectonic units dating prior to the alpine event, in the Palaeozoic and early Mesozoic. The youngest exhumation phases determined by thermochronology age patterns in the central part of the SC are Latest Cretaceous -Middle Eocene. The south-western sector of the SC underwent Oligocene -Miocene exhumation, whereas most of the EC correspond to Early -Middle Miocene phases, except for the SE Carpathians (Curvature Carpathians) which started uplifting in both Miocene and Latest Pliocene -Quaternary exhumation episodes [START_REF] Merten | Thermo-tectonic evolution of a convergent orogen with low topographic build-up: Exhumation and kinematic patterns in the Romanian Carpathians derived from thermochronology[END_REF]. The EC are built on a central Crystalline Unit (correspondent to present Rodna, Maramureș, Rarău and Hășmaș Mts.), Cretaceous Flysch (Ceahlău and Ciucaș Mts, extending towards the SC in Bucegi Mts) and Palaeogene Flysch (Table 1, Fig. 1b). Internal volcanism during the Miocene led to the formation of EC volcanic massifs while the Pliocene -Quaternary comprised both a rapid uplift of 500-1000 m (which led to the formation of the most recent depression areas), and the strongest volcanic activity in the area [START_REF] Săndulescu | Geotectonica României[END_REF][START_REF] Linzer | Kinematic evolution of the Romanian Carpathians[END_REF][START_REF] Mutihac | Gologia României[END_REF].

The SC are comprised of three major Crystalline Units: the Getic Overthrust Nappe (Șureanu, Căpățânii, Lotru, Cindrel, Godeanu Mts), the Supragetic Overthrust Nappe (Făgăraș and Iezer Mts) and the Danubian Nappes (Retezat and Parâng Mts), the latter being formed by granitic and granodioritic batholiths in their central areas and marginal limestone massifs (Fig. 1a, b for location, Table 1). [START_REF] Mîndrescu | Geomorfometria circurilor glaciare din Carpații românești[END_REF] Present neo-tectonic movements show a differential uplift trend of the Carpathian orogeny with mean values of 1-3 mm/yr, higher values up to 3-5 mm/yr reported in the Eastern Făgăraş, Bucegi Mts. and the Curvature Carpathians which are associated with the activity from Vrancea seismic region [START_REF] Hoeven | Observation of present-day tectonic motions in the Southeastern Carpathians: Results of the ISES / CRC-461 GPS measurements[END_REF].

The past glacial activity in the Carpathian area is expressed by well-preserved glacial cirques, valleys and associated rock walls (Mîndrescu et al., 2010), most of which were modelled during

LGM and Late Glacial cold phases (Popescu et al. 2017a).

The [START_REF] Micu | Climate of the Romanian Carpathians[END_REF]. In-situ RW thermal monitoring above 2200 m a.s.l. (Vasile and Vespremeanu-Stroe, 2017) exhibits prolonged seasonal frost (140-150 days/season with potential frost penetration depths reaching 2 m), and mean annual rock temperatures (MART) of 0.5°C on the North-exposed slopes and 3-4°C higher MART on the southern slopes, where daily temperature oscillations prevail, and continuous freezing rarely sets within the -3…-8°C freezing window.

METHODS AND DATA

Rock wall mapping

RW were mapped based on the available time records of Google Earth satellite imagery.

Because some images were not clear enough for a good differentiation between the RW and the adjacent geomorphological units, a comparison with higher resolution air photography was undertaken (orthophoto images available for view only from the National Agency for Cadastre and Land Legislation -ANCPI at 1-5 m resolution). Further, the 25 m resolution EU-DEM digital surface model (EEA) was used to check slopes inclination and the inflection points within longitudinal profiles at the contact with the talus or at the top of the glacial cirques. The term rock wall refers herein to steep, bare and compact rock surfaces, with angles usually higher than 37-40 degrees [START_REF] Gruber | A mass-conserving fast algorithm to parameterize gravitational transport and deposition using digital elevation models[END_REF]. We took into account RW with areas larger than 200 m 2 in order to avoid patchy rock surfaces partially covered with vegetation or sporadic discontinuous outcrops. Considering these constraints, the analysis resumes to 21 mountain units in the EC and the SC (Fig. 1a for location), where rock surfaces matching these criteria were identified on satellite imagery. The geological map of Romania (scale 1: 200 000, Geological Institute of Romania) was used to determine the rock type of each mapped RW, assuring a complete spatial coverage over the entire range.

Rock wall morphometry

Mean RW area, total coverage, altitude, relative height, orientation and slope values were computed for each mountain unit (Table 2) and per lithology types, using the values from all units developed on the same rock class (sedimentary, igneous, metamorphic or volcanic). Each RW was vectorized and the resulting polygons were overlaid on the DEMs and then used for calculating the morphometric parameters in ArcGIS software. Mean RW orientation was determined by averaging raster direction. The RW area was calculated for a 2D projection of the RW polygons on the EU-DEM.

Statistical analysis

The statistical analysis was performed in RStudio (R version 3.4.0), and consisted in three stages. First, a data normality check was performed using the Shapiro-Wilk normality test [START_REF] Shapiro | An Analysis of Variance Test for Normality (Complete Samples)[END_REF], which indicated the non-normality of the data. Then, for each mountain unit, a Kruskal-Wallis one-way analysis of variance test [START_REF] Kruskal | Use of Ranks in One-Criterion Variance Analysis[END_REF] was performed in order to check if there are any statistically significant differences between groups of quantitative parameters, namely the relations between exposures and morphometry, and between lithology and morphometry. The Kruskal-Wallis Test is the non-parametric alternative to ANOVA (one-way analysis of variance), which checks if the analysed groups are subsets from the same population. The test computes the rank variance of the interest variable for the combined groups, and then calculates the H statistic (Equation 1)

𝐻 = 𝑁-1(𝑔𝑛 𝑛 (𝑡 𝑖 -𝑇 𝑗 ) 2 𝑔𝑠𝑛 𝑛 (𝑡 𝑗 -𝑇 𝑖 ) 2 ) , Equation (1)
where nn is the sample size of the group, g is the sum of the group n, snn is the sum of the corresponding group n, ti is the average observed rank sums for the group, tj is the observed rank for a value in the corresponding group, and Ti is the observed total average rank sums.

( [START_REF] Mckight | Kruskal-Wallis Test[END_REF]. The computed H statistic then indicates whether the groups come from the same population by comparing it to a critical value, which for our analysis corresponds to a 95% confidence or a p-value < 0.05. For H values beyond the critical threshold the Kruskal-Wallis Test indicates strong differences between analysed groups. Finally a post-hoc Dunn's

Test [START_REF] Dunn | Multiple Comparisons among Means[END_REF] for multiple comparisons was performed in order to identify which mountain units have significantly different values of each pair of the analysed parameters.

Surface exposure ages

24 boulders (20 from Făgăraș and four from Bucegi Mts) were sampled for cosmogenic 10 Be exposure dating. The samples are part of an extensive study regarding deglaciation and RSF, which counts more than 120 rock surfaces (unpublished data). During sample processing, the abundance of post-Younger Dryas resulting ages raised questions about RSF triggering the detachment of such boulders, as documented in other European mountain ranges (synthesis of studies in [START_REF] Pánek | Cosmogenic age constraints on post-LGM catastrophic rock slope failures in the Tatra Mountains (Western Carpathians)[END_REF]. The samples included in the present study were collected along 7 glacial valleys in Făgăraș Mts and one valley in Bucegi Mts, ranging from 1205 to 2287 m a.s.l., on metric-size boulders from both valley/cirque centre and peripheral. Sample size varied from 2 to 3 cm thick and sampled rock surfaces were vegetation free. Additionally, we accounted other 14 post-Younger Dryas boulder ages documented in the literature in Rodna [START_REF] Gheorghiu | Testing climate synchronicity between Scotland and Romania since the last glacial maximum[END_REF], Parâng [START_REF] Gheorghiu | Deglaciation constraints in the Parâng Mountains, Southern Romania, using surface exposure dating[END_REF] and Retezat massifs [START_REF] Reuther | Late Pleistocene glacial chronology of the Pietrele Valley, Retezat Mountains, Southern Carpathians constrained by 10Be exposure ages and pedological investigations[END_REF][START_REF] Ruszkiczay-Rüdiger | Limited glacial erosion during the last glaciation in mid-latitude cirques (Retezat Mts, Southern Carpathians, Romania)[END_REF] which were considered as outliers in studies regarding deglaciation chronology, raising the database to 38 values from 16 valleys/cirques.

Secondary, we aimed for absolute dating of a rock glacier (Doamnei RG) surface in Făgăraș

Mts., where we sampled four boulders from the RG body on a longitudinal profile but also the source RW area above; for comparing a North/South RW exposure, the corresponding Southface of Doamnei RW was also sampled.

Dating procedure for Făgăraș and Bucegi samples is described below.

The sample purification followed the procedure of [START_REF] Merchel | An update on radiochemical separation techniques for the determination of long-lived radionuclides via accelerator mass spectrometry[END_REF]. Samples were crushed and sieved to the 0.25 -1 mm fraction. Magnetic separation was performed on all samples with a magnetic separator "Frantz LB-1". The other minerals that are embedded in samples were eliminated with mixtures of HCl and H2SiF6. Then atmospheric 10 Be was eliminated by HF (48%) dissolutions. Before the total dissolution, 150 mg of a 9 Be carrier solution (concentration 3025 ± 9 μg/g; [START_REF] Merchel | Towards more precise 10 Be and 36 Cl data from measurements at the 10 À14 level: Influence of sample preparation[END_REF] manufactured in-house from a phenakite crystal were added to the samples. The total dissolution of quartz was performed with HF 48% (3.6 mL per g of quartz and 30 mL in excess). The resulting solutions were evaporated until dryness and samples were recovered with hydrochloric acid. Subsequently samples were precipitated with ammonia before successive separations through an anion exchange column (Dowex 1X8) to remove iron and a cation exchange column (Dowex 50WX8), and to discard boron and recover Be [START_REF] Merchel | An update on radiochemical separation techniques for the determination of long-lived radionuclides via accelerator mass spectrometry[END_REF]. BeO targets were prepared by mixing Niobium powder with the BeO oxide for AMS measurements.

All samples were chemically performed at Laboratoire National des Nucléides Cosmogéniques (LN2C) at CEREGE (Aix en Provence, France) and targets of purified BeO were prepared for AMS measurement at ASTER, the French National AMS Facility (CEREGE, Aix en Provence).

The measurements were calibrated against an In-House standard (STD11) Braucher et al., 2015 standard, using an assigned 10 Be/ 9 Be ratio of (1.191) × 10 -11 (1.09%). Analytical uncertainties (reported as 1ϭ) included for all samples. The 10 Be half-life of (1.387±0.01) × 10 6 years [START_REF] Chmeleff | Determination of the 10Be half-life by multicollector ICP-MS and liquid scintillation counting[END_REF] was used.

Production rates were scaled following Stone, 2000 with a sea level high latitude production rate of 4.02±0.36 atoms/g SiO2/yr [START_REF] Borchers | Geological calibration of spallation production rates in the CRONUS-Earth project[END_REF]. Rock density of 2.5 g/cm 3 was used for all samples. Topographic shielding was calculated using the CosmoCalc 2.

RESULTS

RW distribution

From the 21 units considered in this study, 11 preserve glacial landforms (Table 1). In the EC, only Rodnei, Maramureş, and Călimani massifs present visible glacial landforms (cirques, valleys and moraine deposits), most of which are in Rodnei Mts. [START_REF] Mîndrescu | Geomorfometria circurilor glaciare din Carpații românești[END_REF]. The best preserved landforms are in Făgăraş Mts., where 207 glacial cirques were mapped (Mîndrescu et al. 2010), in Retezat Mts., which show extensive moraine deposits, glacial lakes and complex glacial valleys [START_REF] Urdea | Munții Retezat. Studiu Geomorfologic[END_REF], and in Parâng Mts., which keep the largest glacial cirques in the Romanian Carpathians [START_REF] Iancu | Masivul Parâng. Studiu Geomorfologic[END_REF]. The distribution of glacial cirques in the SC particularly reflects the main ridge orientation, with moderate differences of frequencies between North and South-exposed cirques (17% on North and 11% on South respectively, reported to 45 degrees bins) and the most favourable conditions of cirque glaciers formation on the East-exposed slopes of the valleys (19.5%) due to the strong eastward aeolian snow-transport acting on the crests and plateaus [START_REF] Vespremeanu-Stroe | The wind regime of Romania -Characteristics, trends and North Atlantic oscillation influences[END_REF]. In terms of area and height, the North exposures preserve generally larger and wider cirques [START_REF] Mîndrescu | Geomorfometria circurilor glaciare din Carpații românești[END_REF].

A total of 791 RW were identified and considered as individual features, most of which were mapped in the SC. In most of the mountain units considered here, the main ridges follow East-West or NE-SW direction (Table 1), tracking the principal fault lines (Fig. 1a,b). The distribution of the RW is further presented, based on mean orientation and RW altitude. The units from the EC count 88 mapped RW in total. The RW found in the schist-prevailing massifs are distributed mainly on NE (23%) and secondary on North and East (Fig. 2a) with an average altitude of 1950 m a.s.l., the North-exposed ones being situated at slightly lower altitudes. The andesitic and basaltic rock outcrops mapped in the EC, are largely grouped on the North and NE (82%) similarly with the metamorphic ones.

In the sedimentary units, RW extend on all orientations, with a maximum frequency on the South (almost 40%) whilst just a few (7%) were mapped on the northern slopes (Fig. 2b), but are limited to altitudes lower than 1800 m, reaching an average of 1634 m a.s.l. which is considerably lower (> 300 m) than the metamorphic and volcanic RW.

Compared to the EC, the number of RW mapped in the mountain units from the SC built on metamorphic rocks is much larger, rising to 331 from which 275 are distributed in two large massifs (Făgăraș and Iezer), and the remaining are spread in 6 units characterized by gentler topography and lower altitude (Fig. 3a, Table 2). The northern direction clearly dominates in both clusters, summing 41% from the total RW number, with similar frequencies on the NW and NE bins (Fig. 3d). The second highest RW frequency corresponds to the eastern orientation (104 RW), followed by the western slopes (77 RW). RW mapped on the southern slopes are scarce and represent 5.4% from the total number. In the large metamorphic massifs of Făgăraș and Iezer the highest RW density is in the range of 2100 -2400 m with a mean altitude of 2200 m (Table 2), but 14 mapped surfaces are higher than 2400 m. Parâng and Retezat Mts. in the SC (Fig. 1 for location) are two examples of mixed lithology, being composed mainly of granitoids plus granodiorites intrusions and of crystalline schist, micaschist, amphibolite (around 23% of the mapped surfaces) (Tables 1,2) (Fig. 3b).

Approximately 54% of the RW from Parâng and Retezat Mts. are North-oriented. Eastern and western exposures account for 25% and 15% of the total number, and only 5.7% of the mapped surfaces from the two massifs were identified on the southern slopes. More than half of the RW are concentrated in the 2100 -2300 m interval and almost 18% extend above 2300 m, the northern ones reaching the highest altitudes. The RW on metamorphic rocks in the two massifs range at slightly lower altitudes and occur almost evenly on the East, North and West-oriented slopes. In the limestone-prevailing units from the SC, (Bădescu and Tîrlă, 2020; see Fig. 1 for location), 69 RW were mapped (Fig. 3c), of which almost 40% are exposed towards West directions, 26%

to the East, 23% are on northern slopes, and 11% on the southern ones (Fig. 3d). In terms of altitudinal distribution, there are also major discrepancies between the massifs, imposed by the structural characteristics of each (Table 1). The statistical analysis shows the clear dominance of West-exposed RWs for most of the considered morphometric parameters compared to the other main orientations (Table 3). This asymmetry is however case-specific and imposed by the particular orientation of many of the sedimentary units, respectively the NNE-SSW-oriented ridges and plateaus of Piatra Craiului and Buila hogbacks, Bucegi, Ciucaș and Hășmaș synclines [START_REF] Mutihac | Structura geologică a teritoriului României[END_REF] or Cerna Valley half-graben [START_REF] Povară | Water flow system within the Cerna Valley graben structure (SW of the Southern Carpathians, Romania)[END_REF], which together with the eastward and southward dipping strata contribute to the larger occurrence of West and North-exposed RWs (Fig. 1b). 3: Results of the statistical analysis (Post-Hoc Dunn's Test) of morphometric parameters for pairs of main exposures. Numbers represent the z-Score, which indicates whether the tested parameter pair has a value above the rank mean (positive value), or below (negative value). For example, the mean area of N/E pair has a z-Score = 2.38, meaning that North has a greater mean area than East. In a similar way, for z-Scores below the rank mean the comparison is read inversely, as in the altitude for S/N pair which has a z-Score = -6.05, meaning that North has a higher altitude than South * p-value (0.05, 0.01]; ** p-value (0.01, 0.001]; *** p-value < 0.001. A p-value < 0.05 indicates a strong statistical significance at 95% confidence level In summary, our observations indicate that RW in metamorphic and igneous units are generally restricted to altitudes higher than 2100 m, show a high density on the North-exposed slopes and lack from the southern ones. In comparison, RW distribution in sedimentary units is spread over a larger range of altitudes but highly dependent on the structure and strata dip-direction which result in prevalent West and North exposures but not as large (asymmetric) as for the metamorphic and igneous rocks.

Ratio Mean area

Ratio

RW morphometry

Although the most numerous RW were identified on metamorphic rocks, this group covers only 21% of the total mapped surface. An even smaller cumulated area corresponds to igneous rocks which represent 12% of the total area. Comparatively, although counting less RW, the sedimentary massifs cumulate 67% of the total mapped area due to their large mean surface (Fig. 4a,b). Consequently, igneous and metamorphic RW are 5 to 6 times smaller than limestone, conglomerates and sandstone RW (Fig. 5a). This visible heterogeneity is imposed by Bucegi and Piatra Craiului, with RW of up to 250 m height. This is also consistent with the wider range of mean RW altitude, as limestone surfaces extend at much lower altitudes (Fig. 5c). Summing up, the high values of the z-scores in Table 4 (maximum -9.01 for schist/limestone area, 10.9 for mean altitude) support the high morphometric differences that exist between the sedimentary units and the metamorphic and igneous ones, which are more similar (-0.77 for schist/granite area, -4.3 for altitude). Although more numerous, metamorphic and igneous RW in the study area present typically small-size surfaces perched to the highest stands of formerly glaciated valleys and cirques. The sedimentary (limestone and conglomerates prevailing) RW are wider, steeper and cover larger areas than all the other lithological groups independent of glacial landmarks and with apparent homogeneity in respect to slope orientation. These characteristics thus reflect major differences between control factors over RW morphometry and distribution depending on geology.

Absolute ages

Sample locations, altitudes, 10 Be concentrations and 10 Be surface exposure ages determined in this study and selected from the literature are presented in Table 5. Multiple values (3 to 5 per valley) available in 8 of the valleys presented here, allowed an intra and inter-massif analysis of age distributions. Values range between 0.97 ± 0.08 ka and 11.3 ± 1.0 ka (Table 5), clustering between 11.3 and 9.1 ka (21 values), thus immediately following the Younger Dryas (12.9-11.7 ka, Rasmussen et al., 2006;12.6-11.4 ka, Tǎmaş et al., 2005) and within time lags of up to 2.6 ka after. A second high frequency cluster was found between 9.0 and 7.0 ka, while younger ages were identified with a frequency of 1.75 values/1000 yrs. There is not a clear correlation between the absolute ages of the rock surfaces and the altitude. Such an attempt would be hindered by the unequal sample distribution, given the fact that most of the boulders are situated above 1800 m a.s.l. Nevertheless, most of the youngest ages (0.94 ± 0.08 ka, 4.51 ± 0.18 ka, 2.19 ± 0.12 ka, 1.27 ± 0.22 ka) were found in the lower or mid-sectors of Dejani, Bâlea and Sâmbăta valleys (Făgăraș Mts.), between 1200 and 1500 m a.s.l., respectively. Although few ages from the onset of the Holocene were also determined below 1400 m in the Retezat and Parâng Mts., (10.7 ± 1.4 ka, Ruszkiczay-Rüdiger et al., 2021, 10.5 ± 0.9 ka, [START_REF] Gheorghiu | Deglaciation constraints in the Parâng Mountains, Southern Romania, using surface exposure dating[END_REF], most boulders dating between 11.5 and 7.0 ka are situated in the highest sector of the valleys and on cirque floors (1800-2200 m a.s.l.). This distribution pattern is similar when comparing both valleys within the same massif (e.g., Făgăraș) and valleys from the five different massifs (Table 5). [START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF] and a spallation production rate of 4.49 ± 0.39 atom (g SiO2) -1 a-1 [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10 Be and 26 Al measurements[END_REF] 5. DISCUSSION

Structure and lithology influence on RW distribution

The RW inventory and morphometric analysis results have emphasized the significant influence of the lithology and geological structure on the characteristics of metamorphic, igneous and sedimentary RW in the Romanian Carpathians, the first two rock categories producing much smaller RW, but developed at higher altitudes, than sedimentary massifs which account for the greatest rock surface coverage overall. Our results also showed that the distribution of RW in the Carpathians is very particular in respect to orientation, with an obvious asymmetry between North and South exposures especially for metamorphic and igneous rocks, yet with extent RW on all orientations in limestone or conglomerate massifs.

Despite the North / South balanced distribution of the glacial cirques [START_REF] Mîndrescu | Geomorfometria circurilor glaciare din Carpații românești[END_REF], we showed that RW present a high asymmetry in the metamorphic and igneous mountain units from the both SC and EC (Fig. 2a,4b,d), where South-exposed slopes are much less frequent whilst the total covered area is almost 30 times higher on the North-exposed RW.

Correspondently, the present rock glaciers distribution in these units accounts for 58% of the mapped rock glaciers in the northern quadrant, and only 13% in the southern one (Onaca et al., 2017), which also suggests a more intense/frequent debris accumulation on the North exposed slopes during Younger Dryas and Early Holocene when presumably most of the rock glaciers formed [START_REF] Onaca | Assessement of internal structure of periglacial landforms from southern carpathians (Romania) using dc resistivity tomography[END_REF]. A similar distribution is described in the Adamello-Presanella massif (Italian Alps), where the main ridge follows the NE-SW direction of the North-bordering fault and valleys radiate from the main ridge, covering all the cardinal directions. Here, based on the inventory of 216 rock glaciers mostly consisting of intrusive granodioritic and tonalitic rocks, [START_REF] Baroni | Distribution and behaviour of rock glaciers in the Adamello-Presanella massif (Italian Alps)[END_REF] show that both active and relict rock glaciers are predominantly facing the North, NW and NE compared to the southern quadrant (which counts 18% of the active / inactive, and 15% of the relict ones), and argument, by comparing front altitudes, that local topoclimate makes northern slopes more favourable to rock glacier formation and preservation.

We further consider that RW preservation conditions are also more restrictive on the southern slopes in igneous and metamorphic slopes, as commented below.

Both granite and schist RW in the study area are characterized by metric joints networks which fit well with the dimensions of the boulders enclosed into the adjacent debris deposits (Vasile and Vespremeanu-Stroe, 2017), supporting intense slope modelling which could have led to the formation of large debris deposits, talus cones and rock glaciers. In a simulation of moisture availability in alpine RWs, [START_REF] Rode | Simulation of moisture content in alpine rockwalls during freeze-thaw events[END_REF]) highlight that the preconditions of water saturation and temperature required for ice segregation are often recorded on the North-exposed slopes but just rarely met on warmer South-exposed rock surfaces, which implies that the latter are not prone for large-size debris production. Thus, South-exposed rock slopes would be subject to small-scale flake and granular rock shattering under the effect of both superficial freezing during snow melting intervals [START_REF] Matsuoka | Frost weathering and rockwall erosion in the southeastern Swiss Alps: Long-term (1994-2006) observations[END_REF] and of diurnal insolation thermal stress during snow-free intervals [START_REF] Eppes | Deciphering the role of solar-induced thermal stresses in rock weathering[END_REF]. We assume that South-exposed RW in the Romanian Carpathians were generally more stable than the other exposures whilst the northern were the most active due to longer permafrost preservation. This is also supported by the old age (52.46 ± 0.89 ka) age yielded by the South-exposed ridge outcrop above Doamnei rock glacier, which was apparently unaffected by LGM, when it most probably stood as nunatak. Humification process (i.e. humus formation in soil profiles) has been inferred to be more intense on Southexposed mountain slopes, where warmer conditions intensify oxidation and create a more optimal environment for microorganisms that degrade organic matter [START_REF] Egli | The effects of exposure and climate on the weathering of late Pleistocene and Holocene Alpine soils[END_REF], compared to North-exposed mountain slopes which incorporate undecomposed or weakly degraded organic matter and are subject to mineral leaching due to colder and wetter conditions. [START_REF] Savi | Efficiency of frost-cracking processes through space and time: An example from the eastern Italian Alps[END_REF] reconstruct frost-cracking intensity and debris production during the Holocene in the Eastern Italian Alps, and emphasizes that high debris accumulation occurred during Early Holocene and also during Atlantic and Subatlantic periods when positive MAAT would have promoted continuous superficial (up to 100 cm deep) frost cracking in the highest peaks (around 3000 m a.s.l.). A similar pattern is supported in the SC by the surface exposure ages that sustain production of large debris in all massifs during Early Holocene and secondary debris production in subsequent phases.

The cumulated effect of these processes could explain a generally faster cover with soil and vegetation on the sunny slopes of both metamorphic and igneous units from this study, whereas on the North, colder thermal regime and the production of large-size boulders led simultaneously to a better preservation of RW, which is reflected in present-day distribution and morphometry.

The large synclines represented by Bucegi and Ciucaș Mts., with main North to South dip direction of the conglomerate and sandstones bedding planes uplifted large RW on the Northfacing cuesta fronts, typical for sedimentary units formed as synclines, perched synclines, or hogbacks which are generally dominated by the compactness and steepness of the cuesta escarpments [START_REF] Huggett | Fundamentals of Geomorphology[END_REF]. This is also the case in the NE-SW-dipping limestone massifs in SC which enhanced the development of the largest RW on their western slopes. Limestone RW in both EC and SC generally lack dense superficial joint networks which, along with increased permeability, limits water availability within the first centimetres of rock and implicitly turns ice segregation less probable, which implies reduced RW modelling by frost shattering and debris accumulation in sedimentary massifs [START_REF] Johnson | The effect of topography, latitude, and lithology on rock glacier distribution in the Lemhi Range, central Idaho, USA[END_REF]. This is reflected by the low number of rock glaciers formed / identified in the SC on sedimentary rocks (only 15 on limestone from a total statistical population of 306 rock glaciers; Onaca et al. 2017). Therefore, we consider the lithology and structure to play the major role in imposing the orientation-related homogeneity which accounts as the primary control in RW distribution and dimensions in the sedimentary units from this study. Secondary, post-glacial RW relaxation would have led to the detachment of massive limestone and conglomerate blocks as sustained by the absolute ages obtained in Bucegi Mts (Table 5). However, in specific cases, such as Piatra Craiului limestone hogback, large debris deposits have accumulated at the base of the main tectonic slopes. For such cases, we assume that the absence of transversal valleys and of the Pleistocene glaciers, both caused by topography, could have created the conditions for the long-term (e.g. Middle to Late Quaternary) debris accumulation, the formation of which is still to be deciphered.

Holocene dynamics of RW depicted by rock-slope failures

Following Last Glacial Maximum deglaciation (19-14.5 ka), small glaciers re-occupied the highest cirques (> 2050 m) during Younger Dryas excepting the southern ones [START_REF] Gheorghiu | Deglaciation constraints in the Parâng Mountains, Southern Romania, using surface exposure dating[END_REF]Popescu et al., 2017a;Pascal et al., 2018). The largest Younger Dryas glaciers are likely to have lasted more than a millennium during Early Holocene (e.g. 10.2 ± 0.9 kyrs, [START_REF] Gheorghiu | Deglaciation constraints in the Parâng Mountains, Southern Romania, using surface exposure dating[END_REF]. Therefore, we consider that most of the boulders from high altitudes (> 2000 m) of early Holocene age dated in this study have originated by RW destabilization as response of ice retreat which occurred mainly in the upper valley/cirques sectors where glaciers were restricted (Fig. 6). High frequency of such events occurring several thousand years after Younger Dryas period are well documented in the European Alps [START_REF] Soldati | Landslides and climate change in the Italian Dolomites since the Late glacial[END_REF][START_REF] Cossart | Slope instability in relation to glacial debuttressing in alpine areas (Upper Durance catchment, southeastern France): Evidence from field data and 10Be cosmic ray exposure ages[END_REF][START_REF] Hormes | 10Be exposure ages of a rock avalanche and a late glacial moraine in Alta Valtellina, Italian Alps[END_REF][START_REF] Prager | Age distribution of fossil landslides in the Tyrol (Austria) and its surrounding areas[END_REF][START_REF] Ivy-Ochs | Surface exposure dating of the Flims landslide, Graubünden, Switzerland[END_REF], Tatra Mts. [START_REF] Pánek | Cosmogenic age constraints on post-LGM catastrophic rock slope failures in the Tatra Mountains (Western Carpathians)[END_REF], in Scotland [START_REF] Ballantyne | Enhanced rock-slope failure following ice-sheet deglaciation: Timing and causes[END_REF] and Scandinavia [START_REF] Mercier | The Höfճ ahólar rock avalanche (sturzström): Chronological constraint of paraglacial landsliding on an Icelandic hillslope[END_REF][START_REF] Hilger | Multiple rock-slope failures from Mannen in Romsdal Valley, western Norway, revealed from Quaternary geological mapping and 10Be exposure dating[END_REF][START_REF] Hilger | Permafrost as a first order control on long-term rock-slope deformation in (Sub-)Arctic Norway[END_REF][START_REF] Vick | Evolution and temporal constraints of a multiphase postglacial rock slope failure[END_REF], but also in Karakorum [START_REF] Shroder | The role of mass movements on landscape evolution in the Central Karakoram: Discussion and speculation[END_REF] or the Andes (Fauqué et al., 2009). Many of these sites record re-activations or secondary clusters during the Sub-Boreal period [START_REF] Hermanns | Rapid rock-slope failures[END_REF].

In Fig. 6 we compare frequency curves of the post Younger Dryas boulders dated in the Romanian Carpathians and of the RSF ages compiled from these studies after excluding the mountain ranges which are influenced by excess of humidity/dryness and correspondently by their variability in time (e.g., Himalaya, Atlas or Cascade Mts). Overall, the Romanian Carpathians show a similar general trend with the other world-wide catenae but with an apparently more rapid and accentuated response to the Early Holocene warming and more humid conditions, so that almost 3/4 of the dated RSF occurred before 8 ka with the highest frequency window positioned during 11.6-9 ka. Conversely, the multi-sites curve reflects a higher sensitivity of RW (deduced via RSF occurrence) to the Holocene Climatic Optimum, which can be explained both by delays in local deglaciation momentum and topoclimatic conditions. Given the relatively small number of samples used in our study (38), this first attempt to assess the RSF evolution in the Romanian Carpathians might also be biased towards the Early Holocene by the high frequency of the high-altitude samples (71% of the samples are >1700 m). As future research, it is necessary to expand the RSF dating by including more cases from the mid and low levels and to compare their histories in order to disentangle the influence of deglaciation, permafrost thawing, thermal and humidity variation. However, some of the European studies describe some similar results with the newly-obtained in the Romanian Carpathians, such as [START_REF] Hermanns | Rapid rock-slope failures[END_REF] which give an estimation of Holocene RSF magnitude in Storfjörden, Norway, showing that the earliest events (12.5 to 10 kyrs) generated by far the largest detached volumes (Fig. 6), compared to the ones dating after 8 kyrs. Similarly, reconstructed magnitudes of large landslides from the Alps [START_REF] Soldati | Landslides and climate change in the Italian Dolomites since the Late glacial[END_REF][START_REF] Ivy-Ochs | Surface exposure dating of the Flims landslide, Graubünden, Switzerland[END_REF], place such events in the first millennia of the Holocene.

Independently of the absolute exposure ages used to assess the RSF probability occurrence during Holocene (Fig. 6), the other four surface ages from Doamnei rock glacier, in the central Făgăraș massif, indicate multiple phases of debris accumulation, and, in the same time, the high magnitude of early-Holocene debris production, which supplied at least the lower half of the rock glacier body, demonstrated by the rapid accumulation of massive boulders as well as their displacement between 12.97 and 9.91 ka ago. RW permafrost decay would have further enhanced subsequent rockfall or rock avalanches of smaller magnitude during the following warm episodes of the Holocene as also described in the Swiss Alps by [START_REF] Nagelisen | Post-glacial rock avalanches in the Obersee Valley, Glarner Alps, Switzerland[END_REF].

The particularly large boulders in the lower half of Doamnei rock glacier are incompatible with frost-cracking intensities estimated for the Holocene [START_REF] Savi | Efficiency of frost-cracking processes through space and time: An example from the eastern Italian Alps[END_REF] and were most probably produced by similar high magnitude (as described by [START_REF] Hermanns | Rapid rock-slope failures[END_REF] slope failures induced by slope relaxation, permafrost degradation and overall weakened slope coherence. We consider this to be the ultimate process/interval of intense modelling in the alpine area of the Carpathians which defined the rock wallsrock glaciers/debris systems preserved until present. Carpathians (blue line), and multi-site composed RSF distribution (red line) using absolute ages from the French Alps [START_REF] Cossart | Slope instability in relation to glacial debuttressing in alpine areas (Upper Durance catchment, southeastern France): Evidence from field data and 10Be cosmic ray exposure ages[END_REF]), Swiss Alps (Ivy-Ochs et al., 2009), Italian Alps [START_REF] Hormes | 10Be exposure ages of a rock avalanche and a late glacial moraine in Alta Valtellina, Italian Alps[END_REF]), Central Andes (Fauque et al., 2009), Northern Iceland [START_REF] Mercier | The Höfճ ahólar rock avalanche (sturzström): Chronological constraint of paraglacial landsliding on an Icelandic hillslope[END_REF], Central Karakoram [START_REF] Shroder | The role of mass movements on landscape evolution in the Central Karakoram: Discussion and speculation[END_REF], Scotland and NW Ireland [START_REF] Ballantyne | Enhanced rock-slope failure following ice-sheet deglaciation: Timing and causes[END_REF]. Grey-filled columns represent estimated volume of Holocene RSF in Norway (after [START_REF] Hermanns | Rapid rock-slope failures[END_REF]

CONCLUSIONS

The distribution of RW mapped in the Southern and Eastern Carpathians depends mainly on the lithology, structure but also weathering processes. In the metamorphic and igneous units, it ultimately relates to geomorphological context, being mostly associated with glacial cirques and valleys headwalls. Most RW are therefore restricted to the highest ridge sectors, while their low heights and areas are explained as a consequence of the lithological predisposition to debris production, especially in permafrost degrading conditions during warming phases of the Early Holocene. The North/South asymmetry in rock glaciers distribution (also signalled in other midlatitude ranges) is reflected by the lack of South-exposed RW. We assume more stable conditions prone to fine debris and soil formation on the southern slopes due to insolation and warmer conditions.

For the sedimentary RW, tectonics and the geological structure are the main controls to explain the occurrence of the large (wide and high) limestone and conglomerate RW in the Romanian Carpathians. Except for Bucegi Mts, which were high and large enough to host complex glaciers during the last glaciation, most of the sedimentary units from SC and EC were not subject to glacial erosion during LGM due to either steep topography (e.g., hogback ridges) or lower altitude, although RW permafrost was probably widespread.

Absolute exposure ages confirm that an intense rock slope degradation via rock-slope failures took place in the Carpathian metamorphic and igneous units in Early Holocene, similar with other European sites, reaching the highest magnitudes 11.6 -9 ka ago especially above 1800 m altitude. We associate the present distribution of RW with this periglacially-active period which was the last time of rock surfaces substantial reshape.

Figure 1 :

 1 Figure 1: (a) Location of Eastern and Southern Carpathians and Apuseni Mts., and outline of the units included in the RW inventory (digital elevation by 1 arc-second resolution ASTER GDEM). (b) Simplified geological map of the Carpathians (modified and adapted from Vaida and Verniers, 2005, and Merten, 2011. Small insert shows the Carpathian Arc position in Romania (please refer to the colour version)

Figure 2 :

 2 Figure 2: Direction and mean altitude of the RW mapped in the EC: (a) the metamorphic schistsprevailing units; the andesitic and basaltic outcrops are represented as volcanic RW; (b) the sedimentary units (limestone and conglomerates prevailing). The radius of the graphs represents the altitude values, the general direction is expressed in sexagesimal degrees and each dot represents a RW (please refer to the colour version)

Figure 3 :

 3 Figure 3: Direction and mean altitude of the RW mapped in the SC: (a) Făgăraș and Iezer Mts., with the secondary cluster of RW mapped in the lower altitude SC metamorphic units; (b) distribution of igneous and metamorphic rock surfaces in Parâng and Retezat Mts.; (d) distribution of limestone, conglomerate and sandstone RW mapped in SC; (d) Cumulated RW distribution on the four main rock categories represented on 30 degrees direction bins (please refer to the colour version)

Figure 4 :

 4 Figure 4: General RW morphometry on the main rock types (counting all features corresponding to the same rock type): (a) Mean RW altitude vs mean RW area; (b) Cumulated area of RW with the same lithology grouped on the four main cardinal directions; (c) Mean RW height and area derived for each massif (please refer to the colour version) Metamorphic and igneous units are relatively similar in terms of RW vertical extension (around 70 m in average) and show little inter-site variation whereas the sedimentary units both for

Figure 5 :

 5 Figure 5: Variability of the morphometry of the RW grouped in the main lithological classes in respect to the main orientations: (a) mean RW area; (b) mean height; (c) mean altitude. The boxes display the median values, the 25-75 quartiles (lowerupper) and the caps show the minimum and maximum values (1.5 IQR)

Figure 6

 6 Figure 6 Cumulated distribution of surface exposure ages attributed to post-Younger Dryas RSF in the Romanian

Table 1 :

 1 Main lithological, structural and morphographic characteristics of the mountain units in which rock walls were

		Unit name (max. altitude)	Lithology/ Structure*	Extent / Direction	Nr. of mapped RW	Glaciation**
		Maramureş (1956 m)	Crystalline schist with peripheral limestone and sandstones Volcanic intrusions (basalts)	15 km long / NW-SE ridge	11	Yes
	Eastern Carpathians	Rodna (2303 m) Călimani (2100 m) Ceahlău (1969 m) Hăşmaş (1973 m)	Crystalline schist, micaschists and paragneiss Horst, Dragoş Vodă Fault Andesites (volcanic) Eroded craters Conglomerates and flysch Suspended syncline Massive limestone West oriented syncline	40 km long / E-W ridge Volcanic cone 15 km / N-S 3.5 km long / NW-SE ridge	19 5 10 6	Yes Yes No No
		Ciucaş (1954 m)	Conglomerates and sandstones, flysch East oriented syncline	Two separated ridges 7 km / SW-NE and 3 km / NW-SE	37	No
			Conglomerates, sandstones
		Bucegi (2505 m)	Limestone with radiolarites N-S syncline, East oriented cuesta	Reversed U-shape 30 km long ridge	49	Yes
			front slope
		Piatra Craiului (2238 m)	Limestone with radiolarites Hogback	25 km long / NNE-SSW ridge	13	No
		Iezer (2459 m)	Micaschists and paragneiss Supragetic overthrust nappe	20 km long / SW-NE ridge	27	Yes
			Micashists and paragneiss,
		Făgăraş (2544 m)	amphibolite Supragetic overthrust nappe	70km long W-E ridge, and multiple secondary N-S ridges	248	Yes
			Northern-Făgăraş Fault Line
	Southern Carpathians	Cozia (1668 m) Buila-Vânturariţa (1885 m) Parâng (2519 m) Şureanu/ Cindrel/ Căpăţânii/ Lotrului (2130 -2244 m)	Gneiss Horst Massive Limestone Hogback Granitoids Amphibolite (Danubian Unit) Micaschists and paragneiss, amphibolite (Getic Unit)	~70 km 2 surface 14 km long / SW-NE ridge 25 km long / E-W ridge 15-25 km long / E-W ridges	3 23 40 8 / 10/ 4 / 5	No No Yes Yes
		Retezat (2509 m)	Granodiorite intrusions Crystalline schist, amphibolite (Danubian Unit)	15 km long main / W-E ridge 2-5 km secondary N-S ridges	187	Yes
		Ţarcu (2196 m)	Conglomerates, sandstones, crystalline limestone Crystalline schist	20 km long / N-S then NE-SW ridge	26	Yes
		Godeanu/ Piule Iorgovanul (2291 m)	Micaschists and paragneiss, amphibolite / Recifal limestone	20 km long / NE-SW ridge	26 / 13	Yes
		Cerna Valley (1200 m)	Recifal limestone Graben	80 km long valley / NS	21	No
	mapped
	* according to the Geological Map of Romania, scale 1:200 000 (Geological Institute of Romania)
	** according to the Map of Glacial Cirques in the Romanian Carpathians

  2 Excel add-in of[START_REF] Vermeesch | CosmoCalc: An Excel add-in for cosmogenic nuclide calculations[END_REF]. Air pressure used is 1013 mBar. There was no quantitative information on the snow cover during the surface exposure duration, hence, no corrections for potential effects of snow cover or denudation were applied to the ages. 10 Be exposure ages were calculated following Equation (2) using muogenic contributions of[START_REF] Braucher | Production of cosmogenic radionuclides at great depth: A multi element approach[END_REF]

	𝑁(𝑥, 𝜀, 𝑡) =	𝑃𝑠𝑝 * exp (-	𝑥 𝐿𝑛	)(1 -exp (-𝑡 ( 𝐿𝑛 𝜀 𝜀 𝐿𝑛 + 𝜆	+ 𝜆)
	+	𝑃𝜇𝑠𝑙𝑜𝑤 * exp (-	𝑥 𝐿𝜇𝑠𝑙𝑜𝑤 𝜀 𝐿𝜇𝑠𝑙𝑜𝑤 )(1 -exp (-𝑡( 𝐿𝜇𝑠𝑙𝑜𝑤 𝜀 + 𝜆	+ λ)
	+	𝑃𝜇𝑓𝑎𝑠𝑡 * exp (-	𝑥 𝐿𝜇𝑓𝑎𝑠𝑡) )(1 -exp (-𝑡( 𝐿𝜇𝑓𝑎𝑠𝑡 𝜀 𝜀 𝐿𝜇𝑓𝑎𝑠𝑡 + 𝜆	+ 𝜆)
						Equation (2)
	where:

N (x, ε, t) is the nuclide concentration function of depth x (g/cm 2 ), denudation rate ε (g/cm 2 /y) and exposure time t(y). Psp, Pμslow, Pμfast and Ln, Lμslow, Lμfast are the production rates and attenuation lengths of neutrons, slow muons and fast muons, respectively. Ln, Lμslow, Lμfast values used are 160, 1500 and 4320 g/cm 2 , respectively

[START_REF] Braucher | In situ produced 10Be measurements at great depths: implications for production rates by fast muons[END_REF]

. λ is the radioactive decay constant. Pμslow, Pμfast are based on

[START_REF] Braucher | Production of cosmogenic radionuclides at great depth: A multi element approach[END_REF]

Table 2 :

 2 The averaged values of the morphometric parameters and total cumulated area of the mapped RW in EC and SC

  The highest mean altitude is recorded in Piatra Craiului (1970 m) while Cerna Valley reaches the lowest (771 m).Bucegi is the highest sedimentary massif (2507 m) from Carpathians, with most RW developed on conglomerates and sandstones, and 15% on limestone outcrops. RW are distributed mainly on the northern slopes (37%) and only 16% are South-facing. Maximum RW density is between 2100 and 2300 m for the North and East-exposed slopes. East, South and West RW are situated at lower elevations and occur on a wider altitudinal range (1800 -2200 m).

Table 4 :

 4 Results of the statistical analysis (Post-Hoc Dunn's Test) of morphometric parameters for pairs of main exposures. Numbers represent the z-Score, which indicates whether the tested parameter pair has a value above the

		Area	Height	Altitude
	Conglomerate/ Limestone	-1.78 (*)	-3.9 (***)	2.93 (**)
	Granite/ Limestone	-7.54 (***) -6.9 (***)	12.89 (***)
	Schist/ Limestone	-9.01 (***) -8.36 (***) 10.9 (***)
	Granite/ Conglomerate	-5.99 (***) -2.79 (**)	10.31 (***)
	Schist/ Conglomerate	-7.46 (***) -3.87 (***) 7.91 (***)
	Schist/ Granite	-0.77 (•)	-0.83 (•)	-4.53 (***)
	rank mean (positive value), or below (negative value)
	* p-value (0.05, 0.01]; ** p-value (0.01, 0.001]; *** p-value < 0.001. A p-value < 0.05 indicates a strong statistical
	significance at 95% confidence level

Table 5 :

 5 Sampling locations, 10 Be concentrations, and 10 Be surface exposure ages for post-Younger Dryas dated boulders in the Romanian Carpathians

	a Analytical uncertainties (reported as 1-σ) included for all samples. No corrections for potential effects of snow cover or denudation
	were applied to the ages
	b Exposure ages calculated using Cronus-Earth 10Be -26Al exposure age calculator v. 2.2 (http://hess.ess.washington.edu/). They
	assume zero erosion, scaling factors according to

  . Exposure ages are presented with the external uncertainties The measured 10Be/9Be AMS ratios were corrected for full processed blank ratios: (3.30±0.50) x 10 -15 . Age uncertainties: the 1st number is the internal uncertainty (AMS measurement, weighting, carrier, blank and half-life; 1-σ). Every reported age was corrected for topographic-and self-shielding

	RETEZAT
									Ruszkiczay-
	Lăpușnicu Valley	Re15-29	45.31	22.78	1167	-	10.30 ± 0.13	10.70 ± 1.40	Rüdiger et al.,
									(2021) c
	Pietrele Valley	Pt-03-02	45.28	22.88	1902	-	23.9	11.40 ± 1.30	Reuther et al. (2007) d
	RODNA
	Pietroasă Valley	RD 30	47.61	24.64	1379	25.11	16.32 ± 0.50	10.50 ± 0.90	Gheorghiu, (2012) ee
	Zănoaga Mare cirque	RD 04	47.6	24.64	1669	28.95	21.65 ± 0.63	11.50 ± 1.00
	Zănoaga Mare cirque	RD 06	47.6	27.63	1767	24.24	22.83 ± 0.63	11.30 ± 1.00
	Zănoaga Mare cirque	RD 07	47.6	27.63	1767	25.85	22.50 ± 0.65	11.10 ± 1.00
	Zănoaga Mare cirque	RD 05	47.6	27.63	1753	24.36	10.89 ± 0.40	5.70 ± 0.50
	Buhăiescu Valley	RD 19	47.58	24.65	1718	23.01	21.61 ± 0.63	10.40 ± 0.90

c

d Exposure age corrected for the effect of topographic shielding and surface geometry e Exposure ages calculated using Cronus-Earth 10Be -26Al exposure age calculator v. 2.2 (http://hess.ess.washington.edu/). They assume zero erosion, scaling factors according to

[START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF]

Acknowledgements: The research financing these results received funding from the Norway Grants 2014-2021 (project code RO-NO-2019-0415 / contract no. 30/2020), from University of Bucharest grants (project code C1.2.PFE-CDI.2021-587) and from the Research Institute of the University of Bucharest fellowship for young researchers (2018-2019 grants, awarded to MV).