Emmanuel Giner

come

A new form of transcorrelated Hamiltonian inspired by range-separated DFT

A new form of transcorrelated Hamiltonian inspired by range-separated DFT Emmanuel Giner a) The present work introduces a new form of explicitly correlated factor in the context of the transcorrelated methods. The new correlation factor is obtained from the r 12 ≈ 0 mathematical analysis of the transcorrelated Hamiltonian, and its analytical form is obtained such that the leading order in 1/r 12 of the scalar part of the effective two-electron potential reproduces the long-range interaction of the range-separated density functional theory. The resulting correlation factor exactly imposes the cusp and is tuned by a unique parameter µ which controls both the depth of the coulomb hole and its typical range in r 12 . The transcorrelated Hamiltonian obtained with such a new correlation factor has a straightforward analytical expression depending on the same parameter µ, and its physical contents continuously change by varying µ : one can change from a non divergent repulsive Hamiltonian at large µ to a purely attractive one at small µ. We investigate the convergence of the ground state eigenvalues and right-eigenvectors of such new transcorrelated Hamiltonian as a function of the basis set and as a function of µ on a series of two-electron systems. We found that the convergence towards the complete basis set is much faster for a quite wide range values of µ. We also propose a specific value of µ which essentially reproduce the results obtained with the frozen Gaussian geminal introduced by Ten-No [CPL-330,169 (2000)].

I. INTRODUCTION

One of the most challenging problem in computational chemistry is the accurate simulation of the electronic structure of atomic and molecular systems, which implies solving the Schroedinger equation for a general N-body system. At the heart of the complexity of such a task relies the rapidly prohibitive computational cost of the wave function methods (WFT), which is induced by the exponential growth with the system size of the Hilbert space involved in the linear eigenvalue problem to be solved.

In usual WFT, the computational bottleneck is mainly determined by two factors: the level of complexity of the wave function which is imposed by the amount of strong correlation effects in the system, and the convergence of the computed quantities with respect to the size of the one-electron basis used to project the Hamiltonian into a finite eigenvalue problem. Typically, strong correlation effects appear when the two-electron coulomb interaction can no longer be considered as a small perturbation over a mean field Hamiltonian. Treating such effects in a black box way involves the use of very flexible wave functions such as selected CI [1][2][3][4][5][6][7][8][9][10][11] , full configuration Quantum Monte Carlo [12][13][14][15][16] (FCI-QMC) or matrix product states 17,18 . These methods all have one thing in common: they manage to select, although through different mathematical frameworks, the part of the exponentially growing Hilbert space which dominates the wave function.

Even if these schemes are extremely efficient to obtain a qualitative description of the wave function through a meaningful selection in an exponentially growing Hilbert space, a quantitative description is often out of reach as it implies the use of very large one-electron basis sets which are needed to recover the slowly converging correlation effects near the electron-electron coalescence point. These short-range correlation effects are mainly due to the diverging character of the coulomb potential, which induces a non analytical behaviour of the wave function at small inter-electronic distance (r 12 ≈ 0): a) Electronic mail: eginer@lct.jussieu.fr the famous electronic cusp originally derived by Kato 19 . For a clear and general derivation of the cusp conditions up to second order for Coulombic systems, see Ref. 20. As shown by the seminal works of Hylleraas 21 and further developed by Kutzelnigg and coworkers, [22][23][24] the slow convergence of the correlation energy with respect to the quality of the basis set is mainly due to the impossibility to represent the cusp in a finite one-electron basis set. As the cusp conditions originate from the divergence of the Coulomb potential at r 12 = 0, an alternative approach would be to develop theories dealing with a smooth potential. Such a path was followed by three distinct branches: the explicitly correlated methods (F12), transcorrelated methods (TC) and the range-separated density functional theory (RS-DFT).

In RS-DFT [START_REF] Savin | Recent Advances in Density Functional Theory[END_REF][START_REF] Toulouse | [END_REF] , the electron-electron interaction is split into a smooth long-range part and a complementary short-range diverging part, the former being treated explicitly by a wave function and the latter by a density functional. Such a splitting is done through the function erf(µr 12 ) where µ is the range separation parameter which allows one to continually move from Kohn-Sham DFT (µ = 0) to pure WFT (µ = ∞). As the effective Hamiltonian obtained in RS-DFT is smooth at r 12 = 0 for a finite value of µ, the corresponding eigenfunction has no electron-electron cusp 27 and therefore the convergence of the results with respect to the basis set are exponential 28 . A number of approximate RS-DFT schemes have been developed involving single-reference [29][30][31][32][33][34][35][36] and multi-reference [37][38][39][40][41][42] WFT methods. Nevertheless, there are still some open issues in RSDFT, such as remaining fractional-charge and fractional-spin errors in the short-range density functionals 43 or the dependence of the quality of the results on the value of the range-separation parameter µ.

Another approach to tackle the problems of WFT related to the electron cusp have been proposed by the F12 methods which introduce a function explicitly depending on the interelectronic coordinate [44][45][46][47][48][49] (called a Jastrow factor) to describe short-range correlation effects which are absent from the finite basis set. Thanks to strong orthogonality, all redundant correlation effects between the Jastrow factor and the basis set are projected out of the wave function, and therefore the Jastrow factor only deals with correlation effects orthonormal to the basis set. The resulting F12 theories, mainly based on perturbation theory and coupled cluster theory, are therefore extremely close to their parent WFT theory, with additional contributions involving two-and three-body terms resulting from the use of the Jastrow factor. It is important to notice that because of the orthogonalization between the Jastrow and the basis set, the complexity of the N-body problem within a given basis set is essentially the same between a F12 theory and the parent theory in WFT. Therefore, no strong compression of the electronic wave function is obtained within a given basis set by the F12 theory, but a fast convergence toward the complete basis set limit (CBS) is nevertheless obtained for correlation energies and other properties sensitive to the basis set.

An alternative point of view using a Jastrow factor is the so-called transcorrelated (TC) approach where the full effect of the Jastrow factor is incorporated into the calculation through a similarity transformation of the original Hamiltonian by the Jastrow factor. Seminal equations were derived by Hirschfelder 50 who obtained an effective non hermitian operator based on a specific form of Jastrow factor, and later on Boys, Handy and co-workers [51][52][53] derived the equations of the TC Hamiltonian for a general exponential form of the Jastrow factor which insures size extensivity. The obtained TC Hamiltonian contains certain new features with respect to the usual Hamiltonian: the TC Hamiltonian contains an additional effective two-and three-electron scalar potential together with is a non hermitian two-electron differential operator. Because of its non-hermitian nature, the TC Hamiltonian looses the variational principle, and Boys and Handy managed to derive equations to optimize both orbitals and Jastrow parameters for a single Slater determinant. Nevertheless, such a constrained form for the Slater part of the wave function together with the loose of variational principle makes it difficult to find a monotonic convergence of computed energies 54 .

Later on, Ten-No 55 proposed to significantly change the paradigm of Boys and Handy: instead of using a single Slater determinant and optimizing both the orbitals and the Jastrow factor, one uses a more elaborate many body theory with a frozen universal Jastrow factor. This strategy was developed using as Møller-Plesset at second order (MP2) in Refs. 55 and 56 and a linearised coupled cluster ansatz in later work 57 . The Jastrow factor introduced by Ten-No, developed as a linear combination of gaussian functions and referred to as the frozen gaussian geminal (FROGG), only depends on the interelectronic coordinate and is optimized such that the scalar effective two-electron potential cancels on average the coulomb potential near r 12 = 0. The use of gaussian functions to represent the Jastrow factor enabled Ten-No to develop a numerical scheme to evaluate exactly the three-electron integrals present in TC Hamiltonian 58 , which is an alternative to the densityfitting method proposed by the same author in the general context of F12 methods 59 . In Ref. 56, Ten-No et. al. introduced a biorthogonal approach, which allows to treat the non-hermitian TC Hamiltonian in a mathematical form more suited to develop approaches based on many-body perturbation theory. It is also noteworthy that the eigenfunctions of the TC Hamiltonian are invariant by orbital rotation, and therefore in the limit where one gives full flexibility to the wave function, the biorthogonal and usual orthogonal approaches give the same results. Nevertheless, as approximations and truncations are necessary done in the wave function ansatz for realistic systems, the biorthogonal approach might improve the quality of the results in practice.

In later works, Umezawa et. al. 60,61 introduced a scheme to couple TC equations with the usual variational Monte Carlo (VMC) scheme in order to optimise the Jastrow factor through a variance miminization. An attempt to make TC equations variational have been proposed by Luo 62,63 based on the empirical experience that in the original single determinant TC approach of Boys and Handy, the orbital optimization seemed to be the source of loss of the variational property. Luo therefore proposed to replace the original TC equations to optimize the orbitals by that of a general Jastrow-Slater ansatz in VMC. An alternative to the non hermitian nature of TC theory have been proposed by Yanai et. al. 64 where the Jastrow factor is expressed as an anti-hermitian operator, allowing one to obtain an hermitian effective Hamiltonian after the similarity transformation, which is then efficiently truncated through the use the canonical transformation 65 . Applications of such method in the context of quantum computing have been recently published 66 .

Further developments of the TC method towards the treatment of solid state systems have been carried by Ochi et. al. [67][68][69][70] , which includes both ground states calculations within a single determinant wave function 67,70 or at the MP2 level 69 , together with excited states within configuration interaction 68 . The TC framework has also been used in the context of DFT to develop new approximations of density functionals [71][72][73] .

More recently, Alavi et. al 74 applied the TC equations with the use of elaborate Jastrow factors which explicitly takes into account the electron-electron together with electron-electronnucleus coordinates. The Jastrow factor used in such work were obtained from the optimized Jastrow factors for He to Ne of Moskowitz et. al. 75 in the context of quantum Monte Carlo. In contrast to previous works where a constrained form was given to the wave function (such as coupled cluster for instance 57 ), Alavi et. al. allow a full flexibility to the Slater part of the wave function to adapt to the elaborated Jastrow factors: they fully solve the N-electron problem corresponding to the TC Hamiltonian within a given basis set. This is done through the use of the FCI-QMC [12][13][14][15][16] which is a projective technique allowing one to obtain the ground state of an operator, hermitian or not, through a stochastic sampling of the corresponding Hilbert space. The use of a TC Hamiltonian in the context of FCI-QMC has an important advantage: because the TC Hamiltonian already contains the effect of an elaborate Jastrow factor, the right-eigenvectors of the TC Hamiltonian are more compact and therefore the FCI-QMC procedure converges faster. This approach share common points with that of Ten-No in the sense that they both use a fixed Jastrow factor and use a flexible form of wave functions, even though the FCIQMC is certainly more flexible than the linearised coupled cluster ansatz of Ten-No used in Ref. 56. In the present work, we derive a new form of Jastrow factor such that the leading 1/r 12 terms of the corresponding TC Hamiltonian reproduce the RS-DFT effective Hamiltonian. As in RS-DFT, the new TC Hamiltonian and Jastrow factor are tuned by a unique range separation parameter µ. The aim of this paper is to establish the analytical form of the TC Hamiltonian and perform a numerical study on the ground state eigenvalues and eigenvectors on a set of two-electron atomic and molecular systems: the helium isoelectronic series from H -to Ne 8+ and H 2 molecule. In the context of the TC theory the interesting features of the present work are that i) the analytical form of the TC Hamiltonian is explicitly known and it turns out that all two-electron integrals can be performed analytically and the three-electron integrals can be efficiently computed in mixed numerical-analytical way, ii) the new Jastrow factor and corresponding TC Hamiltonian are tuned by a unique parameter µ which allows some flexibility, iii) there always exists a regime of µ which significantly improves the basis set convergence of WFT, iv) we propose schemes to find a value of µ which automatically adapts to the system and improves basis set convergence.

