Julien Paquier
email: julien.paquier@lct.jussieu.fr

Emmanuel Giner

Julien Toulouse
email: toulouse@lct.jussieu.fr

Relativistic short-range exchange energy functionals beyond the local-density approximation

We develop relativistic short-range exchange energy functionals for four-component relativistic rangeseparated density-functional theory using a Dirac-Coulomb Hamiltonian in the no-pair approximation. We show how to improve the short-range local-density approximation exchange functional for large range-separation parameters by using the on-top exchange pair density as a new variable. We also develop a relativistic short-range generalized-gradient approximation exchange functional which further increases the accuracy for small rangeseparation parameters. Tests on the helium, beryllium, neon, and argon isoelectronic series up to high nuclear charges show that this latter functional gives exchange energies with a maximal relative percentage error of 3%. The development of this exchange functional represents a step forward for the application of four-component relativistic range-separated density-functional theory to chemical compounds with heavy elements.

I. INTRODUCTION

Range-separated density-functional theory (RS-DFT) (see, e.g., Refs. 1 and 2) is an extension of Kohn-Sham densityfunctional theory (DFT) [3] which rigorously combines a wave-function method accounting for the long-range part of the electron-electron interaction with a complementary shortrange density functional. RS-DFT has a faster basis convergence than standard wave-function theory (WFT) [4] and can improve over usual Kohn-Sham density-functional approximations (DFAs) for the description of strong-correlation effects (see, e.g., Refs. 5 and 6) or weak intermolecular interactions (see, e.g., Refs. 7 and 8).

For the description of compounds with heavy elements, RS-DFT can be extended to a four-component relativistic framework [9][10][11]. In particular, in Refs. 9 and 10, second-order Møller-Plesset perturbation theory and coupled-cluster theory based on a no-pair Dirac-Coulomb Hamiltonian with longrange electron-electron interaction were combined with shortrange non-relativistic exchange-correlation DFAs and applied to heavy rare-gas dimers. One limitation, at least in principle, in these works is the neglect of relativity in the short-range density functionals. It is thus desirable to develop appropriate short-range relativistic exchange-correlation DFAs for fourcomponent RS-DFT in order to quantify the error due to the neglect of relativity in the short-range density functionals and possibly increase the accuracy of these approaches. As a first step toward this, in Ref. 11 some of the present authors developed a short-range relativistic local density-functional approximation (srRLDA) exchange functional based on calculations on the relativistic homogeneous electron gas (RHEG) with the Coulomb and Coulomb-Breit electron-electron interactions in the no-pair approximation.

In the present work we test this srRLDA exchange functional on atomic systems, namely the helium, beryllium, neon, and argon isoelectronic series up to high nuclear charges Z, using a four-component Dirac-Coulomb Hamiltonian in the no-pair approximation. We reveal that, for these relativistic ions with large Z, the srRLDA exchange functional is quite inaccurate even for large values of the range-separation parameter µ. We show how this functional can be improved by using the on-top exchange pair density as a new variable. Finally, we further improve the short-range relativistic exchange functional by constructing a generalized-gradient approximation (GGA), achieving a 3% maximal relative energy error.

The paper is organized as follows. In Section II we lay out the formalism of RS-DFT for a four-component relativistic Dirac-Coulomb Hamiltonian in the no-pair approximation. In Section III we give the computational details for our calculations. In Section IV we test the srRLDA exchange functional and discuss its limitations. In Section V we improve the sr-RLDA exchange functional by using the on-top exchange pair density. In Section VI, we construct and test short-range relativistic exchange GGAs. Finally, Section VII contains our conclusions. In the Appendices, we derive the uniform coordinate scaling relation for the relativistic short-range exchange density functional and the expression of the on-top exchange pair-density in a four-component no-pair relativistic framework.

II. RELATIVISTIC RANGE-SEPARATED DENSITY-FUNCTIONAL THEORY

Let us first discuss the choice of the relativistic quantum many-particle theory on which to base relativistic RS-DFT and the general strategy that we follow in this work. Clearly, since RS-DFT combines WFT and DFT, we need a relativistic framework which is convenient for both of them. Relativistic Kohn-Sham DFT has been formulated based on quantum electrodynamics (QED) [12][START_REF] Engel | Density Functional Theory[END_REF][START_REF] Engel | Relativistic Electronic Structure Theory, Part 1: Fundamentals, Theoretical and Computational Chemistry[END_REF], even though the no-pair approximation [START_REF] Sucher | [END_REF]16] is normally introduced at a later stage for practical calculations. As regards WFT, the best tractable relativistic framework is the recently developed effective QED Hamiltonian [17][18][19][20][21] incorporating all QED effects obtained with non-retarded two-particle interactions. In the present work, we will stick however to the most common choice of the four-component Dirac-Coulomb Hamiltonian in the nopair approximation, which can easily be used for both WFT and DFT. This relativistic framework can be derived in several ways (see, e.g., Refs. 17,18,[20][21][22].

In the effective QED approach, the Hamiltonian is written in second quantization with normal ordering with respect to the vacuum state of empty positive-energy oneparticle states and completely filled negative-energy oneparticle states, while incorporating charge-conjugation symmetry. This Hamiltonian has a stable vacuum state and is physically meaningful. In this approach, the no-pair approximation, corresponding to projecting this Hamiltonian onto the many-electron subspace generated by positive-energy oneparticle states, is just a convenient approximation (but in principle not necessary), akin to the idea of restricting the orbitals entering the wave function to an active orbital subspace in the complete-active-space self-consistent field method. The no-pair Dirac-Coulomb Hamiltonian is then obtained by further neglecting the effective QED one-particle potential corresponding to vacuum polarization and electron self-energy. By contrast, in the configuration-space approach, the Hamiltonian is written either in first quantization or, equivalently, in "naive" second quantization (i.e., without normal ordering with respect to a stable vacuum state). The resulting Hamiltonian has thus an unstable vacuum state, corresponding to empty positive-energy one-particle states and empty negativeenergy one-particle states, it has no bound states (the electronic states that should be normally bound being embedded in the continuum of excitations to positive-energy states and deexcitations to negative-energy states), and hence is per se unphysical. However, by projecting this Hamiltonian onto the many-electron subspace generated by the positive-energy oneparticle states, we recover the same physically relevant no-pair Dirac-Coulomb Hamiltonian as the one obtained by starting with the effective QED approach.

One drawback of the no-pair approximation is that the projector onto the subspace of electronic states depends on the separation between positive-energy and negative-energy oneparticle states, and therefore depends on the potential used to generate these one-particle states. If the projector is applied to the Hamiltonian, the whole resulting projected Hamiltonian is thus dependent on this potential. As mentioned in Ref. 14, this is problematic for formulating DFT since we cannot isolate, as normally done, an universal part of the Hamiltonian, and we thus cannot define universal density functionals. However, instead of thinking of the projector as being applied to the Hamiltonian, we can equivalently think of the projector as being applied to the considered many-electron state and optimize the projector simultaneously with the wave function. In this way, we can introduce universal density functionals, similarly to non-relativistic DFT, defined such that for a given density a constrained-search optimization of the projected wave function will determine alone the optimal projector without the need of pre-choosing a particular potential, at least for systems for which positive-energy and negativeenergy one-particle states can be unambiguously separated. Again, both the effective QED approach or the configurationspace approach can a priori be used for doing so. In the effective QED approach, the projector would be optimized (by rotations between positive-energy and negative-energy one-particle states) using an energy minimization principle. In the configuration-space approach, the projector is optimized using a minmax principle [23][24][25][26][START_REF] Saue | Theoretical Chemistry and Physics ofHeavy and Superheavy Elements[END_REF][START_REF] Almoukhalalati | [END_REF]. In the present work, we follow the configuration-space approach and leave for future work the alternative formulation based on the effective QED approach.

We thus start with the Dirac-Coulomb electronic Hamiltonian (see, e.g., Refs. 29 and 30)

Ĥ = TD + Vne + Ŵee , (1)
where the kinetic + rest mass Dirac operator TD , the nucleielectron interaction operator Vne , and the Coulomb electronelectron interaction operator Ŵee are expressed using fourcomponent creation and annihilation field operators ψ † (r) and ψ(r) without normal reordering with respect to a stable vacuum state. We thus write TD as

TD = ψ † (r) c (α • p) + β mc 2 ψ(r) dr, (2)
where p = -i∇ r is the momentum operator, c = 137.036 a.u. is the speed of light, m = 1 a.u. is the electron mass, and α and β are the 4 × 4 Dirac matrices

α = 0 2 σ σ 0 2 and β = I 2 0 2 0 2 -I 2 , (3)
where σ = (σ x , σ y , σ z ) is the 3-dimensional vector of the 2×2 Pauli matrices, and 0 2 and I 2 are the 2 × 2 zero and identity matrices, respectively. Similarly, we write Vne and Ŵee as

Vne = v ne (r) n(r) dr, (4)
where v ne (r) is the nuclei-electron potential, and Ŵee = 1 2 w ee (r 12 ) n2 (r 1 , r 2 ) dr 1 dr 2 ,

where w ee (r 12 ) = 1/r 12 is the Coulomb electron-electron potential, and n(r) = ψ † (r) ψ(r) and n2 (r 1 , r 2 ) = ψ † (r 1 ) ψ † (r 2 ) ψ(r 2 ) ψ(r 1 ) are the density and pair density operators, respectively.

Introducing a set of orthonormal 4-component-spinor orbitals {ψ p (r)} which are eigenfunctions of a one-particle Dirac Hamiltonian with some potential, and assuming that this set of orbitals can be partitioned into a set of positive-energy orbitals and a set of negative-energy orbitals, {ψ p (r)} = {ψ p (r)} ε p >0 ∪ {ψ p (r)} ε p <0 , the no-pair [START_REF] Sucher | [END_REF]16] relativistic ground-state energy of a N-electron system can be defined using a minmax principle [23][24][25][26][START_REF] Saue | Theoretical Chemistry and Physics ofHeavy and Superheavy Elements[END_REF][START_REF] Almoukhalalati | [END_REF], that we will formally write as,