The paper is organized as follows. In the first section, inspired by the RS-DFT effective Hamiltonian, we derive the analytical form of a new Jastrow factor tuned by a unique parameter µ, and discuss the physical content of such a new Jastrow factor. Then, we show that the corresponding TC Hamiltonian has a straightforward analytical structure, which is briefly reviewed as a function of the µ parameter. We perform a numerical study of the behaviour of the ground state eigenpair of H[µ] in the case of the helium atom: in the first part we analyze the speed of convergence of the ground state energy with respect to the basis set and the parameter µ, then we study the behaviour of the right eigenvectors in real space and compare it with the numerically exact ground state wave function. We show that the rapid convergence of the total energies coincides with a rapid convergence of the right eigenvectors in real space, and how to recover a very good approximation of the exact ground state wave function. Based on such encouraging results, we study the ground state energy of the helium isoelectronic series from H -to Ne 8+ and compare different schemes to obtain a value of µ which automatically adapts to the system and still provides a fast convergence of the total energies. Eventually, we study the H 2 molecule and compare different regimes of correlation. Numerical comparison with the TC Hamiltonian obtained with the FROGG is also performed for all systems studied.

II. A NEW FORM OF JASTROW FACTOR FOR TRANSCORRELATED HAMILTONIANS

This section is dedicated to the derivation of a new form of Jastrow factor mimicking at short range the effective Hamiltonian of RS-DFT. For the sake of simplicity, we derive the main equations in the case of the helium atom which have an explicit form in terms of the r 12 coordinate. Then, we derive the analytical form of the corresponding new transcorrelated Hamiltonians H[µ] in the case of a general N-electron system, and study its physical content as a function of the parameter µ.

H = h c + 1 r 12 , (1)
where

h c = - 1 2 2 i=1 ∂ 2 ∂r i 2 + 2 r i ∂ ∂r i + 2Z r i - ∂ 2 ∂r 12 2 + 2 r 12 ∂ ∂r 12 - r 1 r 1 • r 12 r 12 ∂ ∂r 1 + r 2 r 2 • r 21 r 21 ∂ ∂r 2
.

(

) 2
Now let us consider the transcorrelated Hamiltonian H[u] obtained by the similarity transformation of the usual Hamiltonian by a Jastrow factor u(r 12 ) depending only on r 12 :

H[u] = e -u(r 12 ) He u(r 12 ) . (3)
With respect to H, the only additional terms arising in H[u] are those coming from the action of the differential operator in r 12 , ,

T [u] = -e -u(
W[u] = - 2 r 12 ∂u(r 12 ) ∂r 12 , (5)
[u] = - ∂ 2 u(r 12 ) ∂r 12 2 - ∂u(r 12 ) ∂r 12 2 , (6) w
one can write the Similarity transformed Hamiltonian as

H[u] = h c + 1 r 12 + W[u] + w[u] + t[u]. (8)
Therefore, one can see from Eq. ( 8) that the two electrons experience a modified scalar potential given by

which is equivalent to

∂u(r 12 , µ) ∂r 12 = 1 -erf(µr 12 ) 2 . ( 10
)
The solution to Eq. ( 10) is

u(r 12 ; µ) = 1 2 r 12 1 -erf(µr 12 ) - 1 2 √ πµ e -(r 12 µ) 2 , (11)
which defines the new Jastrow factor u(r 12 ; µ) which depends on a unique parameter µ.

The main differences of the present approach with the FROGG introduced by Ten-No 55 is that instead of optimizing the Jastrow factor as a linear combination of Gaussian functions such that 1

r 12 + W[u] + w[u]
≈ 0 when r 12 = 0 through a least-square fit with a specific weighting function, we impose the analytical structure of u(r 12 ; µ) such that i) the sum of all terms in 1/r 12 (i.e. W[u] + 1 r 12 ) exactly provides the non divergent long-range interaction erf(µr 12 )/r 12 for all r 12 , ii) the new Jastrow factor u(r 12 ; µ) is more flexible than the FROGG as it is tuned by a unique parameter µ which determines both the depth of the hole and its typical range in r 12 .

A few properties of u(r 12 ; µ)

The new Jastrow u(r 12 ; µ) factor is tuned by a unique parameter, µ, which have the unit of the inverse of a distance. From Eq. ( 11) we can notice that lim

r 12 →∞ u(r 12 ; µ) = 0, (12)
which means that lim

r 12 →∞ e u(r 12 ;µ) = 1, (13)
and therefore u(r 12 ; µ) impacts only the small r 12 behaviour, which is expected since the equation determining u(r 12 , µ) is obtained from the small r 12 limit of the Schroedinger equation (see Eq. ( 9)).

To better understand the physical content of at r 12 ≈ 0, let us Taylor expand the function u(r 12 ; µ)

u(r 12 ; µ) = - 1 2 √ πµ + 1 2 r 12 - µ 2 √ π r 2 12 + O(r 4 12 ). (14)
A first remark is that u(r 12 ; µ) exactly restores the cusp as lim

r 12 →0 ∂ ∂r 12 u(r 12 ; µ) = 1 2 , (15)
implying that the right-eigenvectors of the TC Hamiltonian obtained with u(r 12 ; µ) are necessary cusp less (see Appendix IX D). The Taylor expansion of Eq. ( 14) inserted in the exponential correlating Jastrow factor reads e u(r 12 ,µ) ≈ e

-1 2 √ πµ e 1 2 r 12 e -µ 2 √ π r 2 12 . (16)
From Eq. ( 16) one can analyze each term: i) as e Coming now to the behaviour of u(r 12 ; µ) as a function of the parameter µ, we report in Figure 1 its dependence in r 12 for a set of values of µ, and we compare it with the FROGG. From Figure1, one can see that: i) the larger the µ, the more short-range is u(r 12 ; µ), ii) the larger the µ, the closer u(r 12 ; µ) is from 0 and therefore the less impact the correlating factor has on the wave function, iii) the smaller the µ, the deeper and broader is the hole imposed by e u(r 12 ,µ) .

As the function u(r 12 ; µ) contains a unique parameter µ, one can try to find the value µ FROGG such that u(r 12 ; µ FROGG ) qualitatively reproduces the FROGG. To do so, we impose that µ FROGG is such that the two Jastrow factors coincide at coalescence, i.e.

u(r 12 = 0; µ FROGG ) = u TN (r 12 = 0), (17)
where u TN (x) is the FROGG introduced by Ten-No

u TN (x) = - 6 i=1 c i e -α i x 2 , (18)
with the coefficients c i and exponents α i are given in Table 1 of Ref. 55. Eq. ( 17) gives

µ FROGG = - 1 2 √ πu TN (0) , (19)
which gives µ FROGG ≈ 0.87. We represent on Fig. 1 the comparison between u(r 12 ; µ FROGG ) and the FROGG, and it is clearly apparent that u(r 12 ; µ FROGG ) reproduces almost perfectly the FROGG, although it decays slightly faster at large r 12 . Therefore, one can see the FROGG as a special case of the new Jastrow factor u(r 12 ; µ) with µ = µ FROGG .

B. The general analytical form of H[µ]

and some of its properties

General analytical form of H[µ] for a N-electron system

Having established the analytical expression of u(r 12 ; µ), one can derive the expression of the associated transcorrelated Hamiltonian for a general N-electron system. Using the general form of the transcorrelated Hamiltonian of Eq. (2) in Ref. 74 and using u(r 12 , µ) as Jastrow factor gives We also compare the FROGG (see Eq. ( 18)) and the value of µ = 0.865 which equals the FROGG at r 12 = 0 (see Eq. ( 17)).

H[µ] ≡ e -τ µ Ĥe τµ = H + H, τµ + 1 2 H, τµ , τµ , (20)
with τµ = i< j u(r i j , µ) and

Ĥ = i -1 2 ∇ 2 i + v(r i ) + i< j 1 r i j . Therefore, the i-th right-eigenvectors Φ ex i [µ]
depends on µ and satisfy the effective Schroedinger equation

H[µ] Φ ex i [µ] = E i Φ ex i [µ] , (21)
where E i is the exact i-th eigenvalue of the physical Hamiltonian. The fact that H[µ] and H share the same spectrum originates from the properties of the similarity transformation. Nevertheless, once projected in an incomplete basis set, H[µ] and H might not have the same eigenvalues. Eq. ( 20) leads to the following transcorrelated Hamiltonian

H[µ] = H - i< j K(r i , r j , µ) - i< j<k L(r i , r j , r k , µ), (22)
where the effective two-and three-body operators K(r 1 , r 2 , µ) and L(r 1 , r 2 , r 3 , µ) are defined as

K(r 1 , r 2 , µ) = 1 2 ∆ 1 u(r 12 , µ) + ∆ 2 u(r 12 , µ) + ∇ 1 u(r 12 , µ) 2 + ∇ 1 u(r 12 , µ) 2 +∇ 1 u(r 12 , µ) • ∇ 2 + ∇ 2 u(r 12 , µ) • ∇ 1 , (23)
L(r 1 , r 2 , r 3 , µ) =∇ 1 u(r 12 , µ) • ∇ 1 u(r 13 , µ) +∇ 2 u(r 21 , µ) • ∇ 2 u(r 23 , µ) +∇ 3 u(r 13 , µ) • ∇ 3 u(r 32 , µ). ( 24
)
As shown in the Appendix IX A, the operator K(r 1 , r 2 , µ) has a straightforward analytical expression which reads

K(r 1 , r 2 , µ) = 1 -erf(µr 12 ) r 12 - µ √ π e -µr 12 2 + 1 -erf(µr 12 ) 2 4
erf(µr 12 ) -1 ∂ ∂r 12 ,

and similarly the analytical expression of L(r 1 , r 2 , r 3 , µ) reads

L(r 1 , r 2 , r 3 , µ) = 1 -erf(
By injecting the form of Eqs. ( 25) and ( 26) in Eq. ( 22), one obtains

H[µ] =h c + i< j Wee (r i j , µ) + t[µ] - i< j<k L(r 1 , r 2 , r 3 , µ) (27)
where Wee (r i j , µ) is the scalar effective two-body potential of

H[µ] Wee (r i j , µ) = erf(µr i j ) r i j + µ √ π e -µr i j 2 - 1 -erf(µr i j ) 2 4 , (28
) and t[µ] is a non hermitian two-body differential operator

t[µ] = erf(µr i j ) -1 ∂ ∂r i j . ( 29
)
As the present study is devoted to two-electron systems, the physics of H[µ] can be qualitatively understood by looking at Wee (r i j , µ), which is illustrated in Fig. 2 for different values of µ and compared with the usual long-range interaction of the RS-DFT framework. From Fig. 2 we can see that the scalar effective potential Wee (r i j , µ) is non divergent, increases when µ increases and globally resembles the long-range interaction used in RS-DFT at large r 12 . Nevertheless, as observed from Fig. 2, it is significantly different from the long-range interaction at small r 12 , and it is not monotonic and attractive (i.e. lim r 12 →0 ∂ ∂r 12 Wee (r 12 , µ) > 0). This difference with respect to the long-range interaction can be understood as the equation obtained to derive the Jastrow factor only takes into account the first-order derivative of u(r 12 , µ) (see Eqs. ( 9) and ( 10)). If one desires to impose that Wee (r 12 , µ) = erf(µr = 1erf(µr 12 ), (30) whose solution is unknown to the best of the author knowledge.

Limit cases for H[µ] with µ

Regarding now the variation of H[µ], as

lim µ→∞ u(µ, r 12 ) = 0, (31)
one recovers the usual Hamiltonian in the µ → ∞ limit:

lim µ→∞ H[µ] = H, ( 32
)
although it is only in the sense of distributions (see Appendix IX E for a more detailed explanation). In the µ → 0 limit, the transcorrelated Hamiltonian H[µ] becomes simply

lim µ→0 H[µ] = h c - 1 4 - i< j ∂ ∂r i j - 1 4 i< j<k r i j r i j • r ik r ik + r ji r ji • r jk r jk + r ki r ki • r k j r k j , (33)
which means that all scalar electron-electron repulsive interaction have been replaced by a differential operator and a three electron interaction. It should be stressed that even if the µ → 0 limit of H[µ] is very unphysical, it conserves exactly the same spectrum than the usual Hamiltonian because it originates from a similarity transformation.

A special case of H[µ] which mimics the FROGG

As the Jastrow factor u(r 12 , µ) for µ = µ FROGG is very close to the FROGG, we can also compare the scalar effective potential Wee (r 12 , µ FROGG ) with that obtained with the transcorrelated Hamiltonian obtained with the FROGG. The latter is simply defined as

HTN = e -τ TN He τTN , (34)
with τTN = i< j u TN (r i j ) and u TN (r i j ) given in Eq. ( 18). This leads to the following transcorrelated Hamiltonian

HTN = H + i< j WTN ee (r i j ) + tTN - i< j<k LTN (r 1 , r 2 , r 3 ) (35)
where

WTN ee (r 12 ) = 1 r 12 -∆ 1 u TN (r 12 ) -∇ 1 u TN (r 12 ) 2 , (36)
tTN

=∇ 1 u TN (r 12 ) • ∇ 2 + ∇ 2 u TN (r 12 ) • ∇ 1 , (37)
LTN (r 1 , r 2 , r 3 ) =∇ 1 u TN (r 12 ) • ∇ 1 u TN (r 13 )

+∇ 2 u TN (r 21 ) • ∇ 2 u TN (r 23 ) +∇ 3 u TN (r 31 ) • ∇ 3 u TN (r 32 ). ( 38
)
As apparent from Fig. 3 where WTN ee (r 12 ) and Wee (r 12 , µ FROGG ) are represented, while the latter is perfectly smooth for all r 12 , the former strongly oscillates near r 12 = 0 between positive and negative values and is not bounded for r 12 = 0. Such behaviour of WTN ee (r 12 ) near r 12 = 0 comes from the fact that, as it is made of gaussian functions, u TN (r 12 ) is cusp-less i.e. u TN (r 12 ) = 0, (39) and as a necessary condition to obtain a non divergent scalar potential is (see Eqs. ( 92) and (93) in Appendix IX D) lim

r 12 →0 ∂ ∂r 12 u(r 12 ) = 1 2 , (40)
WTN ee (r 12 ) is necessary not bounded at r 12 = 0. Nevertheless, as it will be shown by numerical computations, the effective potentials WTN ee (r 12 ) and Wee (r 12 , µ FROGG ) are almost identical in average, as they numerically agree for r 12 > 1 and that WTN ee (r 12 ) essentially oscillates around Wee (r 12 , µ FROGG ) for small values of r 12 . 36)).

III. DETAILED NUMERICAL STUDY OF H[µ] IN THE CASE OF THE HELIUM ATOM

The first part of our numerical study study concerns the convergence of the eigenvalues of H[µ] with respect to the basis set, and more precisely the dependence on µ of such convergence.

Considering a given basis set B and the corresponding projector P B , one can define the transcorrelated Hamiltonian projected on B by

H[µ] B = P B H[µ]P B , (41)
and its eigenvalues

E B i [µ] and right eigenvectors Φ B i [µ] satisfy H[µ] B Φ B 0 [µ] = E B i [µ] Φ B i [µ] . (42)
As long as the basis set B is incomplete, E B i [µ] might not coincide with E i , and because of its non hermitian nature, H[µ] B looses the variational principle and therefore its eigenvalues

E B i [µ]
have no reasons to be bounded from below.

A. µ dependence of the convergence of the ground state eigenvalue with respect to the basis set

Computational details

We implemented all integrals needed for the computation of the matrix elements of H[µ] on a usual Gaussian orbital basis set (see Appendix IX B for more details) together with the modification of the Slater rules due to the presence of the non hermitian term in H[µ]. We also implemented all necessary integrals to use the FROGG derived by Ten-No in Ref. 55 in the context of the TC Hamiltonian. All implementation was done as a plugin of the program quantum package 77 . We use the restricted Hartree-Fock (RHF) molecular orbitals to build all matrix elements of H[µ], and then we solve the two-body problem giving full flexibility to

Φ B i [µ] (i.e. Φ B i [µ]
is expressed as a linear combination of all possible Slater determinants within a given basis set B) with a non hermitian eigensolver to obtain both right and left eigenvectors together with the eigenvalues.

Numerical results: ground state energies

We report in Table I the ground state eigenvalue E B 0 [µ] of H[µ] B in the AVXZ basis sets (X=D,T,Q,5) for different values of µ, and we report in Figure 4 the graphical representation of such data. We also report in the tables and figures the results obtained with the FROGG. Several observations can be done from these data. First, because of the loose of the variational principle, E B 0 [µ] can take values below the exact non relativistic energy, and the smaller the µ, the more pronounced is such effect. Nevertheless, one can also observe that, for each µ, the error with respect to the exact energy gets smaller when enlarging the basis set, which illustrates that in the limit of a complete basis set the eigenvalues of any TC Hamiltonian coincide with that of the physical Hamiltonian.

Qualitatively, one can observe that there are two regimes of µ: µ ∈ [0. Also, it is striking to see that, in a given basis set, the ground state eigenvalue of the TC Hamiltonian obtained by the FROGG and that obtained by H[µ] with µ = µ FROGG coincide with less than 1 mH. The latter results confirms that as u(r 12 , µ FROGG ) ≈ u TN (r 12 ) (see Fig. 1 and Eqs. ( 17) and ( 19)), the TC Hamiltonians obtained with the two Jastrow factor are very similar, as illustrated for instance by the effective potential in Fig. 3.

Coming now to the speed of convergence of E B 0 [µ] with respect to the basis set, it is always faster than that of E B 0 [FCI] for the whole range of values of µ selected here. To quantify better such observations, one can identify the basis set from which the error with respect to the exact non relativistic energy is smaller, in absolute value, than 1mH. For 0.3 ≥ µ ≥ 0.4 one can see that the error with respect to E 0 is never larger, in absolute value, than 1mH whatever the basis set used, which represents a strong improvement with respect to the FCI. For 0.5 ≥ µ ≥ 1.0 the 1mH error threshold is reached at the augcc-pVTZ level and the accuracy are comparable to the FCI in the aug-cc-pV5Z basis set. For 1.6 ≥ µ ≥ 3.0, such accuracy is reached at the aug-cc-pVQZ level, whenever such accuracy is reached at the aug-cc-pV5Z level using the regular FCI Hamiltonian.

From the study of the ground state of the helium atom, one can conclude that i) the use of a simple Jastrow factor such as u(r 12 ; µ) within the transcorrelated framework can strongly improve the speed of convergence of the total energy of the helium atom with a quite wide range of values of µ, ii) the results obtained with µ = µ FROGG reproduce essentially that obtained with the FROGG, iii) for µ ∈ [0.35, 0.7] the results obtained are always better than that obtained with the FROGG.

B. Study of the ground state eigenvector in real space

In Sec. III A 2 we shown the improvement of the speed of convergence of the ground state energies E B 0 [µ] of the helium atom, which was effective for a quite large range of µ starting from µ = 0.2 to µ = 1.6. Nevertheless, one can wonder if this good behavior results from a kind of error cancellation or if there is indeed a deeper mathematical reason for such an improvement of speed of convergence. To bring insights to these observations, we investigate the behaviour in real space of the ground state eigenvectors of the TC Hamiltonian and the usual physical Hamiltonian in the case of the helium atom.

For a two-electron system, the exact ground state eigenvector

Φ ex 0 [µ] of H[µ]
is directly related to the exact ground state wave function Ψ ex 0 by the Jastrow factor u(r 12 , µ) through

Ψ ex 0 (r 1 , r 2 ) = 1 √ N µ e u(r 12 ,µ) Φ ex 0 [µ](r 1 , r 2 ), (43)
where the normalization factor N µ is simply

N µ = Φ ex 0 [µ]|e 2u(r 12 ,µ) |Φ ex 0 [µ] . (44)
Therefore, in a given basis set B, the ground state right eigenvector

Φ B 0 [µ] of H[µ]
B can be used to estimate the exact ground state wave function Ψ ex 0 through

Ψ B 0 [µ](r 1 , r 2 ) = 1 N µ B e u(r 12 ,µ) Φ B 0 [µ](r 1 , r 2 ), (45)
where the normalization factor N µ B is simply

N µ B = Φ B 0 [µ]|e 2u(r 12 ,µ) |Φ B 0 [µ] . (46)
Then, the quality of a given basis set B to represent H[µ] can be also quantified from the vicinity between Ψ B 0 [µ](r 1 , r 2 ) and the exact wave function Ψ ex 0 (r 1 , r 2 ). To illustrate these ideas, we plotted a cutting of Ψ ex 0 (r 1 , r 2 ) and Ψ B 0 [µ](r 1 , r 2 ) with different basis sets and values of µ. The cutting used is the following: we set two electrons on a circle of radius r = 0.5 a.u. from the nucleus and look at the wave functions as a function of the angle θ 12 between the r 1 and r 2 . For this study, we approximate the exact ground state wave function Ψ ex 0 (r 1 , r 2 ) of the helium atom by the wave function developed by Umrigar et. al. 79 which contains explicitly the r 12 coordinate and which provides a total energy accurate up to twelve digits.

We report in Figs 5 and 6 the plot of the right eigenvector Φ B 0 [µ](r 1 , r 2 ), estimated exact wave function Ψ B 0 [µ](r 1 , r 2 ) for µ = 0.3 and µ = 1.0 (respectively) and compared to the exact wave function Ψ ex 0 (r 1 , r 2 ) and FCI wave function for r = 0.5 a.u..