E 0 = minmax Ψ + Ψ + | TD + Vne + Ŵee |Ψ + . (6)
In this equation, we search over normalized wave functions of the form |Ψ + = P+ |Ψ , where P+ is the projector on the N-electron-state space generated by the set of positive-energy orbitals {ψ p (r)} ε p >0 and |Ψ is a general N-electron antisymmetric wave function, and the notation minmax Ψ + = min Ψ max P+ = max P+ min Ψ means a minimization with respect to Ψ and a maximization with respect to P+ . This maximization must be done by rotations of the positive-energy orbitals {ψ p (r)} ε p >0 with its complement set of negative-energy orbitals {ψ p (r)} ε p <0 . Here, we have assumed that the optimum of the minmax is a saddle point in the wave-function parameter space (which can be calculated with a multiconfiguration self-consistent-field (MCSCF) algorithm [START_REF] Almoukhalalati | [END_REF][START_REF] Jensen | [END_REF]32]), so that the same energy is obtained whatever the order of min Ψ and max P+ . Note that, in the nonrelativistic limit (c → ∞), the energy gap between positiveand negative-energy orbitals of order 2mc 2 goes to infinity and the maximization over P+ becomes useless, and thus the minmax principle properly reduces to the non-relativistic minimization principle. Now we attempt to formulate a relativistic DFT within this no-pair approximation. Following the spirit of the constrained-search formulation of non-relativistic DFT [33,34], we propose to define the no-pair relativistic universal density functional as

F[n] = minmax Ψ + →n Ψ + | TD + Ŵee |Ψ + = Ψ + [n]| TD + Ŵee |Ψ + [n] , (7)
where the minmax procedure is identical to that in Eq. ( 6) except for the additional constraint that Ψ + yields the density n, i.e. Ψ + |n(r)|Ψ + = n(r). In Eq. ( 7), Ψ + [n] is the optimal wave function for the density n. We will again assume that the optimum of the minmax is a saddle point in the densityconstrained wave-function parameter subspace. Of course, this functional is only defined for densities which come from a wave function of the form of Ψ + , which we will refer to as Ψ + -representable densities. Note that, consistently with neglecting the Breit electron-electron interaction, we will only consider functionals of the density and not of the density current. The no-pair relativistic ground-state energy of Eq. ( 6) can be in principle obtained from F[n] as a stationary point with respect to variations over Ψ + -representable densities

E 0 ∈ stat n F[n] + v ne (r) n(r) dr , (8)
where we have introduced the notation stat n to designates the set of stationary energies with respect to variations of n. Due to the minmax principle in Eqs. ( 6) and ( 7), we can only assume a stationary principle in Eq. ( 8), instead of the usual non-relativistic minimization principle over densities. This situation is in fact similar to the problem of formulating a pure-state time-independent variational extension of DFT for excited-state energies [35,36].

We now define a no-pair relativistic long-range universal density functional, similarly to Eq. ( 7), as

F lr, µ [n] = minmax Ψ + →n Ψ + | TD + Ŵlr, µ ee |Ψ + = Ψ µ + [n]| TD + Ŵlr, µ ee |Ψ µ + [n] , (9)
with the long-range electron-electron interaction operator Ŵlr, µ ee = (1/2) w lr, µ ee (r 12 ) n2 (r 1 , r 2 ) dr 1 dr 2 where w lr, µ ee (r 12 ) = erf(µr 12 )/r 12 is the long-range electron-electron potential and µ is the range-separation parameter. In Eq. ( 9)

Ψ µ + [n]
is the optimal wave function for the density n and range-separation parameter µ. We can thus decompose the density functional F[n] as

F[n] = F lr, µ [n] + Ēsr, µ Hxc [n], (10)
which defines the complement relativistic short-range Hartree-exchange-correlation density functional Ēsr, µ Hxc [n]. Plugging Eq. ( 10) into Eq. ( 8), we conclude that the no-pair relativistic ground-state energy of Eq. ( 6) corresponds to a stationary point of the following range-separated energy expression over Ψ + wave functions

E 0 ∈ stat Ψ + Ψ + | TD + Vne + Ŵlr, µ ee |Ψ + + Ēsr, µ Hxc [n Ψ + ] , (11)
where n Ψ + is the density of Ψ + . For practical calculations, we will assume that the no-pair relativistic ground-state energy corresponds in fact to the minmax search over