From these figures, one can notice that: i) the right eigenvector Φ B 0 [µ](r 1 , r 2 ) converges faster than the FCI wave function with respect to the basis set, ii) that for µ = 0.3 the wave function Φ B 0 [µ](r 1 , r 2 ) presents a maximum at r 12 = 0, whereas for µ = 1.0 it has a minimum at coalescence just as the FCI wave function, iii) the wave function Φ B 0 [µ](r 1 , r 2 ) is always larger that the FCI one when r 12 ≈ 0. iv) that Ψ B 0 [µ](r 1 , r 2 ) provides a remarkably good approximation to Ψ ex 0 (r 1 , r 2 ) from the AVTZ basis set.

A qualitative understanding of such observations can be provided by some physical considerations regarding H[µ] and its dependence with µ.

i) The shape of the Jastrow factor u(r 12 ; µ) leading to the transcorrelated Hamiltonian is such that as long as µ < ∞, the effective potential Wee (r 12 , µ) in H[µ] is non divergent, and therefore its eigenfunctions do not have to fulfill the cusp condition which leads to a faster convergence with B. Nevertheless, one can also remark that for µ > 0.5, the speed of convergence deteriorates as µ increases, even though it remains finite. This can be explained by the fact that even if Wee (r 12 , µ) remains bounded from above, when µ is large such operator is poorly represented in a finite basis set.

ii) The fact that the wave function Φ B 0 [µ](r 1 , r 2 ) presents a maximum at r 12 = 0 for µ = 0.3 whereas it provides a minimum for µ = 1.0 can be explained by the shape of the effective potential Wee (r 12 , µ): the latter is much more attractive when r 12 ≈ 0 for µ = 0.3 than for µ = 1.0 (see Fig. 7 for a graphical representation).

iii) Similarly, as for any µ < ∞ one has that Wee (r 12 , µ) < 1 r 12 , the value of the eigenvector of Φ B 0 [µ](r 1 , r 2 ) is necessary larger than the FCI wave function at r 12 = 0 as in the latter the interaction is more repulsive. This implies that the on-top pair density (i.e. the probability of finding two electrons at the same position) is necessary larger for the eigenvectors of H[µ] than that for the eigenvectors of H.

IV. STUDY OF THE HELIUM ISOELECTRONIC SERIES

Having investigated in the case of the helium atom the behaviour of both the eigenvalues and eigenvectors of H[µ] with both µ and the basis set used, we now study the ground state energies of the isoelectronic series of the helium atom from H - to Ne 8+ . As seen from the previous study on the helium atom, the quality of the eigenvalues of H[µ] depends on the value of µ but there is a quite wide regimes of µ for which they converge significantly faster with the basis set with respect to the usual Hamiltonian. Of course, the 'optimal' range of µ might depend on the system and we investigate here different approaches to systematically find a reasonable value of µ and test it on this isoelectronic series which consist in weakly correlated systems We also compare with the results of the TC obtained with the u(r 12 ) = u TN (r 12 ) (see Eq. ( 18)), referred here after as FROGG. covering a quite wide range of densities. As in the study of the helium atom of Sec III, in a given basis set B and for a given system, we use RHF molecular orbitals and the full flexibility is given to the eigenvectors Φ B 0 [µ] . In order to find a systematic way to determine a reasonable value of µ, one must keep in mind the physical meaning of such quantity: it has the unit of the inverse of a distance, and determines, at leading order, the depth of the hole imposed by the Jastrow factor. Therefore, it must have the typical scale of the inverse of correlation effects, which of course strongly depends on the system through its density for instance. In the context of RS-DFT, Toulouse et. al. 80 have investigated different flavour of µ varying in space through the density of the system in a given point in space. More specifically, they introduced (see Eq. ( 12) of 80) a range separation parameter typical for correlation effects in the uniform electron gas (UEG)

µ r s,c (r) = 2 √ α/π √ r s (r) , (47)
with α = (9π/4) -1/3 , r s = (n4π/3) -1/3 and n being the density. Such function µ r s,c (r) depends on the system and on the position in space through the density n(r) of the system at a given point r. Therefore, we propose here to use the average value of µ r s,c (r) over the Hartree Fock electronic density to define the value of µ for a specific system in a given basis set:

µ r s,c = 1 N e drµ r s,c (r)n HF (r) (48)
where n HF (r) is the HF density and N e the number of electrons in the system. We report in Table II the convergence with the basis set B of the error with respect to the exact ground state energies of the isoelectronic series of the helium atom using µ r s,c in H[µ]. We also report in Table II Toulouse et. al. 80 which was proportional 1 r s (r) and found that it leads to values of µ that are too large and therefore do not lead to significant improvement of the basis set convergence.

Still with the aim of finding good values of µ, we propose here the derivation of another approach based on the mapping of the on-top pair density of the UEG and a model built with the Jastrow factor defined by (11). Assuming a single Slater determinant ansatz for a Jastrow Slater wave function with a Jastrow factor defined in (11) and neglecting the corresponding normalization factor, the on-top pair density is

n µ 2 (r) = 1 2 n(r) 2 e -1 √ πµ . (49)
Also, the exact on-top pair density can be estimated from the UEG through

n UEG 2 (r) = n(r) 2 g 0 (n(r)) ( 50
)
where g 0 (n) is the structure factor of the UEG at a given density n. Therefore, one can then find the value µ such that the two on-top pair density coincide

n µ 2 (r) = n UEG 2 (r) ⇔ µ UEG (r) = log 2g 0 (n(r)) √ π . (51)
Then, one can define an average value of µ UEG (r) over the HF density

µ UEG = 1 N e drµ UEG (r)n HF (r). ( 52
)
We report in Table II for different basis sets the error with respect to the exact ground state energy of E B 0 [ µ UEG ] for the isoelectronic series studied here. From these data, one can see that the MAD is sensibly the same than that with the µ r s,c but with a MSD of opposite sign with respect to the latter. This correlates with the fact that the values of µ UEG are in general larger than that of µ r s,c , except for H -. Because the MAD are essentially the same but that the MSD are of opposite sign, it means that there exists an optimal value of µ between µ r s,c and µ UEG which might be optimal. Therefore, we propose to define the average between µ r s,c and µ UEG

µ r UEG s,c = µ r s,c + µ UEG 2 , ( 53
)
and the results obtained are represented in Table II. From this data, one can clearly see that µ r UEG s,c

gives a better MAD and MSD as it is below 1 mH from the double-zeta quality basis set and still improve when enlarging the basis set.

We can also compare the results obtained with µ UEG , µ r s,c and µ r UEG s,c to that obtained with the FROGG and µ FROGG , which are also reported in Table II. From these data one can clearly observe several things: i) while for the quite diffuse systems as H -and He, the energy obtained with the FROGG are of comparable quality that that with µ UEG , µ r s,c and µ r UEG s,c , the situation deteriorates in a monotonic way from Li + to Ne 8+ , ii) from Li + to Ne 8+ , the FROGG converges from below the exact energies, and the larger the Z, the more pronounced is the error with respect to the exact ground state energy, specially in the double-zeta quality basis, iii) the results obtained with µ FROGG show a similar behaviour with respect to that of using the FROGG, iv) the MAD obtained with the FROGG is significantly higher, specially for the double-zeta quality, than that obtained with µ UEG , µ r s,c and µ r UEG s,c . These results illustrates that the Jastrow factor must adapt to the typical length-scale of the system, and therefore shows the benefit of using an easily tunable Jastrow factor such as u(r 12 , µ).

V. THE GROUND STATE POTENTIAL ENERGY CURVE H 2

Having established in Sec II B the analytical form of H[µ] for a general molecular system, we apply here the new TC Hamiltonian H[µ] on the study of the ground state potential curve of H 2 . As in the study of the helium atom of Sec III, in a given basis set B and for a given geometry, we use RHF molecular orbitals and the full flexibility is given to the eigenvectors Φ B 0 [µ] . We report in Figs 8 the difference with respect to the estimated exact ground state potential energy curve using several levels of calculations for the aug-cc-pVDZ and aug-cc-pVTZ basis sets. We do not report the values in the aug-cc-pVQZ because the FCI values are already near the CBS.

From the Figs 8, one can observe that i) all methods provide a near chemical accuracy from the AVTZ basis set, ii) the TC Hamiltonian obtained with the FROGG and that with µ FROGG provide very similar results, iii) the values obtained with µ UEG and µ r s,c provide very similar results and improve accuracy over the FROGG for the AVDZ basis set, iv) a value of µ = 0.5 provides a good description in the AVDZ basis set while remaining above the exact energy in the AVTZ basis set, which is not the case for µ = 0.3. It should also be noticed that when the system becomes strongly correlated (i.e. when the bond is stretched), the TC eigenvalues, whatever the Jastrow chosen, do not show any pathological behaviour.

VI. CONCLUSION

In the present paper, we derived a new form of Jastrow factor u(r 12 , µ) for the TC framework and performed the first numerical tests on a set of two electron atomic and molecular systems. In contrast to the FROGG introduced by Ten-No 55 which was obtained as a linear combination of Gaussians from a leastsquare fit with a fixed weighting function, the analytical form of the new Jastrow factor proposed here has been derived from the analysis of the leading order terms in 1 r 12 of the general TC Hamiltonian, such that it mimics the long-range interaction erf(µr 12 )/r 12 used in the RS-DFT framework (see Sec. II A). The Jastrow factor u(r 12 , µ) resulting from such mathematical conditions naturally imposes the cusp and digs the Coulomb hole over a typical length scale tuned by a unique parameter µ: the smaller the µ, the broader and deeper is the Coulomb hole dug by u(r 12 , µ), and in the µ → ∞ limit, u(r 12 , µ) = 0. The TC Hamiltonian H[µ] obtained from the similarity transformation of the usual Hamiltonian with u(r 12 , µ) is analytical and depends on a unique parameter µ. While varying µ, the physical content of H[µ] significantly changes: one can continuously change from a repulsive non divergent Hamiltonian at large µ, to a partially attractive Hamiltonian at intermediate values of µ, and eventually to a purely attractive Hamiltonian at small values of µ (see Sec. II B). Also, in the µ → ∞ limit one recovers the usual physical Hamiltonian, just as in RS-DFT.

The numerical investigation of the ground state of the helium atom indicates that for values of µ ranging from 0.2 to 1.6 one TABLE II. Error (in mH) with respect to the exact non relativistic energies of the ground state eigenvalue of the usual Hamiltonian and H[µ] for the helium isoelectronic series with Dunning basis sets basis sets for different flavour of µ. For H -and He, the basis sets used are the aug-cc-pVXZ series (X=D,T,Q) and for the Li + -Ne 8+ series, the cc-pCVXZ (X=D,T,Q) basis sets with core-valence functions were used. The mean absolute deviation (MAD) and mean signed deviation (MSD) are also reported for each basis set and method. We also report the average value of the µ considered (referred as µ ) in the double-zeta basis. The method referred as FROGG is the frozen gaussian geminal introduced by Ten-No 55 , and the value of µ FROGG = 0.87 is defined in Eq. (17). The exact non relativistic energy of all systems studied here was obtained from Ref. can strongly improve the convergence of the computed energies with respect to the basis set (see Sec. III). An investigation of the right eigenvectors in real space for different values of µ have shown interesting behaviour: while for µ = 1 the right-eigenvectors provide a depletion of the probability of finding electrons near r 12 = 0, for µ = 0.3 one obtains an increase of such a probability. These different behaviours can be understood as the effective potential entering in the TC Hamiltonian becomes more and more attractive as µ lowers. Nevertheless, when the right eigenvectors are multiplied by the corresponding Jastrow factor, one obtains a very good approximation of the numerically exact ground state wave function (see Sec. III B). Also, we found a value of µ ≈ 0.87 which provides a Jastrow factor very close to the FROGG, and which produces essentially the same eigenvalues than the FROGG, indicating that one can see the FROGG as a special case of the more flexible Jastrow factor u(r 12 , µ).

In order to study the performance of the present approach on a broader set of two-electron systems, we investigated the ground state energy of the helium isoelectronic series from H - to Ne 8+ (see Sec. IV). As the density considerably vary among the isoelectronic series, we propose a very simple scheme to provide a value of µ which automatically adapts to the system, the value of which increases with the nuclear charge of the atom considered. The results obtained for this systemdependent value of µ are found to be very accurate as the MAD with respect to the exact ground state energies over this isoelectronic series in double-zeta quality is about 1 mH and keeps lowering while improving the quality of the basis set. In contrast, the MAD obtained with the FROGG is much higher: it is of about 7 mH in the double-zeta quality basis set, which can be understood by the fact that the typical length scale of the correlation effect imposed by the FROGG is less adapted for core densities than for typical valence densities.

The last part of our study is dedicated to the study of the ground state energy curve of H 2 with different schemes, and we show that the present scheme provides an improvement of the convergence of the energy over the whole potential energy curve.

The perspective of this work are mainly to study more realistic systems which will includes the expensive three-body terms, and to explore the possibility of having a value of µ varying in space.

VII. ACKNOWLEDGMENT

The author would like to thank Julien Toulouse and Andreas Savin for stimulating discussions.

VIII. DATA AVAILABILITY

The data that supports the findings of this study are available within the article and from the corresponding author upon reasonable request.

IX. APPENDIX

A. Analytical form for the operator K(r 1 , r 2 , µ)

To compute the operators K(r 1 , r 2 , µ) one needs to compute first the gradient of u(r 12 , µ) with respect to r 1 which reads

∇ 1 u(r 12 , µ) = 1 -erf(µr 12 ) 2r 12 r 1 -r 2 . (54)
Then, the term ∇ 1 u(r 12 , µ) 2 + ∇ 2 u(r 12 , µ) 2 is simply

∇ 1 u(r 12 , µ) 2 + ∇ 2 u(r 12 , µ) 2 = 1 -erf(µr 12 ) 2 2 . ( 55
)
The non hermitian operators in K(r 1 , r 2 , µ) contains

∇ 1 u(r 12 , µ) • ∇ 1 = 1 -erf(µr 12 ) 2r 12 (x 1 -x 2 ) ∂ ∂x 1 + (y 1 -y 2 ) ∂ ∂y 1 + (z 1 -z 2 ) ∂ ∂z 1 , (56)
and as ∇ 1 u(r 12 , µ) = -∇ 2 u(r 12 , µ), the total non hermitian operator in K(r 1 , r 2 , µ) can be written as

∇ 1 u(r 12 , µ) • ∇ 1 + ∇ 2 u(r 12 , µ) • ∇ 2 = 1 -erf(µr 12 ) 2r 12 (x 1 -x 2 ) ∂ ∂x 1 - ∂ ∂x 2 + (y 1 -y 2 ) ∂ ∂y 1 - ∂ ∂y 2 + (z 1 -z 2 ) ∂ ∂z 1 - ∂ ∂z 2 . (57)
One can notice that the form of Eq. ( 57) provides an explicit form to compute analytical integrals with Gaussian types functions. Nevertheless, one can notice that as

∂ ∂r x 12 = 1 2 ∂ ∂x 1 - ∂ ∂x 2 , (58)
one can write

∇ 1 u(r 12 , µ) • ∇ 1 + ∇ 2 u(r 12 , µ) • ∇ 2 = 1 -erf(µr 12 ) r 12 r 1 -r 2 • ∇ r 12 .
(59) Then, introducing the spherical coordinate system for r 1r 2 in terms of e u = r 1 -r 2 r 12 and the angle θ and φ, one can write

∇ r 12 = ∂ ∂r 12 e u + 1 r 12 ∂ ∂θ e θ + 1 r 12 sin(θ) ∂ ∂φ e φ , (60)
and as r 1r 2 = r 12 e u one obtains

r 1 -r 2 • ∇ r 12 = r 12 ∂ ∂r 12
.

Therefore

∇ 1 u(r 12 , µ) • ∇ 1 + ∇ 2 u(r 12 , µ) • ∇ 2 = 1 -erf(µr 12 ) ∂ ∂r 12 , (62)
which is a more compact expression than Eq. (57).

The computation of the Laplacian reads

∆ 1 u(r 12 , µ) + ∆ 2 u(r 12 , µ) = 2 1 -erf(µr 12 ) r 12 - µ √ π e -µr 12 2 . (63)
B. Computation of integrals involving K(r 1 , r 2 )

The only integrals involved of K(r 1 , r 2 ) over Gaussian functions which are not analytical are of the types v kl i j = dr 1 dr 2 φ i (r 1 )φ j (r 2 ) g(r 12 , µ) 2 φ k (r 1 )φ l (r 2 ), (64) with

g(x, µ) = erf(µx) -1. (65)
To make integrals analytical, we first fit the function erfc(x) with a simple Slater-Gaussian function erfc

(x) ≈ h(x, α, β, c) (66)
with

h(x, α, β, c) = e -αx-βx 2 (67)
and α = 1.09529 and β = 0.756023. Then, by posing y = µx, one obtains

g(x, µ) ≈e -αµx-β(µx) 2 =h(x, αµ, βµ 2 ). (68)
Therefore, one can fit g(x) 2 as

g(x, µ) 2 = 1 -erf(µx) 2 = e -αµx e -βµ 2 x 2 2 = e -2αµx e -2βµ 2 x 2 = h(x, 2αµ, 2βµ 2 ). (69)
Then we fit the Slater function as a linear combination of Gaussians

e -X = N s m=1 c m e -ζ m X 2 . (70)
In the present work, we use N s = 20 and the {c m , ζ m } parameters are reported in Table IX B. By posing X = γx one can fit any Slater function as

e -γx = N s m=1 c m e -ζ m γ 2 x 2 . (71)
Eventually, the function g(x, µ) 2 is obtained as a linear combination of Gaussian

g(x, µ) 2 ≈ N s m=1 c m e -2µ 2 2αζ m +β x 2 , (72)
which makes then the integrals analytical

v kl i j ≈ N s m=1
c m dr 1 dr 2 φ i (r 1 )φ j (r 2 )e -2µ 2 2αζ m +β (r 12 ) 2 φ k (r 1 )φ l (r 2 ). ( 73) All numerical tests performed for µ > 0.1 show that this fit is highly accurate. We need to compute the following integral

φ i φ j φ m | L(r 1 , r 2 , r 3 , µ)|φ k φ l φ n = + φ i φ j φ m | L1 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n + φ i φ j φ m | L2 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n + φ i φ j φ m | L3 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n (74)
where L1 (r 1 , r 2 , r 3 , µ) = w µ (r 12 )w µ (r 13 )r 12 • r 13 (75) with

w µ (r 12 ) = 1 -erf(µr 12 ) 2r 12 . (76)
Let us compute the first term of Eq. ( 74):

φ i φ j φ m | L1 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n =
dr 1 dr 2 dr 3 φ i (r 1 )φ j (r 2 )φ m (r 3 )w µ (r 12 )w µ (r 13 )

r 12 • r 13 φ k (r 1 )φ l (r 2 )φ n (r 3 ). (77)
Such matrix element can be further decomposed into

φ i φ j φ m | L1 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n = φ i φ j φ m | Lx 1 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n + φ i φ j φ m | Ly 1 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n + φ i φ j φ m | Lz 1 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n , (78)
where

φ i φ j φ m | Lx 1 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n = dr 1 dr 2 dr 3 φ i (r 1 )φ j (r 2 )φ m (r 3 )w µ (r 12 )w µ (r 13 ) (x 1 -x 2 )(x 1 -x 3 )φ k (r 1 )φ l (r 2 )φ n (r 3 ). (79)
If we define the following function

W x mn (r) = dr φ m (r )φ n (r )w µ (r -r )(x -x ), (80)
it can be calculated easily by

W x mn (r) = x w mn (r) -w x mn (r) (81)
with

w mn (r) = dr φ m (r )φ n (r )w µ (r -r ), (82)
and

w x mn (r) = dr φ m (r )φ n (r )w µ (r -r )x , (83)
which are analytical integrals. Then, one can rewrite the matrix element of Eq. ( 79) as

φ i φ j φ m | Lx 1 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n = drφ i (r)φ k (r)W x mn (r)W x jl (r), (84)
Therefore, the matrix element of the L1 (r 1 , r 2 , r 3 , µ) is simply

φ i φ j φ m | L1 (r 1 , r 2 , r 3 , µ)|φ k φ l φ n = drφ i (r)φ k (r) W x mn (r)W x jl (r) + W y mn (r)W y jl (r) + W z mn (r)W z jl (r) , (85)
which can be numerically evaluated easily. The total matrix element of the full operator is then

φ i φ j φ m | L(r 1 , r 2 , r 3 , µ)|φ k φ l φ n = drφ i (r)φ k (r) W x mn (r)W x jl (r) + W y mn (r)W y jl (r) + W z mn (r)W z jl (r) + drφ u (r)φ l (r) W x mn (r)W x ik (r) + W y mn (r)W y ik (r) + W z mn (r)W z ik (r) + drφ m (r)φ n (r) W x jl (r)W x ik (r) + W y jl (r)W y ik (r) + W z jl (r)W z ik (r) , (86)
which can be computed on the fly through a simple numerical integration.

D. r 12 → 0 limits for different forms of Hamiltonians and related cusp conditions

Having established the analytical form of the transcorrelated Hamiltonian H[u] for a general Jastrow factor u(r 12 ) in the case of the helium atom in Sec. II A 1, one can then study the behaviour of H[u] when r 12 → 0 and compare it to two other types of Hamiltonian: the usual physical Hamiltonian and that entering the RS-DFT framework. The analysis of the r 12 → 0 behaviour leads to the different condition that the eigenvectors of such operators must fulfill. It should be noticed that, although it was established in the case of the helium atom, the leading order analysis of the r 12 → 0 limit that will be carried out in that section are valid for a general N electron system. For the sake of simplicity of the notations we focus on the ground state of these operators but the arguments are valid for any bounded eigenvectors.

In the case of the usual Hamiltonian, the exact ground state wave function Ψ ex 0 must satisfy the Schroedinger equation in real space HΨ ex 0 (r 1 , r 2 ) = E 0 Ψ ex 0 (r 1 , r 2 ) ∀(r 1 , r 2 ). ( 87)

When looking at r 12 ≈ 0, all terms multiplying 1 r 12 in HΨ ex 0 (r 1 , r 2 ) must remain finite, which, according to Eq. ( 2), translates into lim