Ψ + E 0 = minmax Ψ + Ψ + | TD + Vne + Ŵlr, µ ee |Ψ + + Ēsr, µ Hxc [n Ψ + ] .( 12
)
Even though we do not see any guarantee that this is always true, it seems a reasonable working assumption for practical calculations. In fact, it corresponds to what is done in practice in no-pair Kohn-Sham DFT calculations [37][38][39][40][41][42][43][44][45][46], which corresponds to Eq. ( 12) in the special case of µ = 0, i.e.

E 0 = minmax Φ + Φ + | TD + Vne |Φ + + E Hxc [n Φ + ] , (13)
where the wave function can be restricted to a single determinant Φ + and E Hxc [n] is the relativistic Kohn-Sham Hartreeexchange-correlation density functional. Another special case of Eq. ( 12) is for µ → ∞ for which we correctly recover the wave-function theory of Eq. ( 6).

As usual, we can decompose the complement relativistic short-range Hartree-exchange-correlation density functional into separate components

Ēsr, µ Hxc [n] = E sr, µ H [n] + E sr, µ x [n] + Ēsr, µ c [n]. (14)
In this expression, E sr,µ H [n] is the short-range Hartree density functional (which has the same expression as in the nonrelativistic case)

E sr,µ H [n] = 1 2 w sr, µ ee (r 12 ) n(r 1 )n(r 2 ) dr 1 dr 2 , (15)
where w sr, µ ee (r 12 ) = w ee (r 12 )w lr, µ ee (r 12 ) is the short-range electron-electron potential, E sr, µ

x [n] is the relativistic shortrange exchange density functional

E sr,µ x [n] = Φ + [n]| Ŵsr,µ ee |Φ + [n] -E sr,µ H [n], (16)
where

Φ + [n] = Ψ µ=0 + [n]
is the relativistic Kohn-Sham single-determinant wave function and Ŵsr, µ ee = (1/2) w sr, µ ee (r 12 ) n2 (r 1 , r 2 ) dr 1 dr 2 is the short-range electron-electron interaction operator, and Ēsr, µ c [n] is the complement relativistic short-range correlation density functional. In Appendix A, we show that the relativistic short-range exchange density functional E sr,µ

x [n] satisfies a uniform coordi-nate scaling relation [Eq. (A6)] which represents an important constraint to impose in approximations.

Even though the present formulation of relativistic RS-DFT seems reasonable for practical chemical applications, it obviously calls for a closer mathematical examination of its domain of validity. In particular, it is clear that the minmax principle of the no-pair approximation in the configuration-space approach breaks down in the strong relativistic regime (i.e., for nuclear charges Z c). Of course, in the non-relativistic limit (c → ∞), relativistic RS-DFT properly reduces to nonrelativistic RS-DFT.

III. COMPUTATIONAL SETUP

We consider the helium, beryllium, neon, and argon isoelectronic series, up to the uranium nuclear charge Z = 92. The electronic density n(r) naturally increases at the nucleus with Z and can be conveniently measured with k F max , i.e. the maximal value taken at the nucleus by the local Fermi wave vector k F (r) = (3π 2 n(r)) 1/3 . The strength of the relativistic effects can be measured by comparing the local Fermi wave vector k F (r) to the speed of light c ≃ 137.036 a.u. (with = m e = 1 a.u.): very little relativistic effects are expected in regions where k F (r) ≪ c, while strong relativistic effects are expected in regions where k F (r) c.

To test the different functionals, we have first performed four-component Dirac Hartree-Fock (DHF) calculations based on the relativistic Dirac-Coulomb Hamiltonian with point-charge nucleus, using our own program implemented as a plugin of the software Quantum package 2.0 [47]. For the helium series, we use the dyall 1s2.3z basis set of Ref. 28 except for Yb 68+ and U 90+ for which the basis set was not available. For these systems, as well as for the beryllium, neon, and argon series, we construct uncontracted eventempered Gaussian-type orbital basis sets [48], following the primitive structure of the dyall-cvdz basis sets for He, Be, Ne, and Ar [49]. For each system and angular momentum, the exponents of the large-component basis functions are taken as the geometric series

ζ ν = ζ 1 q ν-1 , (17)
where ζ 1 is chosen among the largest exponents from the dyall-cvdz basis set for the given element and angular momentum [49,50], and the parameter q is optimized by minimizing the DHF total energy. The small-component basis functions are generated from the unrestricted kinetic-balance scheme [51]. The basis-set parameters are given in the Supplementary Material.

Using the previously obtained DHF orbitals, we then estimate the short-range exact exchange energy

E sr,µ x = 1 2 w sr,µ ee (r 12 ) n 2,x (r 1 , r 2 )dr 1 dr 2 , (18)
where n 2,x (r 1 , r 2 ) is the exchange pair density

n 2,x (r 1 , r 2 ) = -Tr[γ(r 1 , r 2 )γ(r 2 , r 1 )], (19)
and γ(r 1 , r 2 ) = N i=1 ψ i (r 1 )ψ † i (r 2 ) is the 4 × 4 one-electron density matrix written with the four-component spinor occupied orbitals {ψ i (r)}. This short-range DHF exchange energy is used as the reference for testing the different exchange energy functionals, which are evaluated with the DHF density n(r) = Tr[γ(r, r)] (and the DHF exchange on-top pair density for some of them, see below) using a SG-2-type quadrature grid [52] with the radial grid of Ref. 53.

IV. SHORT-RANGE EXCHANGE LOCAL-DENSITY APPROXIMATIONS

The non-relativistic short-range local-density approximation (srLDA) for the exchange functional has the expression

E sr,LDA, µ x [n] = n(r) ǫ sr,HEG, µ x (n(r)) dr, ( 20
)
where the non-relativistic short-range homogeneous electron gas (HEG) exchange energy per particle ǫ sr,HEG, µ x (n) can be found in Refs. 1, 54, and 55. The relativistic generalization of this functional, referred to as srRLDA, is

E sr,RLDA, µ x [n] = n(r) ǫ sr,RHEG, µ x (n(r)) dr, (21)
where the short-range RHEG exchange energy per particle ǫ sr,RHEG, µ x (n) is given in Ref. 11 with arbitrary accuracy as systematic Padé approximants with respect to the dimensionless variable c = c/k F = c/(3π 2 n) 1/3 (we employ here the Padé approximant of order 6) with coefficients written as functions of the dimensionless range-separation parameter µ/k F . The dependence of ǫ sr,RHEG, µ x (n) on the dimensionless parameters c and µ/k F is a consequence of the uniform coordinate scaling relation of Eq. (A6) which is valid of the RHEG.

The relative percentage errors of the srLDA and srRLDA exchange functionals with respect to the short-range DHF exchange energy, i.e. 100 × (E sr,DFA, µ

x -E sr, µ

x )/|E sr, µ x |, are plotted in Fig. 1 as a function of the dimensionless rangeseparation parameter µ/k F max for three representative members of the neon isoelectronic series (Ne, Xe 44+ and Rn 76+ ). The relativistic effects go from very small for Ne to very large for Rn 76+ .

For µ = 0, the short-range interaction reduces to the fullrange Coulomb interaction, and we observe that both the nonrelativistic and relativistic LDA exchange functionals underestimate (in absolute value) the DHF exchange energy by 5% to 10 %. As previously noted [56], the non-relativistic LDA exchange functional (evaluated with a relativistic density) fortuitously gives exchange energies with lower errors than the relativistic LDA exchange functional for systems with significant relativistic effects (Xe 44+ , and Rn 76+ ). When µ increases, the srLDA and srRLDA exchange functionals show quite different behaviors for these relativistic systems. The relative error of the srLDA exchange energy changes sign with µ and eventually goes to a negative constant for µ → ∞, corresponding to an overestimation in absolute value. By contrast, the relative error of the srRLDA exchange energy always remains positive and goes to a positive constant for µ → ∞, corresponding to an underestimation in absolute value. The more relativistic the system is, the largest this overestimation or underestimation is. While for Ne at large µ both the srLDA and srRLDA exchange functionals have almost vanishing relative errors, for Rn 76+ at large µ the srLDA exchange energy is too negative by a little more than 5% and the srRLDA exchange energy is too positive by almost 20%. Clearly, both the srLDA and srRLDA exchange functionals are not accurate for relativistic systems.

In non-relativistic theory, it is known that the srLDA exchange functional becomes exact for large µ [START_REF] Toulouse | [END_REF], which is one of the key advantages of RS-DFT. As apparent from Figure 1, for relativistic systems, this nice property does not hold anymore for both the srLDA and srRLDA exchange functionals. This observation can be understood by using the distributional asymptotic expansion of the short-range interaction for large µ

[2] w sr, µ ee (r 12 ) = π µ 2 δ(r 12 ) + O 1 µ 3 , (22)
which directly leads to the asymptotic expansion of the shortrange exact exchange energy

E sr,µ x = π 2 µ 2 n 2,x (r, r) dr + O 1 µ 4 , (23)
where n 2,x (r, r) is the on-top exchange pair density. In the non-relativistic theory, considering the case of closed-shell systems for the sake of simplicity, the on-top exchange pair density is simply given in terms of the density as [57]

n NR 2,x (r, r) = - n(r) 2 2 , (24)
and the srLDA exchange functional becomes indeed exact for large µ

E sr,LDA,µ x [n] = π 2 µ 2 n HEG,0 2,x (n(r)) dr + O 1 µ 4 , (25)
with the on-top exchange pair density of the non-relativistic HEG

n HEG,0 2,x (n) = - n 2 2 . (26)
In the relativistic theory, the on-top exchange pair density is no longer a simple function of the density

n 2,x (r, r) = -Tr[γ(r, r) 2 ], (27)
which is not equal to -n(r) 2 /2, except in the special case of two electrons in a unique Kramers pair (see Appendix B). Therefore, for relativistic systems with more than two electrons, we see that the srLDA exchange functional is not exact for large µ [Eq. (25)]. The srRLDA exchange functional is also not exact for large µ. It takes the form

E sr,RLDA,µ x [n] = π 2 µ 2 n RHEG,0 2,x (n(r)) dr + O 1 µ 4 , (28)
with the on-top exchange pair density of the RHEG

n RHEG,0 2,x (n) = - n 2 4 (1 + h(c)), ( 29
)
and the function [11] h

(c) = 9 4 c2 + c4 -c 4 arcsinh 1 c 2 √ 1 + c2 -c2 arcsinh 1 c .( 30
)
For an alternative but equivalent expression for n RHEG,0 2,x (n), see

Eq. (A1) of Ref. 39. The srRLDA exchange functional is in fact not even exact at large µ for two electrons in a unique Kramers pair. In Section V, we show how to impose the largeµ behavior on the srLDA and srRLDA exchange functionals.

V. SHORT-RANGE EXCHANGE LOCAL-DENSITY APPROXIMATIONS WITH ON-TOP EXCHANGE PAIR DENSITY

In order to impose the correct large-µ behavior of the srLDA and srRLDA exchange functionals for relativistic systems, we need to introduce a new ingredient in these functionals, namely the exact (relativistic) on-top exchange pair density n 2,x (r, r), or equivalently the on-top exchange hole

n x (r, r) = n 2,x (r, r) n(r) . (31)
A simple way to use n x (r, r) to correct the srLDA exchange functional is to find, at each position r, the effective density n eff (r) at which the on-top exchange hole of the HEG, n HEG,0 x (n) = n HEG,0 2,x (n)/n = -n/2, is equal to the ontop exchange hole of the inhomogeneous system considered, n x (r, r), i.e.

n HEG,0 x (n eff (r)) = n x (r, r), (32)
which simply gives n eff (r) = -2n x (r, r). We then define the srLDA exchange functional with the on-top exchange pair density (srLDAot) using this effective density as

E sr,LDAot, µ x [n] = n(r) ǫ sr,HEG, µ x (n eff (r)) dr. ( 33
)
This approximation could be considered either as an implicit functional of the density alone since n x (r, r) is an implicit functional of the density through the orbitals, or as an explicit functional of both the density and the on-top exchange hole n x (r, r). This approximation corresponds to changing the transferability criterion in the LDA: at a given point r, instead of taking the exchange energy per particle of the HEG having the same density than the inhomogeneous system at that point, we now take the exchange energy per particle of the HEG having the same on-top exchange hole than the inhomogeneous system at that point. Interestingly, this approximation can be thought of as a particular application of the recently formalized connector theory [START_REF] Vanzini | Auxiliary systems for observables: dynamical local connector approximation for electron addition and removal spectra[END_REF][START_REF] Vanzini | [END_REF]. Similarly, we can correct the srRLDA exchange functional by finding, at each position r, the effective density n R eff (r) at which the on-top exchange hole of the RHEG,

n RHEG,0 x (n) = n RHEG,0 2,x (n)/n = -(n/4)(1 + h(c))
, is equal to the on-top exchange hole of the inhomogeneous system considered, n x (r, r), i.e.

n RHEG,0 x (n R eff (r)) = n x (r, r). ( 34
)
This equation is less trivial to solve than Eq. ( 32) since n RHEG,0 x (n) is a complicated nonlinear function of n (through c). However, at each point r, a unique solution n R eff (r) exists since the function n → n RHEG,0 x (n) is monotonically decreasing and spans the domain ]-∞, 0] in which n x (r, r) necessarily belongs. In practice, we easily find n R eff (r) by a numerical iterative method, and we use it to define the srRLDA exchange functional with the on-top exchange pair density (srRLDAot) as

E sr,RLDAot, µ x [n] = n(r) ǫ sr,RHEG, µ x (n R eff (r)) dr. ( 35
)
Both the srLDAot and srRLDAot exchange functionals now fulfill the exact asymptotic expansion for large µ [Eq. ( 23)].

In fact, restoring the correct on-top value of the exchange hole could be beneficial for any value of µ, given the fact that the accuracy of non-relativistic Kohn-Sham exchange DFAs has been justified by the exactness of the underlying LDA on-top exchange hole (in addition to fulfilling the correct sum rule of the exchange hole) [START_REF] Burke | [END_REF]. Finally, we note that, in the nonrelativistic limit (c → ∞), we have n eff (r) = n R eff (r) = n(r) and all these short-range exchange functionals reduce to the non-relativistic srLDA exchange functional (i.e., srLDAot = srRLDAot = srRLDA = srLDA).