Similarly to Eq. (88), when looking at r 12 ≈ 0, the terms multiplying 1 r 12 in H[u]Φ ex 0 (r 1 , r 2 ) must remain finite, which involves also the term W[u] in H[u] (see Eq. ( 7)). Therefore, imposing the finiteness of the limit r 12 → 0 of Eq. ( 91 Φ ex 0 (r 1 , r 2 ) = c, c < ∞.

Then, the term 1 r 12 1-2 ∂u(r 12 ) ∂r 12

in Eq. (92) looks like an effective electron electron interaction induced by the presence of the Jastrow factor in H[u]. Therefore, as long as one imposes that ∂u(r 12 )

∂r 12 r 12 =0 = 1 2 , (93)
which is nothing but the cusp condition for the Jastrow factor u(r 12 ), the scalar term proportional to 1 r 12 in Eq. (92) vanishes at r 12 = 0, and then one obtains a non-divergent effective electron-electron interaction. Within the condition of Eq. (93), a sufficient condition for Φ ex 0 (r 1 , r 2 ) to fulfill the general condition of Eq. ( 92) is ∂Φ ex 0 (r 1 , r 2 ) ∂r 12

r 12 =0 = 0, (94) which implies that as long as the Jastrow factor contains the cusp conditions, the eigenfunction of H[u] is cuspless. For instance, a simple Jastrow factor of the form u(r 12 ) = 1 2 r 12 releases Φ ex 0 (r 1 , r 2 ) from the constraint of fulfilling the cusp condition, and produces a non diverging effective potential for the TC Hamiltonian.

Another form of effective Hamiltonians leading to cuspless eigenfunctions are those obtained from RS-DFT where the Coulomb interaction 1/r 12 is split into a non-divergent longrange interaction erf(µr 12 )/r 12 and a complementary shortrange interaction erfc(µr 12 )/r 12 , where erf(x) and erfc(x) are the error and complementary error functions, respectively. The parameter which controls such a splitting is the so-called range separation parameter µ: RS-DFT reduces the usual WFT when µ → ∞, and µ → 0 it gives back the usual Kohn-Sham theory. In practice, RS-DFT introduces a self-consistent Schroedingerlike equation which must be fulfilled by an effective wave function Ψ µ . Applied to a two electron system and looking in regions where r 12 ≈ 0, the leading terms of the self-consistent Schroedinger-like equation of RS-DFT reads lim erf(µr 12 )

r 12 = 2µ √ π , (96)
one obtains that ∂Ψ µ (r 1 , r 2 ) ∂r 12

r 12 =0 = 0, (97) which means that the wave functions Ψ µ , just as Φ ex 0 are cuspless, because they deal with an effective non-divergent interaction.

Therefore, one can see that there is a similarity between the eigenvectors of RS-DFT and of the transcorrelated Hamiltonian: they are both cuspless as they originate from non divergent Hamiltonians.

E. Large µ limit of H[µ]

As mentioned in Sec. II B 2, one must recover the usual physical Hamiltonian in the large µ limit. Nevertheless, one can notice that lim µ→∞ Wee (r 12 ) =

which seems to contradict the intuitive limit of Eq. ( 32). Nevertheless, as Eq. (100) deals with operators, the equality must be considered in the sense of distributions. Therefore, considering a generic bounded square integrable wave function ψ(r 1 , . . . , r N ) yielding to a integrable pair density n 2 (r 1 , r 2 ), the equality between H[µ] and H is, in the sense of distributions, defined by the following equality ψ|H|ψ = lim where ñ2 (x) is the function n 2 (r 1 , r 2 ) integrated over (r 1 , r 2 ) with the constraint that |r 1r 2 | = x. As the pair density n 2 (r 1 , r 2 ) remains finite and is integrable, the function ñ2 (x) cannot diverge faster than 1

x when x → 0 and therefore one obtains that dr 12 ñ2 (r 12 )δ(r 12 )r 2 12 = 0, (104) which implies Eq. ( 32).

A

  . A new form of Jastrow factor u(µ, r 12 ) tuned by a single parameter µ1. A simple physical picture: the transcorrelated Hamiltonian for the helium atomLet us write the Hamiltonian of the helium atom using the r i = |r i | and r 12 = |r 1r 2 | coordinates76

-1 2 √ 1 2 r 12 exactly 2 √ π r 2

 211222 πµ < 1, the Jastrow factor digs the Coulomb hole, ii) the term in e imposes the cusp; iii) and the last term e -µ 12 reshapes the curvature of the Coulomb hole.

FIG. 1 .

 1 FIG.1. Shape of u(r 12 ; µ) (left) and of exp u(r 12 ; µ) (right) for different values of µ. We also compare the FROGG (see Eq. (18)) and the value of µ = 0.865 which equals the FROGG at r 12 = 0 (see Eq. (17)).

1/r 12 FIG. 3 .

 123 FIG.3. Shape of Wee (r 12 , µ) for µ = µ FROGG and comparison with the effective potential WTN ee (r 12 ) obtained with the FROGG (see Eq. (36)).

FIG. 5 .FIG. 6 .

 56 FIG. 5. Helium atom, radius r = 0.5 a.u.: approximation of the exact wave function Ψ B 0 [µ](r 1 , r 2 ) and right eigenvector Φ B 0 [µ](r 1 , r 2 ) for µ = 0.3 in the AVDZ and AVTZ basis sets (left) and AVQZ and AV5Z (right). Comparison with the FCI wave function in the same basis sets and the estimated exact wave function Ψ ex 0 (r 1 , r 2 ).

W 3 W1/r 12 FIG. 7 .

 3127 FIG.7. Effective potential Wee (r 12 , µ) for µ = 0.3 and µ = 1.0, compared to the usual 1/r 12 Coulomb interaction.

FIG. 8 .

 8 FIG.8. H 2 molecule: Error (in m a.u.) with respect to the estimated exact ground state energy for the aug-cc-pVDZ and aug-cc-pVTZ basis sets (AVDZ, AVTZ, respectively) as a function of the inter nuclear distance for different values of µ, and comparison with FCI. R eq refers to the equilibrium distance of 1.401 a.u., and the estimated exact non relativistic potential energy curve was taken from Ref.81.

12 lim

 12 1 , r 2 ) = c, c < ∞ (88)or equivalently by multiplying Eq. (88) by r similarity transformed Hamiltonians, the exact ground state eigenvector must also satisfy an eigenvalue equation in real-space, H[u]Φ ex 0 (r 1 , r 2 ) = E 0 Φ ex 0 (r 1 , r 2 ) ∀(r 1 , r 2 ).

r 12 Ψ µ (r 1 , r 2 )

 1212 = c, c < ∞,

2 ,

 2 (99) being the Dirac distribution. Threfore, one could be tempted to write that lim µ→∞H[µ] = H + δ(r 12 )H,

  limit for of H[µ] (see Eq. (100)), it implies thatdr 1 dr 2 n 2 (r 1 , r 2 )δ(r 12 ) = 0. (102)To intuitively show Eq. (102), one can perform a series of change of variable (r 1 , r 2 ) → (r 12 , 1 2 (r 1 + r 2 )) and then use the spherical coordinates for r 12 dr 1 dr 2 n 2 (r 1 , r 2 )δ(r 12 ) = dr 12 ñ2 (r 12 )δ(r 12 )

The working equation for u(µ, r 12 ) inspired by RSDFT

	want now to impose the form of u(r 12 ) such that it mimics,
	at leading order in 1/r 12 , the long-range effective interaction
	entering in the effective Hamiltonian of the RS-DFT (see Eq.
	(95) in Appendix IX D). Mathematically, this condition implies
	that
			W[u] +	1 r 12	=	erf(µr 12 ) r 12
	⇔ -	2 r 12	∂u(r 12 , µ) ∂r 12	+	1 r 12	=	erf(µr 12 ) r 12	,
								1 r 12 + W[u] + w[u], and an additional differential operator t[u], the latter
								making H[u] non-hermitian.
								Having established the form of H[u] in Eq. (8), one can
								notice that the leading order terms in 1/r 12 are 1 r 12 + W[u]. We

2.

  Convergence of the ground state eigenvalue of H[µ] with respect to the AVXZ basis set series (X=D,T,Q,5) for different values of µ.

		-2.888				-2.888
	Energy (a.u.)	-2.9 -2.898 -2.896 -2.89 -2.892 -2.894	µ = 0.2 µ = 0.3 µ = 0.5 µ = 0.7 FROGG FCI Exact		Energy (a.u.)	-2.89 -2.892 -2.898 -2.896 -2.894	µ = 0.7 FROGG µ = 2.0 µ = 3.0 FCI Exact
		-2.904 -2.902				-2.9
		-2.906				-2.902
		-2.908 AV5Z	AV4Z	AV3Z	AV2Z	-2.904 AV5Z	AV4Z	AV3Z	AV2Z
				1/X				1/X
	FIG. 4.

TABLE I .

 I Ground state eigenvalue (in Hartree) of H[µ] for the He atom with the AVXZ (X=D,T,Q,5) basis sets for different values of µ, and error (in mH) with respect to the exact non relativistic energy. We also compare with the results of the TC obtained with the u(r 12 ) = u TN (r 12 ) (see Eq. (18) and Ref.55), referred here as FROGG, and the value of µ FROGG = 0.87 is defined in Eq. (17). The exact non relativistic energy was obtained from Ref.78.

	Basis/µ	0.2	0.3	0.35	0.5	0.7	µ FROGG = 0.87
	AVDZ	-2.906309 (-2.58)	-2.904172 (-0.45)	-2.903317 (0.41)	-2.901420 (2.3)	-2.899516 (4.21)	-2.897927 ( 5.80)
	AVTZ	-2.905278 (-1.55)	-2.904691 (-0.97)	-2.904468 (-0.74)	-2.903969 (-0.24)	-2.903488 (0.24)	-2.903214 ( 0.51)
	AVQZ	-2.904325 (-0.6)	-2.903986 (-0.26)	-2.903879 (-0.15)	-2.903702 (0.02)	-2.903620 (0.1)	-2.903589 ( 0.13)
	AV5Z	-2.904229 (-0.51)	-2.903928 (-0.2)	-2.903834 (-0.11)	-2.903702 (0.02)	-2.903676 (0.05)	-2.903678 (-0.1)
	Basis/µ	1.0	1.6	2.0	3.0	∞ (FCI)	FROGG
	AVDZ	-2.896734 (6.99)	-2.892773 (10.95)	-2.891477 (12.25)	-2.890212 (13.51)	-2.889548 (14.18)	-2.897972 (5.75 )
	AVTZ	-2.903069 (0.66)	-2.902421 (1.3)	-2.901932 (1.79)	-2.901159 (2.56)	-2.900598 (3.13)	-2.903919 (-0.20)
	AVQZ	-2.903558 (0.17)	-2.903359 (0.36)	-2.903244 (0.48)	-2.902951 (0.77)	-2.902534 (1.19)	-2.903833 (-0.11)
	AV5Z	-2.903675 (0.05)	-2.903634 (0.09)	-2.903592 (0.13)	-2.903479 (0.24)	-2.903201 (0.52)	-2.903825 (-0.10)
				Exact non relativistic: -2.903724

  the value of µ r s,c in the double-zeta quality basis for each system, the value in larger basis sets varying by less than 0.1%. As one can observe from TableII, the ground state energies of H[ µ r s,c ] converge much faster than that of the usual Hamiltonian as the MAD for the ground state energies is of 2.64 mH, 0.48 mH and 0.26 mH in the double-, triple-and quadruple zeta quality basis sets (respectively) whereas it is of 12.16 mH, 4.71 mH and 1.64 mH for the usual Hamiltonian. Also, one can observe that, except for the low density systems such as the H -and He species, it converges from below the exact ground state energy. Regarding now the value of µ r

		0.34
		0.32
	Wave functions, µ = 0.3	0.22 0.24 0.26 0.28 0.3
		0.2
		0.18	-4	-3	-2	-1	0	1	2	3	4
							θ 12