The relative percentage errors of the srLDAot and srRL-DAot exchange functionals for Ne, Xe 44+ , and Rn 76+ are reported in Figure 1. The most prominent feature is of course the correct recovery of the large-µ asymptotic behavior for both the srLDAot and srRLDAot exchange functionals. It turns out the srLDAot and srRLDAot exchange functionals give very similar exchange energies for all values of µ. This comes from the fact that when going from srLDA to srLDAot Eq. ( 32) tends to make the LDA exchange hole shallower and when going from srRLDA to srRLDAot Eq. ( 34) tends to make the relativistic LDA exchange hole deeper, making finally for very close descriptions. The absolute relative percentage errors of the srLDAot and srRLDAot exchange functionals are always below 10%, and below about 2% for µ/k F max ≥ 0.5.

VI. SHORT-RANGE EXCHANGE GENERALIZED-GRADIENT APPROXIMATIONS

In order to improve over the short-range LDA exchange functionals at small values of the range-separation parameter µ, we now consider short-range GGA exchange functionals. We start with the non-relativistic short-range extension of the Perdew-Burke-Ernzerhof (PBE) [61] of Refs. 62 and 63, referred to as srPBE,

E sr,PBE, µ x [n] = n(r) ǫ sr,HEG, µ x (n(r)) 1 + f µ x (n(r), ∇n(r)) dr, ( 36
)
with the function

f µ x (n, ∇n) = κ - κ 1 + b( μ)s 2 /κ , ( 37
)
where s = |∇n|/(2k F n) is the reduced density gradient and μ = µ/(2k F ) is a dimensionless range-separation parameter. In this expression, κ = 0.840 is a constant fixed by For µ = 0, this srPBE exchange functional reduces to the standard PBE exchange functional [61], and for large µ it reduces to the srLDA exchange functional.

A simple relativistic extension of this srPBE exchange functional can be obtained by replacing the srLDA part by the srRLDA one while using the same density-gradient correction f µ x (n, ∇n), to which will refer as srRLDA/PBE, (38) which reduces to the srRLDA exchange functional for large µ.

E sr,RLDA/PBE, µ x [n] = n(r) ǫ sr,RHEG, µ x (n(r)) 1 + f µ x (n(r), ∇n(r)) dr,
The srPBE and srRLDA/PBE exchange functionals have the same (incorrect) asymptotic expansions as the srLDA and srRLDA exchange functionals [Eqs. ( 25) and ( 28)], and we can thus use the same effective densities in Eqs. (32) and (34) to restore their large-µ behaviors, which defines the srPBEot and srRLDA/PBEot exchange functionals (39) where ∇n eff (r) = -2∇n x (r, r), and

E sr,PBEot, µ x [n] = n(r) ǫ sr,HEG, µ x (n eff (r)) 1 + f µ x (n eff (r), ∇n eff (r)) dr,
E sr,RLDA/PBEot, µ x [n] = n(r) ǫ sr,RHEG, µ x (n R eff (r)) 1 + f µ x (n R eff (r), ∇n R eff (r)) dr, ( 40
)
where ∇n R eff (r) = [dn RHEG,0

x (n R eff (r))/dn R eff ] -1 ∇n x (r, r). In Figure 2, we report the relative percentage errors of the srPBE, srRLDA/PBE, srPBEot, and srRLDA/PBEot exchange energies for Ne, Xe 44+ , and Rn 76+ . For Ne, where the relativistic effects are very small, all these functionals give almost the same exchange energy, as expected. For Xe 44+ and Rn 76+ , even though the srPBE and srRLDA/PBE exchange functionals are more accurate than the srLDA and srRLDA exchange functionals at µ = 0 (see Figure 1), they eventually suffer from the same large inaccuracy as srLDA and srRLDA as µ increases. This problem is solved by using the effective densities from the on-top exchange pair density, the srPBEot and srRLDA/PBEot exchange functionals giving vanishing errors at large µ. Similarly to what was observed for srLDAot and srRLDAot, the srPBEot and srRLDA/PBEot functionals give very close exchange energies for all values of µ. Interestingly, we see that using the effective densities also reduces the errors of srPBE and srRLDA/PBE at µ = 0, making srPBEot and srRLDA/PBEot quite accurate in this full-range limit. Thus, the srPBEot and srRLDA/PBEot exchange functionals are definitely an improvement over srLDAot and sr-RLDAot. We observe a maximal absolute percentage error of about 3% for Rn 76+ for µ/k F max ≈ 0.2.

In order to further reduce the errors, in particular for intermediate values of µ, we now consider a relativistic correction to the density-gradient term in the srRLDA/PBEot exchange functional. We define a short-range relativistic PBE exchange functional using the on-top exchange pair density, referred to as srRPBEot,

E sr,RPBEot, µ x [n] = n(r) ǫ sr,RHEG, µ x (n R eff (r)) × 1 + f µ x (n R eff (r), ∇n R eff (r)) φ µ (n R eff (r)) dr, (41)
where, in the spirit of the work of Engel et al. [56], we have introduced a multiplicative relativistic correction φ µ (n) to the term f µ x (n, ∇n) of the form

φ µ (n) = 1 + a 1 (µ/c) c2 + a 2 (µ/c) c4 1 + b 1 (µ/c) c2 + b 2 (µ/c) c4 . (42)
Since φ µ (n) only depends on the dimensionless parameters c and µ/c, it does not change the uniform coordinate scaling of the functional which still fulfills the scaling relation of Eq. (A6).

After some tests, we chose to impose a 1 (µ/c) = b 1 (µ/c) to avoid overcorrections in low-density regions which have very small relativistic effects. We started to determine the coefficients for µ = 0 by minimizing the mean squared relative percentage error of the exchange energy with respect to the reference DHF exchange energy for 7 systems of the neon isoelectronic series (Ne, Ar 8+ , Kr 26+ , Xe 44+ , Yb 60+ , Rn 76+ , U 82+ ), giving a 1 (0) = b 1 (0) = 1.3824, a 2 (0) = 0.3753, and b 2 = 0.4096. The resulting relativistic correction factor φ µ=0 (n) can be seen in Figure 3. It correctly tends to 1 in the low-density (k F → 0) or non-relativistic (c → ∞) limit, and remains very close to 1 for k F ≪ c. In regions with very high densities, the relativistic correction factor φ µ=0 (n) induces a slight reduction of the effective density-gradient correction term in the functional, reducing a bit the relative error on the exchange energy for the heaviest systems.

For µ 0, we have searched for coefficients in Eq. ( 42) which reduce the largest errors of the srRLDA/PBEot exchange energy observed at intermediate values of µ (see Figure 2). We chose coefficients depending on µ/c of the form

a 1 (µ/c) = b 1 (µ/c) = a 1 (0)[1 -erf(µ/c)], ( 43
)
a 2 (µ/c) = a 2 (0)[1 -erf(µ/c)], ( 44
) b 2 (µ/c) = b 2 (0)[1 -β erf(µ/c)], (45)
with β = -4.235 which has been found by minimizing the mean squared relative percentage error of the short-range exchange energy for the same 7 systems of the neon isoelectronic series and for 4 intermediate values of the range-separation parameter (µ/k F max = 0.05; 0.1; 0.2; 0.4). The resulting relativistic correction factor φ µ (n) is reported in Figure 3. It still tends to 1 in the low-density limit, but goes down to 0 when µ ≫ c in the high-density limit. The higher the value of µ the faster it decreases as a function of k F .

In Figure 4, we report the relative percentage errors of the srRPBEot exchange functional for systems of the helium, beryllium, neon, and argon series. For µ = 0, this functional achieves an error of at most about 1% for all systems, and it has the correct large-µ limit. The maximum absolute percentage errors, which are found for intermediate values of µ, tend to grow with Z but remain at most about 3% for the heavier systems. The srRPBEot exchange functional represents a significant improvement over the srPBEot and srRLDA/PBEot exchange functionals for the heavier systems.

VII. CONCLUSIONS

In this work, we have tested the srRLDA exchange functional developed in Ref. 11 on three systems of the neon isoelectronic series (Ne, Xe 44+ , and Rn 76+ ) and compared it to the usual non-relativistic srLDA exchange functional. Both functionals are quite inaccurate for relativistic systems and do not have the correct asymptotic behavior for large rangeseparation parameter µ. In order to fix this large-µ behavior, we have then defined the srLDAot and srRLDAot exchange functionals by introducing the exact on-top exchange pair density as a new variable. These functionals recover the correct asymptotic behavior for large µ but remain inaccurate for small values of µ. To improve the accuracy for small values of µ, we have then developed a relativistic short-range GGA exchange functional also using the on-top exchange pair density as an extension of the non-relativistic srPBE exchange functional. Tests on the systems of the isoelectronic series of He, Be, Ne, and Ar up to Z = 92 show that this srRPBEot exchange functional gives a maximal relative percentage error of 3% for intermediate values of µ and less than 1% relative error for µ = 0. Of course, in the non-relativistic limit (c → ∞), all the relativistic functionals introduced in this work properly reduce to their non-relativistic counterparts.

Possible continuations of this work includes further tests on atoms and molecules, extension to the Gaunt or Breit electronelectron interactions, development of the short-range relativistic correlation functionals, and use of a local range-separation parameter.

µ/k F max He (k F max =4.7) Ne 8+ (k F max =26) Ar 16+ (k F max =48) Kr 34+ (k F max =107) Xe 52+ (k F max =189) Yb 68+ (k F max =288) Rn 84+ (k F max =469) U 90+ (k F max =580) - 2
µ/k F max Be (k F max =10) Ne 6+ (k F max =26) Ar 14+ (k F max =49) Kr 32+ (k F max =110) Xe 50+ (k F max =193) Yb 66+ (k F max =301) Rn 82+ (k F max =463) U 88+ (k F max =544) - 2
µ/k F max Ne (k F max =26) Ar 8+ (k F max =49) Kr 26+ (k F max =109) Xe 44+ (k F max =192) Yb 60+ (k F max =301) Rn 76+ (k F max =462) U 82+ (k F max =543) - 2
µ/k F max Ar (k F max =49) Kr 18+ (k F max =110) Xe 36+ (k F max =193) Yb 52+ (k F max =303) Rn 68+ (k F max =467) U 74+ (k F max =550)
FIG. 4. Relative percentage error on the short-range exchange energy calculated with the srRPBEot functional for systems of helium, beryllium, neon, and argon isoelectronic series.

Appendix A: Uniform coordinate scaling relation for the relativistic no-pair short-range exchange density functional Here, we generalize the uniform coordinate scaling relation of the non-relativistic exchange density functional [66] and of the non-relativistic short-range exchange density functional [67] to the case of the relativistic no-pair short-range exchange density functional of Eq. ( 16). Since the scaling relation involves scaling the speed of light c, we will explicitly indicate in this section the dependence on c.

First, we introduce the non-interacting Dirac kinetic + rest mass energy density functional T c s [n] defined by Eq. ( 9) in the special case of a vanishing range-separation parameter, µ = 0,