s,c , one can notice that it increases with the atomic number. This is expected as the density becomes more picked near the nucleus while the nuclear charge increases, which induces necessary that the typical inter electronic distance lowers with Z. In a similar fashion, we tried another definition of µ based on the work of

  78.

				H -						He					Li +
		DZ	TZ			QZ	µ	DZ	TZ		QZ	µ	DZ	TZ		QZ	µ
	FCI	3.72	1.19			0.61	+∞	14.18	3.13		1.19	+∞	10.72	3.35		1.58	+∞
	E B 0 [ µ UEG ] E B 0 [ µ rs,c ] E B 0 [ µ r UEG ] s,c E B 0 [FROGG] E B 0 [µ FROGG ]	1.69 2.30 2.03 3.22 3.12	0.60 0.67 0.64 0.91 0.91	Be 2+	0.38 0.39 0.39 0.45 0.46	0.350 0.479 0.410 --	5.28 4.87 5.07 5.75 5.80	0.44 0.37 0.40 -0.20 0.51	B 3+	0.12 0.12 0.12 -0.11 0.13	0.815 0.771 0.792 --	0.9 -0.82 0.04 -4.08 -1.52	0.09 -0.38 -0.12 -1.24 -0.61	C 4+	-0.02 -0.14 -0.08 -0.56 -0.21	1.274 0.980 1.127 --
		DZ	TZ			QZ	µ	DZ	TZ		QZ	µ	DZ	TZ		QZ	µ
	FCI	11.35	3.62			1.41	+∞	11.91	4.21		1.53	+∞	12.46	4.76		1.67	+∞
	E B 0 [ µ UEG ] E B 0 [ µ rs,c ] E B 0 [ µ r UEG ] s,c E B 0 [FROGG] E B 0 [µ FROGG ]	1.2 -1.5 -0.08 -6.38 -3.81	0.17 -0.53 -0.15 -1.63 -1.02	N 5+	-0.02 -0.17 -0.09 -0.53 -0.32	1.727 2.152 1.440 --	1.41 -2.16 -0.18 -8.24 -6.33	0.43 -0.34 0.09 -1.45 -1.05	O 6+	-0.01 -0.22 -0.11 -0.49 -0.45	2.179 1.300 1.740 --	1.67 -2.69 -0.19 -9.95 -8.48	0.37 -0.4 0.01 -1.44 -1.20	F 7+	0.01 -0.25 -0.1 -0.48 -0.58	2.631 1.434 2.032 --
		DZ	TZ			QZ	µ	DZ	TZ		QZ	µ	DZ	TZ		QZ	µ
	FCI	13.1	5.71			1.79	+∞	13.84	6.47		2.02	+∞	14.68	7.06		2.25	+∞
	E B 0 [ µ UEG ] E B 0 [ µ rs,c ] E B 0 [ µ r UEG ] s,c E B 0 [FROGG] E B 0 [µ FROGG ]	2.06 -2.97 -0.01 -11.35 -9.98	0.45 -0.64 -0.15 -1.74 -1.54	0.01 -0.26 -0.09 -0.39 -0.68 Ne 8+	3.082 1.556 2.318 --	2.54 -3.13 0.27 -12.60 -11.09	0.77 -0.70 -0.04 -1.81 -1.73 MAD -0.02 -0.31 -0.12 -0.34 -0.80	3.533 1.670 2.600 --	3.14 -3.10 0.69 -13.45 -11.69	1.15 -0.54 0.23 -1.55 -1.63 MSD -0.02 -0.36 -0.14 -0.32 -0.90	3.984 1.774 2.879 --
		DZ	TZ			QZ	µ	DZ	TZ		QZ		DZ	TZ		QZ
	FCI	15.66	7.61			2.36	+∞	12.16	4.71		1.64		12.16	4.71		1.64
	E B 0 [ µ UEG ] E B 0 [ µ rs,c ] E B 0 [ µ r UEG ] s,c E B 0 [FROGG] E B 0 [µ FROGG ]	3.89 -2.85 1.26 -13.79 -11.79	1.57 -0.28 0.57 -1.17, -1.41,		0.03 -0.35 -0.1 -0.26 -0.90	4.434 1.874 3.1543 --	2.38 2.64 0.98 8.88 7.36	0.60 0.48 0.24 1.31 1.16		0.07 0.26 0.13 0.39 0.54	--	2.38 -1.21 0.89 -7.09 -5.58	0.60 -0.28 0.15 -1.13 -0.88		0.07 -0.16 -0.03 -0.30 -0.42	--

TABLE III .

 III Set of coefficients c m and exponents ζ m for the fit of e -X

	ζ m	c m
	30573.77073	0.00338925525
	5608.452381	0.00536433869
	1570.956734	0.00818702846
	541.3978511	0.01202047655
	212.4346963	0.01711289568
	91.31444574	0.02376001022
	42.04087246	0.03229121736
	20.43200443	0.04303646818
	10.37775161	0.05624657578
	5.468807545	0.07192311571
	2.973735292	0.08949389001
	1.661441902	0.10727599240
	0.9505256082	0.12178961750
	0.5552868397	0.12740141870
	0.3304336002	0.11759168160
	0.1998230323	0.08953504394
	0.1224684076	0.05066721317
	0.07575825322	0.01806363869
	0.04690146243	0.00305632563
	0.02834749861	0.00013317513

C. Analytical form of L(r 1 , r 2 , r 3 , µ) and computation of related integrals


J Deng

B Camenen

T Drevet

L Pénard

Evaluation of the deposit dynamics on a gravel bar after different hydrologic events

Estimating the fine sediment budget in alpine rivers is an important challenge to better understand the morphodynamics of such rivers. In this paper, we describe the dynamics of fine-sediment deposits over a gravel bar before and after two flushing events in an alpine river based on pictures taken from the side of the river. These pictures are processed using a new image analysis method. It allows to estimate the spatial distribution of fine-sediment deposit on the gravel bar using photos taken at different dates. Photos before and after the flushing events in 2018 and 2019 were analysed. Comparison with the local in-situ measurement of surface grain size distribution has been carried out using the Wolman peddle count method. The image detection method presents a good agreement on the spatial distribution and variations of fine sediment deposits on the gravel bar with an acceptable bias. A discussion is provided on the limits and validation of the proposed method as well as on its applicability on medium term analysis.

INTRODUCTION