T c s [n] = minmax Φ + →n Φ + | T c D |Φ + = Φ c + [n]| T c D |Φ c + [n] , (
Φ c +,γ [n]| T cγ D |Φ c +,γ [n] = γ 2 Φ c + [n]| T c D |Φ c + [n] , (
A3) and the right-hand side is minmax optimal by definition of

Φ c + [n]
. Therefore, we conclude that

Φ c +,γ [n] = Φ cγ + [n γ ]. (A4)
From the definition of the relativistic short-range exchange energy density functional E sr,µ,c

x

[n] = Φ c + [n]| Ŵsr,µ ee |Φ c + [n] - E sr,µ H [n],
we then arrive at the scaling relation

E sr,µγ,cγ x [n γ ] = γE sr,µ,c x [n], (A5)
or, equivalently,

E sr,µ,c x [n γ ] = γE sr,µ/γ,c/γ x [n]. (A6)
This scaling relation is an important constraint which is satisfied by our approximate density functionals. Besides, it shows that the low-density limit (γ → 0) corresponds to the non-relativistic limit (c → ∞), while the high-density limit (γ → ∞) corresponds to the ultra-relativistic limit (c → 0). It also shows that, for a fixed value of the range-separation parameter µ, low-density regions explore the functional in the short-range limit (µ → ∞) and high-density regions explore the functional in the full-range limit (µ = 0).      The on-top exchange pair density has the expression In the non-relativistic limit, each orbital has a definite spin state, i.e. ψ i (r) = (ψ Lα i (r), 0, 0, 0) or ψ i (r) = (0, ψ Lβ i (r), 0, 0), and we recover the well-known expression of the on-top exchange pair density in terms of the spin densities or, for closed-shell systems, n NR 2,x (r, r) = -n(r) 2 /2. However, in the relativistic case, n 2,x (r, r) can no longer be generally expressed explicitly with the density, as seen from the presence of terms mixing different spinor components in Eq. (B5). There are however two exceptions. The first exception is provided by one-electron systems for which it is easy to check that n 2,x (r, r) = -n(r) 2 , as in the non-relativistic case. The second exception is provided by systems of two electrons in an unique Kramers pair, for which n 2,x (r, r) = -n(r) 2 /2, as in the non-relativistic case. Indeed, for closed-shell systems, the one-electron density matrix can be decomposed into Kramers

ψ Lα i (
(r) =                -ψ Lβ i (r) * ψ Lα i (r) * -ψ Sβ i (r) * ψ Sα i (r) *                . ( B7
)
In this case, the density can then be expressed as n(r) = 2Tr[γ + (r, r)], and the on-top exchange pair density as The reason why systems with one electron or two electrons in a single Kramers pair constitute exceptions is that in these systems exchange only represents in fact a self-interaction correction, and we have E sr,µ x [n] = -E sr,µ H [n] for one electron and E sr,µ

x [n] = -(1/2)E sr,µ H [n] for two electrons in a single Kramers pair, as for the non-relativistic theory.

FIG. 1 .

 1 FIG.1. Relative percentage error on the short-range exchange energy calculated with the srLDA, srRLDA, srLDAot, srRLDAot functionals for three representative members of the neon isoelectronic series (Ne, Xe 44+ , and Rn 76+ ).

FIG. 2 .

 2 FIG.2. Relative percentage error on the short-range exchange energy calculated with the srPBE, srRLDA/PBE, srPBEot, srRLDA/PBEot functionals for three representative members of the neon isoelectronic series (Ne, Xe 44+ , and Rn 76+ ).

FIG. 3 .

 3 FIG.3. Relativistic correction factor φ µ (n) to the density-gradient term [Eq. (42)] as a function of k F for several values of µ.

  A1) where Φ c + [n] is the relativistic Kohn-Sham single-determinant wave function. Let us now consider the scaled wave function Φ c +,γ [n] defined by, for N electrons, Φ c +,γ [n](r 1 , ..., r N ) = γ 3N/2 Φ c + [n](γr 1 , ..., γr N ), (A2) where γ > 0 is a scaling factor. The wave function Φ c +,γ [n] yields the scaled density n γ (r) = γ 3 n(γr) and is the minmax optimal wave function of Φ + | T cγ D |Φ + since it can be checked that

n 2 ,

 2 x (r, r) = -Tr[γ(r, r) 2 ] =i (r)| 2 |ψ Lα j (r)| 2 + |ψ Lβ i (r)| 2 |ψ Lβ j (r)| 2 + 2ψ Lα i (r)ψ Lβ i (r) * ψ Lβ j (r)ψ Lα j (r) * +|ψ Sα i (r)| 2 |ψ Sα j (r)| 2 + |ψ Sβ i (r)| 2 |ψ Sβ j (r)| 2 + 2ψ Sα i (r)ψ Sβ i (r) * ψ Sβ j (r)ψ Sα j (r) * +2ψ Lα i (r)ψ Sα i (r) * ψ Sα j (r)ψ Lα j (r) * + 2ψ Lβ i (r)ψ Sβ i (r) * ψ Sβ j (r)ψ Lβ j (r) * (B4)+2ψ Lα i (r)ψ Sβ i (r) * ψ Sβ j (r)ψ Lα j (r) * + 2ψ Lβ i (r)ψ Sα i (r) * ψ Sα j (r)ψ Lβ j (r) * .

2 =

 2 -n α (r) 2n β (r) 2 , (B5)

  contributions γ(r, r) = γ + (r, r) + γ -(r, r),(B6)whereγ + (r, r) = N/2 i=1 ψ i (r) ψ † i (r) and γ -(r, r) = N/2 i=1 ψi (r) ψ † i (r), and ψi (r) is the Kramers partner of ψ i (r) ψi

n 2 ,

 2 x (r, r) = -2 Tr[γ + (r, r) 2 ] + Tr[γ + (r, r)γ -(r, r)] , (B8) where we have used Tr[γ + (r, r) 2 ] = Tr[γ -(r, r) 2 ]. For an unique Kramers pair (i.e. for N = 2), it is easy to check that Tr[γ + (r, r) 2 ] = (Tr[γ + (r, r)]) 2 and Tr[γ + (r, r)γ -(r, r)] = 0, and thus n 2,x (r, r) = -2(Tr[γ + (r, r)]) 2 = -n(r) 2 2 for N = 2. (B9)

  r)ψ Lα i (r) * ψ Lα i (r)ψ Lβ i (r) * ψ Lα i (r)ψ Sα i (r) * ψ Lα i (r)ψ Sβ i (r) * ψ Lβ i (r)ψ Lα i (r) * ψ Lβ i (r)ψ Lβ i (r) * ψ Lβ i (r)ψ Sα i (r) * ψ Lβ i (r)ψ Sβ i (r) * ψ Sα i (r)ψ Lα i (r) * ψ Sα i (r)ψ Lβ i (r) * ψ Sα i (r)ψ Sα i (r) * ψ Sα i (r)ψ Sβ i (r) * ψ + |ψ Lβ i (r)| 2 + |ψ Sα i (r)| 2 + |ψ Sβ i (r)| 2 .

				
	Sβ i (r)ψ Lα i (r) * ψ Sβ i (r)ψ Lβ i (r) * ψ Sβ i (r)ψ Sα i (r) * ψ Sβ i (r)ψ Sβ i (r) *	              	,	(B2)
	which leads to the density
		N
	n(r) = Tr[γ(r, r)] =	i=1	|ψ Lα i (r)| 2 (B3)

SUPPLEMENTARY MATERIAL

See Supplementary Material for the parameters of the eventempered basis sets constructed in this work.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.


Hua Wang

Francisco Felis

Severine Tomas

Fabien Anselmet

Muriel Amielh

An Improved Image Processing Method for Particle Characterization by Shadowgraphy

Keywords: Shadowgraphy, image processing, particle characterization

Shadowgraphy is one of the most popular imaging techniques to characterize moving particles by their size, geometry as well as velocity, due to its simplicity. However, it requires advanced image processing to handle various image defects such as non-uniform illumination, overlapped particles, etc., which are normally only solved for individual applications. This study proposes a robust image processing method for particle shadowgraphy, aiming to process imperfect particle shadow images. The proposed method first detects qualified particles from particle shadow images, and then processes detected particles individually. Therefore different defects from different particles can be handled separately and locally. An overlapped particles detection and separation algorithm is also implemented to improve the accuracy of size and geometry characterization. The proposed method is first proved by synthetic generated particle shadow images, followed by a proof test with shadow images from a transparent dot pattern target. Finally this method is successfully applied to a shadow image acquired from a water spray and proved to be able to handle various issues of shadowgraphy.

Introduction

Over past few years, there has been a growing interest and also increasing availability of commercial image analysis software for atomization characterization by means of shadowgraphy for spherical as well as irregular shape particles. Compared with point measurement techniques such as phase Doppler anemometry, it is easy to setup and the results are more straightforward to understand. Additionally, it can be applied for particles with any shape [START_REF] Bachalo | Spray diagnostics[END_REF][START_REF] Tropea | [END_REF] and it is not affected by multi-scattering effects [3,4]. However, most of imaging process methods for shadowgraphy are made for individual applications and often demand high image quality [5], which cannot be achieved in many cases due to the complexity of the spray or the limitation from optical access. In general, these image processing methods are facing various challenges like too many droplets with various shapes in the field of viewwhich results in overlapped particles [6] or out-of-focus particles, non-uniform illumination, etc. In this study, an advanced shadowgraphy image processing method is proposed to handle shadow images from complex sprays. The principle of this method will be firstly introduced, followed by verification from synthetic images, dot target images with known dot sizes, and then a validation with shadow images from a real spray.

Principle of shadow image processing

This method firstly detects particles based on a global threshold determined by the userswhich depends on individual measurement. Optionally, a 2D Mexican hat filter can be applied to convolute the raw image and consequently a map of intensity gradient can be obtained, which indicates the edge of each particle. This procedure is proved to be more effective to detect particles when the back illumination is not uniform and the particles are out of focus [7]. Then each particle is segmented from the global image (gradient image in the case of optional filtering) with a bounding box which is slightly larger than itself for further process. This process is shown in Figure 1. Then, the particle segmentation is performed from this gradient image. After the particle is segmented from raw image, a local threshold will be applied to the bounding box to extract the particle. The local threshold is calculated by:

𝑇ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑 𝑙𝑜𝑐𝑎𝑙 = 𝑖 𝑚𝑖𝑛 + (𝑖 𝑚𝑎𝑥 -𝑖 𝑚𝑖𝑛 ) * 0.61 (1)
Where imin and imax are the minimal and maximum gray values in the segmented single particle image. The factor of 0.61 is proposed by [8]. This value is the one which allows obtaining, based on point spread function, a polynomial relationship between contrast of particle image and measured particle diameter. Consequently, there is a potential to correct the measured size of those out-of-focus particles or to define a depth-of-field thickness where all particles are focused. None of these solutions is presently implemented in our algorithm. Segmented single particle image uses this threshold to be binarized and particle will be extracted for further geometry characterization. By using a local threshold for each detected particle, the issues coming from nonuniform back illumination can be minimized. In addition, using a fixed factor to determine local threshold for particle diameter evaluation can also avoid subjectivity introduced from different users. More specifically, global threshold selected by the user would only affect the identification of particlesnumber of particles identified; once the particle is identified, its size only depends on the local threshold calculated from formula (1).

After the particle extraction, its diameter and other geometry information will be calculated by standard geometry characterization functions from OpenCV. After geometry characterization, three validation methods can be applied if necessary: a) Particle area (range). b) Particle eccentricity, which is defined as the ratio between the long axis and short axis. c) Particle roundness, also called as inertia moment, which measures the spread of points around the centre of mass. Its definition is given by the formula below:

𝑖𝑛𝑒𝑟𝑡𝑖𝑎 = 𝑛𝑢 20 + 𝑛𝑢 02 (2)
𝑛𝑢 𝑗𝑖 = 𝑚𝑢 𝑗𝑖 𝑚 00 𝑖+𝑗 2 +1
(3)

𝑚𝑢 𝑗𝑖 = ∑ (𝑥 -𝑥̅ ) 𝑗 • (𝑦 -𝑦 ̅) 𝑖 𝑥,𝑦 (4)
Where x and y are coordinates of pixels within the particle, 𝑥̅ and 𝑦 ̅ are the centroid coordinates of the particle.

All the three validations above are aiming to reject noise caused by dirt in the imaging system, which are often very small in the field of view and have different shapes from expected particles.

In addition, a method to detect and separate overlapped particles is also implemented here. This method assumes that overlapped particles (seen in Figure 2) can be separated into two (or more) smaller particles by shrinking from its boundary, which represents the kernels of those overlapped particles. If two (or more) smaller particles can be identified after shrinking, these particles will be identified as overlapped particles, and a reexpand process will be applied to those kernels to calculate the size of each separated particle. The amount of reexpand is identical as the amount of shrink. The process is described in Figure 2. This overlapped particles detection and correction is based on the assumption mentioned above. In some applications, particles shape can be highly irregular and therefore they are identified as overlapped particles. In this case, the proposed method is no longer suitable.

Results and discussion

The proposed shadow image processing method has firstly been evaluated by synthetic particle image, which is shown in Figure 3, together with the contour provided by shadow image processing. The synthetic particle image is a binary image made from Microsoft Paint program, with particles having diameters from 10 to 80 pixels. The diameter results provided by shadow image processing are given in Table 1. The comparison between measured diameter from shadow processing and ground truth shows a very good agreement between each other. This proves the accuracy of particle edge detection and diameter evaluation. The overlapped particles detection and separation is also firstly tested by synthetic particles image. The left side of Figure 4 shows a group of synthetic particles with slight overlapping and the right side of the same figure shows the particle contour given by shadow processing. Overlapped particles have been well detected and recovered by this method. This method has been also tested by synthetic particles images with severe overlapping but it failed to identify the overlapping. Figure 6 shows the diameter histogram from shadow image processing, which indicates a very narrow distribution of each size group. Table 2 shows the comparison of target specifications and the results measured from shadow processing, which also shows a very good agreement. All small dots and middle dots are detected except only one middle dot which is missed because it is connected to a dirt on the target (shown in Figure 5.b). The missed big dots from shadow processing are due to the fact that they are very close to the frame border and are therefore rejected.

Different global thresholds have been also tested to study the influence of this parameter. Tests with various global thresholds show very similar result (variance < 1%). The slight difference is introduced by different numbers of particles identified with different global thresholds and therefore the statistical results are changed. Finally the proposed method has been tested with a shadow image acquired from a water spray, which is shown in Figure 7, together with the contour determined by shadow image processing. Typical challenges of shadow measurements on water spray can be all identified here, such as non-spherical particles, out-of-focus particles, overlapped particles, non-uniform back illumination, and so on. Most of the particles can be detected and segmented correctly from this image with the proposed method. Particularly, that in the right bottom side (on top of a big out-of-focus particle) has been correctly detected. In addition, the two overlapped particles next to this big out-of-focus particle have also been detected and separated correctly. There are several sharp particles which are not identified, but this is because they are close to the image boundary and therefore are rejected. Several out-of-focus particles are also detected, which can cause a bias on the size result. Correction to these out-of-focus particles is not included in this study. Fdida and Blaisot proposed a method to correct the out-of-focus particles or to define a depth-of-field thickness where all particles are focused, based on the point spread function [8]. These methods will be implemented and investigated in the next stage of this project.

Conclusions

A new shadow processing method is proposed with advanced particle detection and segmentation algorithm, together with an overlapping particles detection and separation method. This method is firstly qualified by synthetic particle images and then further evaluated with shadow images of a transparent target with dot patterns.

Later the new method is applied to shadow images from a water spray.

All test results show that this new shadow processing method can handle various issues of shadowgraphy measurement on water sprays. Particularly, this method is robust and less sensitive to parameters compared with existing shadow processing methods. In addition, this method can handle overlapped particles in the shadow image and provide more accurate results.

Figure 1 .

 1 Figure 1. Particle detection and segmentation: (a) general routine, (b) Optional Mexican hat filter to be used instead of global threshold

Figure 2 .

 2 Figure 2. Method to detect and correct overlapped particles

Figure 3 .

 3 Figure 3. Synthetic particles image and contours (red) determined by shadow processing

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Left, synthetic overlapped particles; Right, contour of each particle detected by shadow image processing

Figure 7 .

 7 Figure 7. Shadow images acquired from a water spray, together with red contour of each detected particle given by shadow image processing

Table 1 .

 1 Comparison of the ground truth and diameter measured by shadow processing on the synthetic particles image

	Ground truth (pix)	Diameter measured by shadow processing (pix)
	10	9.933
	20	19.93
	30	29.821
	40	39.960
	50	49.861
	60	59.869
	70	69.901
	80	79.869

Table 2 .

 2 Summary of measurement results from shadow imaging processing on the transparent dot target shadow image

	Dot Size (μm)	Number of Dots	Dots detected	Average size measured	Standard deviation of size
				by shadow processing	measured by shadow
				(μm)	processing (μm)
	62.5	1681	1681	63.98	0.92
	125	1240	1239	127.75	1.08
	250	400	371	255.73	0.90

Acknowledgements

This work is supported by a collaboration project among IRSTEA, IRPHE and Dantec Dynamics A/S, aiming to improve characterization of moving irregular particles by shadowgraphy measurements.


est destinée au dépôt

  , Jérôme Gaillardet, Isabelle Braud, Olivier Bour, Nathalie Dörfliger, C. Galy, Fatima Laggoun-Défarge, Guillaume Nord, D. Six, Tiphaine Tallec To cite this version: F. Hankard, Jérôme Gaillardet, Isabelle Braud, Olivier Bour, Nathalie Dörfliger, et al.. OZCAR: The French Critical Zone initiative: connecting the Earth's skin observatories. EGU General Assembly 2018, Apr 2018, Vienna, Austria. pp.1, 2018. hal-02608274