A significant amount of fine sediments can be found in piedmont gravel bed rivers. Gravel bars represent one of the main zones of storage [START_REF] Bridge | Description and interpretation of fluvial deposits: a critical perspective[END_REF]. It remains however very difficult to estimate the total amount of fine sediment stock on gravel bars since these sediments are easily washed out and can also infiltrate into the gravel matrix. Indeed, it often represents a nonnegligible part of the sediment flux passing through the system [START_REF] Misset | The dynamics of suspended sediment in a typical alpine alluvial river reach: Insight from a seasonal survey[END_REF]. Moreover, very different dynamics can be observed from one class of sediments to another, i.e. washload of silts versus graded-suspension of sands [START_REF] Camenen | Estimation of the volume of a fine sediment deposit over a gravel bar during a flushing event[END_REF]. It is thus a real challenge to provide fine sediment budgets over a system. The purpose of this paper is to present some field experiments describing the dynamics of fine-sediment deposits during two flushing events in an alpine river based on pictures taken from the side of the river [START_REF] Camenen | Estimate of fine sediment deposit dynamics over a gravel bar using photography analysis[END_REF]. Using a new image analysis method based on clustering (presented in section 3, [START_REF] Lloyd | Least squares quantization in pcm[END_REF] , this paper provides an analysis of fine sediment dynamics (budget and spatial distribution) over a gravel bar after two dam flushing events in 2018 and 2019 (presented in section 4). A discussion of the results is then presented using local measurements of surface grain size distribution before and after a flushing event in 2019.

FIELD MEASUREMENTS

Location of the field measurements

The Arc en Maurienne River is an engineered alpine river characterised by a nival hydrologic regime and marked by an intense input of fine sediments from the catchment. The total area of the catchment is 1957 km 2 at the confluence with the Isère River. The area of interest is a single gravel bar located within a system of alternate bars 9 km downstream of St-Jean-de-Maurienne (see Figure 1). This reach, which has a mean slope of 0.6%, has been studied since 2005 [START_REF] Jaballah | Alternate bar development in an alpine river following engineering works[END_REF]. One point of interest of this site is that flushing events of the three run-of-the-river dams are conducted yearly in June allowing the performance of in-situ experiments at a specified date known in advance. Topographic, bathymetric and flow measurements are regularly carried out during these flushing events. A focus has been made on the fine sediments dynamics during the flushing events [START_REF] Antoine | Estimation des flux de matières en suspension lors des chasses hydrauliques de l'Arc de 2006 à 2011[END_REF][START_REF] Antoine | Dynamique des matériaux en suspension le long de rivières aménagées de montagne[END_REF][START_REF] Camenen | Estimation of the volume of a fine sediment deposit over a gravel bar during a flushing event[END_REF][START_REF] Camenen | Estimation of sand suspension in a secondary channel of an alpine river[END_REF] and the impact of the gravel bar on fine sediment exchange between the bed and the flow appeared to be a significant issue [START_REF] Camenen | Field experiment on the dynamics of fine sediments over a gravel bar in an alpine gravel-bed river[END_REF].

Hydro-sedimentary measurements

Generally, a flushing event lasts around 12 hours and follows a designed hydrograph with two plateaus, one at Q ≈ 80 m 3 /s and a second one at Q ≈ 150 m 3 /s, both approximately lasting 4 hours. The studied gravel bar is totally inundated for a flow discharge Q ≈ 100 m 3 /s. During both flushing events in 2018 and 2019, the inundation was clearly photographed by the camera implemented on the side of the river (see Figure 2(b)).

The field campaigns included topographic measurements of the gravel bar before and after the event. A hydro-sedimentary station located 9 km upstream the study site (at St-Jean-de-Maurienne) monitored the fine sediment fluxes. During the previous flushing events, very large concentrations of SSM (Suspended Sediment Matter) were often observed, up to C max = 30 g/l [START_REF] Antoine | Dynamique des matériaux en suspension le long de rivières aménagées de montagne[END_REF]. In 2018, a peak of concentration of 32 g/l has been observed during the flush event (see Figure 2(c)). However, in 2019, a high concentration of SSM has been observed 4.5 hours before the flushing event with a peak of 31.5 g/l, and the peak concentration of the flushing event reached only C max = 8.7 g/l (see Figure 2(d)).

Local sediment measurement

Wolman pebble count method [START_REF] Wolman | A method of sampling coarse river-bed material[END_REF] was carried out to measure the surface grain size distribution on the studied gravel bar before and after the flushing event in 2019. Following the instructions of the Wolman method, 5 patches (see Figure 3(a)) were firstly decided. An operator went all over each patch along several backwards and forwards stochastic trajectories, picked up the particle beneath the tip of the toe of the boot and measured the intermediate axis b of this particle by the use of gabarit from the FISP (Federal Interagency Sedimentation Project) US SAH97 model.

A grain size distribution (GSD) analysis from Wolman method samples has been carried out for the flushing event in 2019. In Figure 3(b), the cumulative frequency of GSD for the samples were plotted. Looking at the fine sediment fraction (d<10 mm), the frequency GSD shows a decrease of frequency for patch 1 and an increase of frequency for patch 2. This may be due to the fact that the higher elevation of patch 2 causes a low water depth and low velocities which results in a low bed shear stress and a fine sediment deposition. On the contrary, some fine sediment erosion happened in patch 1 because it is close to the secondary channel where local bed shear stress is relatively higher (exceed 10 Pa for the 2015 flush according to, [START_REF] Camenen | Estimation of the volume of a fine sediment deposit over a gravel bar during a flushing event[END_REF].

IMAGE DETECTION METHOD

Data collection and pre-processing

A new image detection method developed for the estimation of the surface deposition of fine sediment on a gravel bed from an image was carried out. A preliminary version was presented by [START_REF] Pénard | Image-based fine sediment detection on gravel bars surface[END_REF] Images were firstly acquired by a time-lapse camera attached on a tree trunk on the top of river bank. They were collected automatically with a frequency of 1 image per 30 minutes during the flushing event and 2 images per day on the other days. The resolution of the camera is approximately 5 cm for each pixel. The camera was calibrated by fixing the intrinsic/extrinsic parameters. Some field measurements were then carried out including Ground Reference Points (GRPs) and image corresponding points selection to prepare for the image orthorectification. Finally images were orthorectified according to the calibration parameters of camera, GRPs and a Digital Elevation Model (DEM).

Bar segmentation

Once data had been collected and pre-processed, the bar segmentation was conducted to estimate the surface area of fine sediment deposits. The following method is illustrated by Figure 4.

The first step of bar segmentation consisted in sharpening the image to emphasize fine details in the image and to get the pixels representing the fine sediments boosted. Starting from the orthorectified image, the image with RGB (Red, Green, Blue) color space was firstly converted into the HSV (Hue, Saturation, Value) color space. Secondly the Fast Fourier Transform (FFT) was applied to obtain the frequency of the image. Then a high-pass filter on the frequency domain of image was implemented to emphasize the pixels of the fine sediment deposits. Finally the inverted FFT was done on the filtered image and the image was switched back to the HSV image domain.

The second step consisted in segmentation of the image into several clusters minimizing withincluster variances to classify the pixels to represent the deposits of fine sediment. In this step, an algorithm called K-means clustering [START_REF] Lloyd | Least squares quantization in pcm[END_REF] was used. Firstly, user decided the number of clusters K. Secondly, the K-means algorithm created random K values representing the centroid value of each cluster and assigned each pixel to the closest centroid after that the pixels form K clusters. Next, the algorithm computed the variance of each cluster and replace the centroid of each cluster until the algorithm minimized their variance. Finally, the pixels in the first cluster were selected as the fine sediment deposit region. It is noted that the algorithm classified the pixels only taking into account the value of each pixel, whereas the link between the value of pixels in the image and the actual size of sediment was unclear due to the limitation of the image resolution. So this method might lead to a systematic bias when calculating the area of deposits.

The last step was the morphological opening and closing of the resulting region. This step aimed at reshaping the selected region by applying morphological dilation and erosion to have a pertinent region.

Validation with Wolman pebble count method

The Wolman pebble count method was considered as a reference method to validate the reliability of the image detection method. Comparison between the results of these two methods has been conducted for the flush in 2019. It is noted that the Wolman pebble count method provided a volume GSD whereas image detection method provided an areal GSD. To be able to compare results of the two methods, firstly, the Wolman GSD was transformed into the areal GSD by the following equation: where Φ i W olman = volumetric percentage of class i in Wolman; D i = diameter of class i; and S i = areal percentage of class i. Secondly, the areal percentage estimated by the image detection method in the Wolman patches was calculated. Finally, the comparison was done for the class

Φ i Area = S i S = Φ i W olman D i S (1)
D i ≤ 2 mm.
Results are plotted in Figure 5. The spatial variation of fine sediment deposits seems to be successfully produced by both two methods. However, the area percentages of fine sediment deposits for all patches are underestimated by the image detection method. It may be due to the selection of the first class of K-means cluster only, that resulted in neglecting the fine sediment deposits from other classes. In addition, the limitation of the pixel size could cause the lack of representation of fine sediment deposits close to the large particles (pebbles or boulders). The pixels of the fine sediment deposition around the pixels of the coarse particles could also be eliminated by the dilation and erosion of mathematical morphological analyse during the image processing. Looking at the percentage from different photos, a non-negligible difference was also noted. This may be due to the time when the photo had been taken causing a difference on the brightness of the photos. This difference on the brightness may lead to a bias when K-means clustering. The pixels would be classified in the same cluster because of the brightness but they actually represented different sediments. A future work could be done on the sensitivity of the light and weather conditions when the photo is taken. Concerning Wolman pebble count method, even if it was a reference method, it produced a relatively high statistic bias which led to a statistical percent error (±20% for a sample size of 400 particles) for the estimation of D 5 [START_REF] Bunte | Sampling surface and subsurface, particle-size distributions in wadable gravel-and cobble-bed streams for analyses in sediment transport, hydraulics, and streambed monitoring[END_REF], i.e. a significant error on the estimation of fine sediment content. Thus, the bias observed in Figure 5 could either come from the Wolman pebble count method or from image detection method.

ESTIMATION OF FINE SEDIMENT DEPOSIT FOR TWO FLUSHING EVENTS

The image detection method has been implemented to evaluate the deposit dynamics for the 2018 and 2019 flushing events. One photo before and after flush were selected for each flush. The pixel size for both orthorectified photos is 5 cm and the bar has been photographed under a similar light condition to minimize the brightness bias. In Figure 6, the direction of the river flow is from the bottom right of the photo to the upper left, the main channel of the river is on the upper right of the photo and the secondary channel is located on the bottom left of the photo.

The surface areas of deposits are presented in Table 1 for the four studied photos. The bar has lost 5.3% of fine sediment surface area after the 2018 flush and 6.1% after the 2019 flush, respectively. Figure 6 shows the corresponding maps of distribution of the deposits. Before the 2018 flush, fine sediment deposits were mostly distributed at the tail of the bar and slightly distributed in the zone close to the secondary channel. After the flush, fine sediments at these two positions were mostly washed out and new deposits were found in the middle of the bar. A similar result could be found for the 2019 flush.

To explain this phenomenon, one should mention the existence of a transverse channel delivering water from secondary channel to main channel at the tail of the bar. During the flood, the water depth in the transverse channel grew with the increase of the discharge; high bed shear stresses can thus be observed leading to some erosion in this transverse channel [START_REF] Camenen | Estimation of the volume of a fine sediment deposit over a gravel bar during a flushing event[END_REF]. However, in the middle of the bar, due to lower bed shear stresses (low water depth), fine sediment deposits could not be resuspended. For other zones on the bar, it was difficult to explain the dynamic of fine sediment deposits by the bed shear stress. Indeed, a phase-lag between the flushing flood wave and the velocity peak propagation was found in the previous flushing events. This phase-lag caused a complex combination of water depth and velocity distribution, which led to a variation in time of local bed shear stress on the bar and thus complicated the dynamic of fine sediment deposits. Furthermore, a difference between the deposition after the 2018 flush and the deposition before the 2019 flush can also be observed. This could be due to the fact that a natural flood has been taken place on 23 October 2018 with a peak discharge of 163 m 3 /s and a peak concentration of 14.7 g/l.

CONCLUSION

A new image detection method has been validated using local field measurements and implemented on a gravel bar to evaluate the fine sediment deposit dynamic after two hydrological flushing events. A validation with Wolman pebble count method showed that image detection method produced a good agreement on the estimation of spatial distribution of fine deposits on the bar. However, an averaged 40% underestimation was found. This underestimation was, from the side of the image detection method, due to the bias of the image processing and the limitation of pixel size, and from the side of the Wolman pebble count method, due to a relatively high statistical percent error when estimating the fine sediment content. A variation of the results from different photos taken a priori while the surface deposit distribution was the same indicated that the method might be sensitive to the brightness of the input photos. The image detection method has been used to evaluate the dynamics of fine sediment after two flushing events (2018 and 2019). The results showed that some erosion took place at the tail of the bar and on the right side of the secondary channel for both two flushing events. As observed by [START_REF] Camenen | Estimation of the volume of a fine sediment deposit over a gravel bar during a flushing event[END_REF], the phase-lag between water depth and velocities on the bar makes the understanding of the fine sediment dynamics complicated to explain by just looking at the local bed shear stresses. One perspective will be the implementation of the image detection method on the photo time series (during several years) to be able to compare the dynamics of fine sediments for different hydrological events and also to link the dynamic of fine sediment with other factors such as rainfall intensity and wind, which could have an impact on the fine sediment deposits on the bar.

Figure 1 .

 1 Figure 1. Arc Catchment and location of experimental site (green points correspond to run-of-the-river dams)

Figure 2 .

 2 Figure 2. Photo of the bar before the flushing event in 2019 under a discharge of 47.4 m 3 /s (a), photo of the inundation over the bar during the flushing event (under a discharge of 129 m 3 /s) (b) and hydrographs and sediment-graphs measured at St-Jean-de-Maurienne hydro-sedimentary station during the flushing event in 2018 (c) and 2019 (d)

Figure 3 .

 3 Figure 3. Location of local sediment measurement (Wolman) patches (a) and Granulometric curves from Wolman peddle count method (b).

Figure 4 .Figure 5 .

 45 Figure 4. The first step: High-pass filter on the Fourier transform of the Value channel (a), the second step: K-means clustering: segmentation into k clusters minimizing within-cluster variances (b) and the last step: Morphological opening and closing of the resulting regions (c)

Figure 6 .

 6 Figure 6. The map of deposit before (a) and after (b) the flush in 2018 and the map of deposit before (c) and after (d) the flush in 2019 estimated by the image detection method.

Table 1 .

 1 Results of estimation of fine sediment deposition area for the flush 2018 and 2019.

	Year	Photo	S bar (m 2 )	S Dep (m 2 )	S Dep /S bar (%)
	2018	before	575.5	252,0	43.8
		after	581.7	223.8	38.5
	2019	before	594.5	301.1	50.6
		after	587.7	261.4	44.5

ACKNOWLEDGEMENTS

This study was supported by Irstea and EDF-CIH (Electricité De France -Centre d'Ingénierie Hydraulique) and by the French National Research Agency (ANR) under the grant ANR-18-CE01-0020 DEAR Project.