Pierre-François Loos

Xavier Blase

Dynamical Correction to the Bethe-Salpeter Equation Beyond the Plasmon-Pole Approximation

The Bethe-Salpeter equation (BSE) formalism is a computationally affordable method for the calculation of accurate optical excitation energies in molecular systems. Similar to the ubiquitous adiabatic approximation of time-dependent density-functional theory, the static approximation, which substitutes a dynamical (i.e., frequency-dependent) kernel by its static limit, is usually enforced in most implementations of the BSE formalism. Here, going beyond the static approximation, we compute the dynamical correction of the electron-hole screening for molecular excitation energies thanks to a renormalized first-order perturbative correction to the static BSE excitation energies. The present dynamical correction goes beyond the plasmon-pole approximation as the dynamical screening of the Coulomb interaction is computed exactly within the random-phase approximation. Our calculations are benchmarked against high-level (coupled-cluster) calculations, allowing to assess the clear improvement brought by the dynamical correction for both singlet and triplet optical transitions.

I. INTRODUCTION

The Bethe-Salpeter equation (BSE) formalism 1,2 is to the GW approximation 3,4 of many-body perturbation theory (MBPT) [START_REF] Onida | [END_REF][START_REF] Martin | Interacting Electrons: Theory and Computational Approaches[END_REF] what time-dependent density-functional theory (TD-DFT) [START_REF] Runge | [END_REF][START_REF] Casida | Time-dependent density functional response theory for molecules[END_REF] is to Kohn-Sham density-functional theory (KS-DFT), [START_REF] Hohenberg | [END_REF]10 an affordable way of computing the neutral (or optical) excitations of a given electronic system. In recent years, it has been shown to be a valuable tool for computational chemists with a large number of systematic benchmark studies on large families of molecular systems appearing in the literature [11][12][13][14][15][16][17][18][19][20][21] (see Ref. 22 for a recent review).

Qualitatively, taking the optical gap (i.e., the lowest optical excitation energy) as an example, BSE builds on top of a GW calculation by adding up excitonic effects (i.e., the electronhole binding energy E B ) to the GW HOMO-LUMO gap

E GW g = ε GW LUMO -ε GW HOMO , (1)
which is itself a corrected version of the Kohn-Sham (KS) gap

E KS g = ε KS LUMO -ε KS HOMO E GW g ≈ E fund g , (2)
in order to approximate the optical gap

E opt g = E N 1 -E N 0 = E fund g + E B , (3)
where

E fund g = I N -A N (4)
is the fundamental gap, I N = E N-1 0 -E N 0 and A N = E N 0 -E N+1 0 being the ionization potential and the electron affinity of the N-electron system, respectively. Here, E N S is the total energy of the S th excited state of the N-electron system, and E N 0 corresponds to its ground-state energy. Because the excitonic effect corresponds physically to the stabilization implied by a) Electronic mail: loos@irsamc.ups-tlse.fr b) Electronic mail: xavier.blase@neel.cnrs.fr the attraction of the excited electron and its hole left behind, we have E opt g < E fund g . Due to the smaller amount of screening in molecules as compared to solids, a faithful description of excitonic effects is paramount in molecular systems.

Most of BSE implementations rely on the so-called static approximation, which approximates the dynamical (i.e., frequency-dependent) BSE kernel by its static limit. In complete analogy with the ubiquitous adiabatic approximation in TD-DFT where the exchange-correlation (xc) kernel is made static, one key consequence of the static approximation within BSE is that double (and higher) excitations are completely absent from the BSE spectrum. Indeed, a frequency-dependent kernel has the ability to create additional poles in the response function, which describe states with a multiple-excitation character, and, in particular, double excitations. Although these double excitations are usually experimentally dark (which means that they usually cannot be observed in photo-absorption spectroscopy), these states play, indirectly, a key role in many photochemistry mechanisms, 23 as they strongly mix with the bright singly-excited states leading to the formation of satellite peaks. 24,25 They are particularly important in the faithful description of the ground state of open-shell molecules, [26][27][28][29] and they are, moreover, a real challenge for high-level computational methods. [29][30][31][32] Double excitations play also a significant role in the correct location of the excited states of polyenes that are closely related to rhodopsin, a biological pigment found in the rods of the retina and involved in the visual transduction. [START_REF] Olivucci | Computational Photochemistry[END_REF][START_REF] Robb | A Computational Strategy for Organic Photochemistry[END_REF][START_REF] Manathunga | [END_REF] In butadiene, for example, while the bright 1 1 B u state has a clear HOMO → LUMO single-excitation character, the dark 2 1 A g state includes a substantial fraction of doubly-excited character from the HOMO 2 → LUMO 2 double excitation (roughly 30%), yet with dominant contributions from the HOMO -1 → LUMO and HOMO → LUMO + 1 single excitations. 31,[36][37][38][39][40][41][42] Going beyond the static approximation is difficult and very few groups have been addressing the problem. 2,28, [43][44][45][46][47][48][49][50][51][52][53][54] Nonetheless, it is worth mentioning the seminal work of Strinati (who originally derived the dynamical correction to the BSE) on core excitons in semiconductors, 2, 55,56 in which the dynamical screening effects were taken into account through the dielectric matrix, and where he observed an increase of the binding energy over its value for static screening and a narrowing of the Auger width below its value for a core hole. Following Strinati's footsteps, Rohlfing and coworkers have developed an efficient way of taking into account, thanks to firstorder perturbation theory, the dynamical effects via a plasmonpole approximation combined with the Tamm-Dancoff approximation (TDA). 43,46,47,57 With such a scheme, they have been able to compute the excited states of biological chromophores, showing that taking into account the electron-hole dynamical screening is important for an accurate description of the lowest n → π * excitations. 46,47,57 Indeed, studying PYP, retinal and GFP chromophore models, Ma et al. found that "the influence of dynamical screening on the excitation energies is about 0.1 eV for the lowest π → π * transitions, but for the lowest n → π * transitions the influence is larger, up to 0.25 eV." 47 A similar conclusion was reached in Ref. 46. Zhang et al. have studied the frequency-dependent second-order Bethe-Salpeter kernel and they have observed an appreciable improvement over configuration interaction with singles (CIS), time-dependent Hartree-Fock (TDHF), and adiabatic TD-DFT results. 51 Rebolini and Toulouse have performed a similar investigation in a range-separated context, and they have reported a modest improvement over its static counterpart. 52,[START_REF] Rebolini | Range-Separated Density-Functional Theory for Molecular Excitation Energies[END_REF] In these two latter studies, they also followed a (non-self-consistent) perturbative approach within the TDA with a renormalization of the first-order perturbative correction.

It is important to note that, although all the studies mentioned above are clearly going beyond the static approximation of BSE, they are not able to recover additional excitations as the perturbative treatment accounts for dynamical effects only on excitations already present in the static limit. However, it does permit to recover, for transitions with a dominant singleexcitation character, additional relaxation effects coming from higher excitations. These higher excitations would be explicitly present in the BSE Hamiltonian by "unfolding" the dynamical BSE kernel, and one would recover a linear eigenvalue problem with, nonetheless, a much larger dimension. [START_REF] Loos | [END_REF] Based on a simple two-level model which permits to analytically solve the dynamical equations, Romaniello and coworkers 48,49 evidenced that one can genuinely access additional excitations by solving the non-linear, frequencydependent eigenvalue problem. For this particular system, it was shown that a BSE kernel based on the random-phase approximation (RPA) produces indeed double excitations but also unphysical excitations. 48 The appearance of these spurious excitations was attributed to the self-screening problem. 27 This was fixed in a follow-up paper by Sangalli et al. 49 thanks to the design of a number-conserving approach based on the folding of the second-RPA Hamiltonian, [START_REF] Wambach | [END_REF] which includes explicitly both single and double excitations. By computing the polarizability of two unsaturated hydrocarbon chains, C 8 H 2 and C 4 H 6 , they showed that their approach produces the correct number of physical excitations.

Finally, let us mention efforts to borrow ingredients from BSE in order to go beyond the adiabatic approximation of TD-DFT. For example, Huix-Rotllant and Casida 26,28 proposed a nonadiabatic correction to the xc kernel using the formalism of superoperators, which includes as a special case the dressed TD-DFT method of Maitra and coworkers, 25,36,37,[START_REF] Maitra | Memory: History , initial-state dependence , and doubleexcitations[END_REF] where a frequency-dependent kernel is build a priori and manually for a particular excitation. Following a similar strategy, Romaniello et al. 48 took advantages of the dynamically-screened Coulomb potential from BSE to obtain a dynamic TD-DFT kernel. In this regard, MBPT provides key insights about what is missing in adiabatic TD-DFT, as discussed in details by Casida and Huix-Rotllant in Ref. 62.

In the present study, we extend the work of Rohlfing and coworkers 43,46,47,57 by proposing a renormalized first-order perturbative correction to the static BSE excitation energies. Importantly, our correction goes beyond the plasmon-pole approximation as the dynamical screening of the Coulomb interaction is computed exactly. In order to assess the accuracy of the present scheme, we report singlet and triplet excitation energies of various natures for small-and medium-size molecules. Our calculations are benchmarked against high-level coupledcluster (CC) calculations, allowing to clearly evidence the systematic improvement brought by the dynamical correction. In particular, we found that, although n → π * and π → π * transitions are systematically red-shifted by 0.3-0.6 eV, dynamical effects have a much smaller magnitude for charge transfer (CT) and Rydberg states. Unless otherwise stated, atomic units are used.

II. THEORY

In this Section, following Strinati's seminal work, 2 we first discuss in some details the theoretical foundations leading to the dynamical BSE. We present, in a second step, the perturbative implementation of the dynamical correction as compared to the standard static approximation.

A. General dynamical BSE

The two-body correlation function L(1, 2; 1 , 2 ) -a central quantity in the BSE formalism -relates the variation of the one-body Green's function G(1, 1 ) with respect to an external non-local perturbation U(2 , 2), i.e., iL(1,

2; 1 , 2 ) = ∂G(1, 1 ) ∂U(2 , 2) , (5)
where, e.g., 1 ≡ (x 1 t 1 ) is a space-spin plus time composite variable. The relation between G and the one-body charge density ρ(1) = -iG(1, 1 + ) provides a direct connection with the density-density susceptibility χ(1, 2) = L(1, 2; 1 + , 2 + ) at the core of TD-DFT. (The notation 1 + means that the time t 1 is taken at t + 1 = t 1 + 0 + , where 0 + is a positive infinitesimal.) The two-body correlation function L satisfies the selfconsistent BSE 2 L(1, 2; 1 , 2 ) = L 0 (1, 2; 1 , 2 ) + d3456 L 0 (1, 4; 1 , 3)Ξ(3, 5; 4, 6)L(6, 2; 5, 2 ), (6) where iL 0 (1, 4; 1 , 3) = G(1, 3)G(4, 1 ), (7a) iL(1, 2; 1 , 2 ) = -G 2 (1, 2; 1 , 2 ) + G(1, 1 )G(2, 2 ), (7b) can be expressed as a function of the one-and two-body Green's functions

G(1, 2) = -i N|T [ ψ(1) ψ † (2)]|N , (8a)
G 2 (1, 2; 1 , 2 ) = -N|T [ ψ(1) ψ(2) ψ † (2 ) ψ † (1 )]|N , (8b)
and

Ξ(3, 5; 4, 6) = i δ[v H (3)δ(3, 4) + Σ xc (3, 4)] δG(6, 5) (9)
is the BSE kernel that takes into account the self-consistent variation of the Hartree potential

v H (1) = -i d2 v(1, 2)G(2, 2 + ), (10)
[where δ is Dirac's delta function and v is the bare Coulomb operator] and the xc self-energy Σ xc with respect to the variation of G. In Eqs. (8a) and (8b), the field operators ψ(xt) and ψ † (x t ) remove and add (respectively) an electron to the N-electron ground state |N in space-spin-time positions (xt) and (x t ), while T is the time-ordering operator.

The resolution of the dynamical BSE starts with the expansion of L 0 and L [see Eqs. (7a) and (7b)] over the complete orthonormalized set of N-electron excited states |N, S (with |N, 0 ≡ |N ). 2 In the optical limit of instantaneous electronhole creation and destruction, imposing t 2 = t + 2 and t 1 = t + 1 , and using the relation between the field operators in their timedependent (Heisenberg) and time-independent (Schrödinger) representations, e.g.,

ψ(1) = e i Ĥt 1 ψ(x 1 )e -i Ĥt 1 , (11)
( Ĥ being the exact many-body Hamiltonian), one gets

iL(1, 2; 1 , 2 ) = θ(+τ 12 ) s>0 χ S (x 1 , x 1 ) χS (x 2 , x 2 )e -iΩ S τ 12 -θ(-τ 12 ) s>0 χ S (x 2 , x 2 ) χS (x 1 , x 1 )e +iΩ S τ 12 , (12)
where τ 12 = t 1t 2 , θ is the Heaviside step function, and

χ S (x 1 , x 1 ) = N|T [ ψ(x 1 ) ψ † (x 1 )]|N, S , (13a) χS
(x 1 , x 1 ) = N, S |T [ ψ(x 1 ) ψ † (x 1 )]|N . (13b)
The Ω s 's are the neutral excitation energies of interest (with

Ω s = E N s -E N 0 )
. Picking up the e +iΩ S t 2 component of both L(1, 2; 1 , 2 ) and L(6, 2; 5, 2 ), simplifying further by χS (x 2 , x 2 ) on both sides of the BSE [see Eq. ( 6)], we seek the e -iΩ S t 1 Fourier component associated with the right-hand side of a modified dynamical BSE, which reads

N|T [ ψ(x 1 ) ψ † (x 1 )]|N, S e -iΩ S t 1 θ(τ 12 ) = d3456 L 0 (1, 4; 1 , 3)Ξ(3, 5; 4, 6) × N|T [ ψ(6) ψ † (5)]|N, S θ[min(t 5 , t 6 ) -t 2 ]. ( 14
)
For the neutral excitation energies falling in the fundamental gap of the system (i.e., Ω S < E fund g due to excitonic effects), L 0 (1, 2; 1 , 2 ) cannot contribute to the e -iΩ S t 1 response term since its lowest excitation energy is precisely the fundamental gap [see Eq. ( 4)]. Consequently, special care has to be taken for high-lying excited states (like core or Rydberg excitations) where additional terms have to be taken into account (see Refs. 55 and 56).

Dropping the space/spin variables, the Fourier components with respect to t 1 of L 0 (1, 4; 1 , 3) reads

[iL 0 ](ω 1 ) = dω 2π G ω - ω 1 2 G ω + ω 1 2 e iωτ 34 e iω 1 t 34 , (15)
with τ 34 = t 3t 4 and t 34 = (t 3 + t 4 )/2. We now adopt the Lehman representation of the one-body Green's function in the quasiparticle approximation, i.e.,

G(x 1 , x 2 ; ω) = p φ p (x 1 )φ * p (x 2 ) ω -ε p + iη × sgn(ε p -µ) , (16)
where η is a positive infinitesimal and µ is the chemical potential. The ε p 's in Eq. ( 16) are quasiparticle energies (i.e., proper addition/removal energies) and the φ p (x)'s are their associated one-body (spin)orbitals. In the following, i and j are occupied orbitals, a and b are unoccupied orbitals, while p, q, r, and s indicate arbitrary orbitals. Projecting the Fourier component

L 0 (x 1 , 4; x 1 , 3; ω 1 = Ω S ) onto φ * a (x 1 )φ i (x 1 ) yields dx 1 dx 1 φ * a (x 1 )φ i (x 1 )L 0 (x 1 , 4; x 1 , 3; Ω S ) = φ * a (x 3 )φ i (x 4 )e iΩ S t 34 Ω S -(ε a -ε i ) + iη θ(τ 34 )e i ε i + Ω S 2 τ 34 + θ(-τ 34 )e i ε a - Ω S 2 τ 34 . (17)
More details are provided in Appendix A. As a final step, we express the terms N|T [ ψ(x 1 ) ψ † (x 1 )]|N, S and N|T [ ψ(6) ψ † ( 5)]|N, S from Eq. ( 14) in the standard electronhole product (or single-excitation) space. This is done by expanding the field operators over a complete orbital basis of creation/destruction operators. For example, we have (see derivation in Appendix B)

N|T [ ψ(6) ψ † (5)]|N, S = -e -iΩ S t 65 pq φ p (x 6 )φ * q (x 5 ) N|â † q âp |N, S × θ(τ 65 )e -i ε p - Ω S 2 τ 65 + θ(-τ 65 )e -i ε q + Ω S 2 τ 65 , (18)
with t 65 = (t 5 + t 6 )/2 and τ 65 = t 6t 5 . The N|â † q âp |N, S are the unknown particle-hole amplitudes.

B. Dynamical BSE within the GW approximation

Adopting now the GW approximation 3 for the xc self-energy, i.e.,

Σ GW xc (1, 2) = iG(1, 2)W(1 + , 2), (19)
leads to the following simplified BSE kernel

Ξ(3, 5; 4, 6) = v(3, 6)δ(3, 4)δ(5, 6) -W(3 + , 4)δ(3, 6)δ(4, 5), ( 20
)
where W is the dynamically-screened Coulomb operator. The GW quasiparticle energies ε GW p are usually good approximations to the removal/addition energies ε p introduced in Eq. ( 16).

Substituting Eqs. ( 17), (18), and (20) into Eq. ( 14), and projecting onto φ * a (x 1 )φ i (x 1 ), one gets after a few tedious manipulations the dynamical BSE:

(ε GW a -ε GW i -Ω S )X ia,S + jb κ(ia| jb) -W i j,ab (Ω S ) X jb,S + jb κ(ia|b j) -W ib,a j (Ω S ) Y jb,S = 0, (21) with X jb,S = N|â † j âb |N, S and Y jb,S = N|â † b â j |N, S
, and where κ = 2 or 0 for singlet and triplet excited states (respectively). This equation is identical to the one presented by Rohlfing and coworkers. 43,46,47 Neglecting the anti-resonant terms, Y jb,S , in the dynamical BSE, which are (usually) much smaller than their resonant counterparts, X jb,S , leads to the well-known TDA. In Eq. ( 21), (pq|rs) = drdr φ p (r)φ q (r)v(rr )φ r (r )φ s (r ), (22) are the bare two-electron integrals in the (real-valued) spatial orbital basis {φ p (r)}, and

W pq,rs (Ω S ) = i 2π dω e -iω0 + W pq,rs (ω) ×       1 Ω S ps -ω + iη + 1 Ω S qr + ω + iη       , (23)
is an effective dynamically-screened Coulomb potential, 48 where

Ω S pq = Ω S -(ε GW q -ε GW p ) and W pq,rs (ω) = drdr φ p (r)φ q (r)W(r, r ; ω)φ r (r )φ s (r ). (24)

C. Dynamical screening

In the present study, we consider the exact spectral representation of W at the RPA level consistently with the underlying GW calculation:

W i j,ab (ω) = (i j|ab) + 2 m [i j|m][ab|m] × 1 ω -Ω RPA m + iη - 1 ω + Ω RPA m -iη , (25)
where m labels single excitations, and

[pq|m] = ia (pq|ia)(X RPA m + Y RPA m ) ia (26)
are the spectral weights. In Eqs. ( 25) and (26), Ω RPA

A RPA B RPA -B RPA -A RPA • X RPA m Y RPA m = Ω RPA m X RPA m Y RPA m , (27)
with

A RPA ia, jb = δ i j δ ab (ε a -ε i ) + 2(ia| jb), (28a)
B RPA ia, jb = 2(ia|b j), (28b)
where the ε p 's are taken as the HF orbital energies in the case of G 0 W 0 63,64 or as the GW quasiparticle energies in the case of self-consistent schemes such as evGW. 20,[64][65][66][67][68][69] The RPA matrices A RPA and B RPA in Eq. ( 27) are of size OV × OV, where O and V are the number of occupied and virtual orbitals (i.e., N orb = O + V is the total number of spatial orbitals), respectively, and X RPA m , and Y RPA m are (eigen)vectors of length OV.

The analysis of the poles of the integrand in Eq. ( 23) yields

W i j,ab (Ω S ) = (i j|ab) + 2 m [i j|m][ab|m] ×        1 Ω S ib -Ω RPA m + iη + 1 Ω S ja -Ω RPA m + iη        . ( 29
)
One can verify that, in the static limit where Ω RPA m → ∞, the matrix elements W i j,ab correctly reduce to their static expression

W stat i j,ab ≡ W i j,ab (ω = 0) = (i j|ab) -4 m [i j|m][ab|m] Ω RPA m , (30)
evidencing that the standard static BSE problem is recovered from the present dynamical formalism in this limit. Due to excitonic effects, the lowest BSE excitation energy, Ω 1 , stands lower than the lowest RPA excitation energy, Ω RPA 1 , so that, Ω S ib -Ω RPA m < 0 and W i j,ab (Ω S ) has no resonances. This property holds for low-lying excitations but special care must be taken for higher ones. Furthermore, Ω S ib and Ω S ja are necessarily negative quantities for in-gap low-lying BSE excitations. Thus, we have

|Ω S ib -Ω RPA m | > Ω RPA m .
As a consequence, we observe a reduction of the electron-hole screening, i.e., an enhancement of electron-hole binding energy, as compared to the standard static BSE, and consequently smaller (red-shifted) excitation energies. This will be numerically illustrated in Sec. IV.

D. Dynamical Tamm-Dancoff approximation

The analysis of the (off-diagonal) screened Coulomb potential matrix elements multiplying the Y jb,S coefficients in Eq. ( 21), i.e.,

W ib,a j (Ω S ) = (ib|a j) + 2 m [ib|m][a j|m] ×        1 Ω S i j -Ω RPA m + iη + 1 Ω S ba -Ω RPA m + iη        , (31)
reveals strong divergences even for low-lying excitations when, for example,

Ω S ba -Ω RPA m = Ω S -Ω RPA m -(ε GW a -ε GW b ) ≈ 0.
Such divergences may explain that, in previous studies, dynamical effects were only accounted for at the TDA level. 2,43,[46][47][48][49]51,52 To avoid confusions here, enforcing the TDA for the dynamical correction (which corresponds to neglecting the dynamical correction originating from the anti-resonant part of the BSE Hamiltonian) will be labeled as dTDA in the following. Going beyond the dTDA is outside the scope of the present study but shall be addressed eventually.

E. Perturbative dynamical correction

From a more practical point of view, Eq. ( 21) can be recast as an non-linear eigenvalue problem and, to compute the BSE excitation energies of a closed-shell system, one must solve the following dynamical (i.e., frequency-dependent) response problem 2

A(Ω S ) B(Ω S ) -B(-Ω S ) -A(-Ω S ) • X S Y S = Ω S X S Y S , (32)
where the dynamical matrices A and B have the same OV ×OV size than their RPA counterparts, and we assume real quantities from hereon. Same comment applies to the eigenvectors X S , and Y S of length OV. Note that, due to its non-linear nature, Eq. ( 32) may provide more than one solution for each value of S . [START_REF] Martin | Interacting Electrons: Theory and Computational Approaches[END_REF]48,49 Accordingly to Eq. ( 21), the BSE matrix elements in Eq. ( 32) read

A ia, jb (Ω S ) = δ i j δ ab (ε GW a -ε GW i ) + κ(ia| jb) -W i j,ab (Ω S ), (33a)
B ia, jb (Ω S ) = κ(ia|b j) -W ib,a j (Ω S ). ( 33b
)
Now, let us decompose, using basic Rayleigh-Schrödinger perturbation theory, the non-linear eigenproblem (32) as a zeroth-order static (hence linear) reference and a first-order dynamic (hence non-linear) perturbation, such that

A(Ω S ) B(Ω S ) -B(-Ω S ) -A(-Ω S ) = A (0) B (0) -B (0) -A (0) + A (1) (Ω S ) B (1) (Ω S ) -B (1) (-Ω S ) -A (1) (-Ω S ) , (34)
with

A (0) ia, jb = δ i j δ ab (ε GW a -ε GW i ) + κ(ia| jb) -W stat i j,ab , (35a)
B (0) ia, jb = κ(ia|b j) -W stat ib,a j . (35b)
and

A (1) ia, jb (Ω S ) = -W i j,ab (Ω S ) + W stat i j,ab , (36a)
B (1) ia, jb (Ω S ) = -W ib,a j (Ω S ) + W stat ib,a j . (36b)
According to perturbation theory, the S th BSE excitation energy and its corresponding eigenvector can then be expanded as

Ω S = Ω (0) S + Ω (1) S + . . . , (37a)
X S Y S = X (0) S Y (0) S + X (1) S Y (1) S + . . . . (37b)
Solving the zeroth-order static problem

A (0) B (0) -B (0) -A (0) • X (0) S Y (0) S = Ω (0) S X (0) S Y (0) S , (38)
yields the zeroth-order (static) Ω (0) S excitation energies and their corresponding eigenvectors X (0) S and Y (0) S . Thanks to firstorder perturbation theory, the first-order correction to the S th excitation energy is

Ω (1) S = X (0) S Y (0) S • A (1) (Ω (0) S ) B (1) (Ω (0) S ) -B (1) (-Ω (0) S ) -A (1) (-Ω (0) S ) • X (0) S Y (0) S . (39)
From a practical point of view, if one enforces the dTDA, we obtain the very simple expression

Ω (1) S = (X (0) S ) • A (1) (Ω (0) S ) • X (0) S . (40)
This correction can be renormalized by computing, at basically no extra cost, the renormalization factor which reads, in the dTDA,

Z S =         1 -(X (0) S ) • ∂A (1) (Ω S ) ∂Ω S Ω S =Ω (0) S • X (0) S         -1 . ( 41
)
This finally yields

Ω dyn S = Ω stat S + ∆Ω dyn S = Ω (0) S + Z S Ω (1) S . (42)
with Ω stat S ≡ Ω (0) S and ∆Ω dyn S = Z S Ω (1) S . This is our final expression. As mentioned in Sec. I, the present perturbative scheme does not allow to access double excitations as only excitations calculated within the static approach can be dynamically corrected. We hope to report a genuine dynamical treatment of the BSE in a forthcoming work.

In terms of computational cost, if one decides to compute the dynamical correction of the M lowest excitation energies, one must perform, first, a conventional (static) BSE calculation and extract the M lowest eigenvalues and their corresponding eigenvectors [see Eq. (38)]. These are then used to compute the first-order correction from Eq. (40), which also require to construct and evaluate the dynamical part of the BSE Hamiltonian for each excitation one wants to dynamically correct. The static BSE Hamiltonian is computed once during the static BSE calculation and does not dependent on the targeted excitation.

Searching iteratively for the lowest eigenstates, via Davidson's algorithm for instance, can be performed in O(N 4 orb ) computational cost. Constructing the static and dynamic BSE Hamiltonians is much more expensive as it requires the complete diagonalization of the (OV × OV) RPA linear response matrix [see Eq. ( 27)], which corresponds to a

O(O 3 V 3 ) = O(N 6
orb ) computational cost. Although it might be reduced to O(N 4 orb ) operations with standard resolutionof-the-identity techniques, 70,71 this step is the computational bottleneck in the current implementation.

III. COMPUTATIONAL DETAILS

All systems under investigation have a closed-shell singlet ground state. We then adopt a restricted formalism throughout this work. The GW calculations performed to obtain the screened Coulomb operator and the quasiparticle energies are done using a (restricted) HF starting point. Perturbative GW (or G 0 W 0 ) 63,64,72 quasiparticle energies are employed as starting points to compute the BSE neutral excitations. These quasiparticle energies are obtained by linearizing the frequencydependent quasiparticle equation, and the entire set of orbitals is corrected. Further details about our implementation of G 0 W 0 can be found in Refs. 73 and 74. Note that, for the present (small) molecular systems, G 0 W 0 @HF and evGW@HF yield similar quasiparticle energies and fundamental gap. Moreover, G 0 W 0 allows to avoid rather laborious iterations as well as the significant additional computational effort of evGW. In the present study, the zeroth-order Hamiltonian [see Eq. (34)] is always the "full" BSE static Hamiltonian, i.e., without TDA. The dynamical correction, however, is computed in the dTDA throughout. As one-electron basis sets, we employ the Dunning families cc-pVXZ and aug-cc-pVXZ (X = D, T, and Q) defined with cartesian Gaussian functions. Finally, the infinitesimal η is set to 100 meV for all calculations. It is important to mention that the small molecular systems considered here are particularly challenging for the BSE formalism, 15,73 which is known to work best for larger systems where the amount of screening is more important. 17,18 For comparison purposes, we employ the theoretical best estimates (TBEs) and geometries of Refs. 30-32 from which CIS(D), 75,76 ADC(2), 77,78 CC2, 79 CCSD, 80 and CC3 81 excitation energies are also extracted. Various statistical quantities are reported in the following: the mean signed error (MSE), mean absolute error (MAE), root-mean-square error (RMSE), and the maximum positive [Max(+)] and maximum negative [Max(-)] errors. All the static and dynamic BSE calculations have been performed with the software QuAcK, [START_REF] Loos | QuAcK: a software for emerging quantum electronic structure methods[END_REF] freely available on github, where the present perturbative correction has been implemented.

IV. RESULTS AND DISCUSSION

First, we investigate the basis set dependence of the dynamical correction. The singlet and triplet excitation energies of the nitrogen molecule N 2 computed at the BSE@G 0 W 0 @HF level for the cc-pVXZ and aug-cc-pVXZ families of basis sets are reported in Table I, where we also report the GW gap, E GW g , to show that corrected transitions are usually well below this gap. The N 2 molecule is a very convenient example for this kind of study as it contains n → π * and π → π * valence excitations as well as Rydberg transitions. As we shall further illustrate below, the magnitude of the dynamical correction is characteristic of the type of transitions. One key result of the present investigation is that the dynamical correction is quite basis set insensitive with a maximum variation of 0.03 eV between aug-cc-pVDZ and aug-cc-pVQZ. It is only for the smallest basis set (cc-pVDZ) that one can observe significant differences. We note further that due to its unbound LUMO, the GW gap of N 2 , and to a lesser extent its BSE excitation energies, are very sensitive to the presence of diffuse orbitals. However, the dynamical correction is again very stable, being insensitive to the presence of diffuse orbitals (at least for the lowest optical excitations). We can then safely conclude that the dynamical correction converges rapidly with respect to the size of the one-electron basis set, a triple-ζ or an augmented double-ζ basis being enough to obtain near complete basis set limit values. This is quite a nice feature as it means that one does not need to compute the dynamical correction in a very large basis to get a meaningful estimate of its magnitude.

Tables II and III report, respectively, singlet and triplet excitation energies for various molecules computed at the BSE@G 0 W 0 @HF level and with the aug-cc-pVTZ basis set. For comparative purposes, excitation energies obtained with the same basis set and several second-order wave function methods [CIS(D), ADC(2), CC2, and CCSD] are also reported. The highly-accurate TBEs of Refs. 30-32 (computed in the same basis) will serve us as reference, and statistical quantities [MAE, MSE, RMSE, Max(+), and Max(-)] are defined with respect to these references. For each excitation, we report the static and dynamic excitation energies, Ω stat S and Ω dyn S , as well as the value of the renormalization factor Z S defined in Eq. (41). As one can see in Tables II andIII, the value of Z S is always quite close to unity which shows that the perturbative expansion behaves nicely, and that a first-order correction is probably quite a good estimate of the non-perturbative result. Moreover, we have observed that an iterative, self-consistent resolution [where the dynamically-corrected excitation energies are re-injected in Eq. ( 39)] yields basically the same results as its (cheaper) renormalized version. Note that, unlike in GW where the renormalization factor lies in between 0 and 1, the dynamical BSE renormalization factor Z S defined in Eq. ( 41) can be smaller or greater than unity. A clear general trend is the consistent red shift of the static BSE excitation energies induced by the dynamical correction, as anticipated in Sec. II C.

The results gathered in Tables II and III are depicted in Fig. 1, where we report the error (with respect to the TBEs) for the singlet and triplet excitation energies computed within the static and dynamic BSE formalism. From this figure, it is quite clear that the dynamically-corrected excitation energies are systematically improved upon their static analogs, especially for singlet states. (In the case of triplets, one would notice a few cases where the excitation energies is underestimated.) In particular, the MAE is reduced from 0.64 to 0.50 eV for singlets, and from 0.41 to 0.27 eV for triplets. The MSE and RMSE are also systematically improved when one takes into account dynamical effects. The second important observation extracted from Fig. 1 is that the (singlet and triplet) Rydberg states are rather unaltered by the dynamical effects with a correction of few hundredths of eV in most cases. The same comment applies to the CT excited state of HCl. The magnitude of the dynamical correction for n → π * and π → π * transitions is much more important: 0.3-0.5 eV for singlets and 0.3-0.7 eV for triplets.

��� �

� � � � �� � � � � � � � � �� � � -��� ��� ��� ��� ��� ��� ��� ��� � � � � � �� � � � � � � � � �� � � -��� -��� ���
Dynamical BSE does not quite reach the accuracy of secondorder methods [CIS(D), ADC(2), CC2, and CCSD] for the singlet and triplet optical excitations of these small molecules. However, it is definitely an improvement in terms of performances as compared to static BSE, especially for triplet states, where dynamical BSE reaches an accuracy close to CIS(D), ADC(2), and CC2.

Table IV reports singlet and triplet excitation energies for larger molecules (acrolein H 2 C --CH -CH --O, butadiene H 2 C --CH -CH --CH 2 , diacetylene HC ---C -C ---CH, glyoxal O --CH -CH --O, and streptocyanine-C1 H 2 N -CH --NH 2 + ) at the static and dynamic BSE levels with the aug-cc-pVDZ basis set. We also report the CC3 excitation energies computed in Refs. 30-32 with the same basis set. These will be our reference as they are known to be extremely accurate (0.03-0.04 eV from the TBEs). [30][31][32][START_REF] Loos | [END_REF] Errors associated with these excitation energies (with respect to CC3) are represented in Fig. 2. As expected the static BSE excitation energies are much more accurate for these larger molecules with a MAE of 0.32 eV, a MSE of 0.30 eV, and a RMSE of 0.38 eV. Here again, the dynamical correction improves the accuracy of BSE by lowering the MAE, MSE, and RMSE to 0.23, 0.00, and 0.29 eV, respectively. Rydberg states are again very slightly affected by dynamical effects, while the dynamical corrections associated with the n → π * and π → π * transitions are much larger and of the same magnitude (0.3-0.6 eV) for both types of transitions. This latter observation is somehow different from the outcomes reached by Rohlfing and coworkers in pre- vious works 46,47 (see Sec. I) where they observed i) smaller corrections, and ii) that n → π * transitions are more affected by the dynamical screening than π → π * transitions. The larger size of the molecules considered in Refs. 46 and 47 may play a role on the magnitude of the corrections, even though we do not observe here a significant reduction going from small systems (N 2 , CO, . . . ) to larger ones (acrolein, butadiene, . . . ). We emphasize further that previous calculations 46,47 were performed within the plasmon-pole approximation for modeling the dynamical behaviour of the screened Coulomb potential, while we go beyond this approximation in the present study [see Eq. ( 23)]. Finally, while errors were defined with respect to experimental data in Refs. 46 and 47, we consider here as reference high-level CC calculations performed with the very same geometries and basis sets than our BSE calculations. As pointed out in previous works, [START_REF] Loos | [END_REF][84][85] a direct comparison between theoretical transition energies and experimental data is a delicate task, as many factors (such as zero-point vibrational energies and geometrical relaxation) must be taken into account for fair comparisons. Further investigations are required to better evaluate the impact of these considerations on the influence of dynamical screening. To provide further insight into the magnitude of the dynamical correction to valence, Rydberg, and CT excitations, let us consider a simple two-level systems with i = j = h and a = b = l, where (h, l) stand for HOMO and LUMO. The dynamical correction associated with the HOMO-LUMO transition reads

W stat hh,ll -W hh,ll (Ω 1 ) = -4[hh|hl][ll|hl]       1 Ω RPA hl - 1 Ω 1 hl -Ω RPA hl       , (43)
where the only RPA excitation energy, Ω RPA hl = ε l -ε h +2(hl|lh), is again the HOMO-LUMO transition, i.e., m = hl [see Eq. ( 26)]. For CT excitations with vanishing HOMO-LUMO overlap [i.e., (h|l) ≈ 0] and small excitonic binding energy, [hh|hl] ≈ 0 and [ll|hl] ≈ 0, so that one can expect the dynamical correction to be weak. Likewise, Rydberg transitions which are characterized by a delocalized LUMO state, that is, a small HOMO-LUMO overlap, are expected to undergo weak dynamical corrections. The discussion for π → π * and n → π * transitions is certainly more complex and moleculespecific symmetry arguments must be invoked to understand the magnitude of the [hh|hl] and [ll|hl] terms.

As a final comment, let us discuss the two singlet states of butadiene reported in Table IV. 31,[36][37][38][39][40][41][42] As discussed in Sec. I, these corresponds to a bright state of 1 B u symmetry with a clear single-excitation character, and a dark 1 A g state including a substantial fraction of double excitation character (roughly 30%). Although they are both of π → π * nature, they are very slightly altered by dynamical screening with corrections of -0.12 and -0.03 eV for the 1 B u and 1 A g states, respectively. The small correction on the 1 A g state might be explained by its rather diffuse nature (similar to a Rydberg states). 86

V. CONCLUSION

The BSE formalism is quickly gaining momentum in the electronic structure community thanks to its attractive computational scaling with system size and its overall accuracy for modeling single excitations of various natures in large molecular systems. It now stands as a genuine cost-effective excited-state method and is regarded as a valuable alternative to the popular TD-DFT method. However, the vast majority of the BSE calculations are performed within the static approximation in which, in complete analogy with the ubiquitous adiabatic approximation in TD-DFT, the dynamical BSE kernel is replaced by its static limit. One key consequence of this static approximation is the absence of higher excitations from the BSE optical spectrum. Following Strinati's footsteps who originally derived the dynamical BSE equations, 2,55,56 several groups have explored the BSE formalism beyond the static approximation by retaining (or reviving) the dynamical nature of the screened Coulomb potential 44,48,49 or via a perturbative approach coupled with the plasmon-pole approximation. 43,46,47,57 In the present study, we have computed exactly the dynamical screening of the Coulomb interaction within the random-phase approximation, going effectively beyond both the usual static approximation and the plasmon-pole approximation. Dynamical corrections have been calculated using a renormalized first-order perturbative correction to the static BSE excitation energies following the work of Rohlfing and coworkers. 43,46,47,57 Note that, although the present study goes beyond the static approximation of BSE, we do not recover additional excitations as the perturbative treatment accounts for dynamical effects only on excitations already present in the static limit. However, we hope to report results on a genuine dynamical approach in the near future in order to access double excitations within the BSE formalism. In order to assess the accuracy of the present scheme, we have reported a significant number of calculations for various molecular systems. Our calculations have been benchmarked against high-level CC calculations, allowing to clearly evidence the systematic improvements brought by the dynamical correction for both singlet and triplet excited states. We have found that, although n → π * and π → π * transitions are systematically red-shifted by 0.3-0.6 eV thanks to dynamical effects, their magnitude is much smaller for CT and Rydberg states.

  m ) are RPA neutral excitations and their corresponding transition vectors computed by solving the (static) linear response problem

  FIG.1. Error (in eV) with respect to the TBEs of Refs. 30-32 for singlet (top) and triplet (bottom) excitation energies of various molecules obtained with the aug-cc-pVTZ basis set computed within the static (red) and dynamic (blue) BSE formalism. CT and R stand for charge transfer and Rydberg state, respectively. See Tables II and III for raw data.

TABLE I .

 I Singlet and triplet excitation energies (in eV) of N 2 computed at the BSE@G 0 W 0 @HF level for various basis sets.

	State 1 Π g (n → π * )	Nature Val.	Ω stat S 9.90	cc-pVDZ (E GW g Ω dyn S 9.58	= 20.71 eV) ∆Ω -0.32 dyn S	Ω stat S 9.92	cc-pVTZ (E GW g Ω dyn S 9.53	= 20.21 eV) ∆Ω dyn S -0.40	Ω stat S 10.01	cc-pVQZ (E GW g Ω dyn S 9.59	= 20.05 eV) dyn ∆Ω S -0.42
	1 Σ -u (π → π * ) 1 ∆ u (π → π * ) 1 Σ + g 1 Π u 1 Σ + u 1 Π u	Val. Val. Ryd. Ryd. Ryd. Ryd.	9.70 10.37 15.67 15.00 22.88 a 23.62 a	9.37 10.05 15.50 14.79 22.73 23.51		-0.33 -0.31 -0.17 -0.21 -0.15 -0.11	9.61 10.27 15.04 14.75 19.03 19.15	9.19 9.88 14.84 14.48 18.95 19.04		-0.42 -0.39 -0.21 -0.27 -0.08 -0.11	9.69 10.34 14.72 14.80 16.78 16.93	9.25 9.93 14.43 14.59 16.71 16.85	-0.44 -0.41 -0.21 -0.29 -0.06 -0.09
	3 Σ + u (π → π * ) 3 Π g (n → π * )	Val. Val.	7.39 8.07		6.91 7.65		-0.48 -0.42	7.46 8.14	6.87 7.62		-0.59 -0.52	7.59 8.24	6.97 7.70	-0.62 -0.54
	3 ∆ u (π → π * )	Val.	8.56		8.15		-0.41	8.52	8.00		-0.52	8.62	8.07	-0.55
	3 Σ -u (π → π * )	Val.	9.70		9.37		-0.33	9.61	9.19		-0.42	9.69	9.25	-0.44
	State 1 Π g (n → π * )	Nature Val.	aug-cc-pVDZ (E GW g Ω stat S Ω dyn S 10.18 9.77	= 19.49 eV) ∆Ω dyn S -0.41	aug-cc-pVTZ (E GW g Ω dyn S 10.42 Ω stat S 9.99	= 19.20 eV) ∆Ω dyn S -0.42	aug-cc-pVQZ (E GW g Ω dyn S 10.52 Ω stat S 10.09	= 19.00 eV) dyn ∆Ω S -0.43
	1 Σ -u (π → π * ) 1 ∆ u (π → π * ) 1 Σ + g 1 Π u 1 Σ + u 1 Π u	Val. Val. Ryd. Ryd. Ryd. Ryd.	9.95 10.57 13.72 14.07 13.80 14.22		9.51 10.16 13.68 14.02 13.72 14.19		-0.44 -0.41 -0.04 -0.05 -0.08 -0.04	10.11 10.75 13.60 13.98 13.98 14.24	9.66 10.33 13.57 13.94 13.91 14.21		-0.45 -0.42 -0.03 -0.04 -0.07 -0.03	10.20 10.85 13.54 13.96 14.08 14.26	9.75 10.42 13.52 13.93 14.03 14.23	-0.45 -0.42 -0.02 -0.03 -0.06 -0.03
	3 Σ + u (π → π * ) 3 Π g (n → π * )	Val. Val.	7.75 8.42		7.12 7.88		-0.63 -0.54	8.02 8.66	7.38 8.10		-0.64 -0.56	8.12 8.75	7.48 8.20	-0.64 -0.56
	3 ∆ u (π → π * )	Val.	8.86		8.32		-0.54	9.04	8.48		-0.56	9.14	8.57	-0.56
	3 Σ -u (π → π * )	Val.	9.95		9.51		-0.44	10.11	9.66		-0.45	10.20	9.75	-0.45

a Excitation energy larger than the fundamental gap.

TABLE II .

 II Singlet excitation energies (in eV) for various molecules obtained with the aug-cc-pVTZ basis set computed at various levels of theory. CT stands for charge transfer.

						BSE@G 0 W 0 @HF			Wave function-based methods
	Mol. HCl	State 1 Π	Nature CT	E GW g 13.43	Ω stat S 8.30	Ω 8.19 dyn S	∆Ω -0.11 dyn S	Z S 1.009	CIS(D) 6.07	ADC(2) 7.97	CC2 7.96	CCSD 7.91	TBE 7.84
	H 2 O	1 B 1 (n → 3s)	Ryd.	13.58	8.09	8.00	-0.09	1.007	7.62	7.18	7.23	7.60	7.17
		1 A 2 (n → 3p)	Ryd.		9.79	9.72	-0.07	1.005	9.41	8.84	8.89	9.36	8.92
		1 A 1 (n → 3s)	Ryd.		10.42	10.35	-0.07	1.006	9.99	9.52	9.58	9.96	9.52
	N 2	1 Π g (n → π * )	Val.	19.20	10.42	9.99	-0.42	1.031	9.66	9.48	9.44	9.41	9.34
		1 Σ -u (π → π * ) 1 ∆ u (π → π * ) 1 Σ + g 1 Π u 1 Σ + u 1 Π u	Val. Val. Ryd. Ryd. Ryd. Ryd.		10.11 10.75 13.60 13.98 13.98 14.24	9.66 10.33 13.57 13.94 13.91 14.21	-0.45 -0.42 -0.03 -0.04 -0.07 -0.03	1.029 1.030 1.003 1.004 1.008 1.002	10.31 10.85 13.67 13.64 13.75 14.52	10.26 10.79 12.99 13.32 13.07 14.00	10.32 10.86 12.83 13.15 12.89 13.96	10.00 10.44 13.15 13.43 13.26 13.67	9.88 10.29 12.98 13.03 13.09 13.46
	CO	1 Π(n → π * )	Val.	16.46	9.54	9.19	-0.34	1.029	8.78	8.69	8.64	8.59	8.49
		1 Σ -(π → π * )	Val.		10.25	9.90	-0.35	1.023	10.13	10.03	10.30	9.99	9.92
		1 ∆(π → π * )	Val.		10.71	10.39	-0.32	1.023	10.41	10.30	10.60	10.12	10.06
		1 Σ +	Ryd.		11.88	11.85	-0.03	1.005	11.48	11.32	11.11	11.22	10.95
		1 Σ +	Ryd.		12.39	12.37	-0.02	1.003	11.71	11.83	11.63	11.75	11.52
		1 Π	Ryd.		12.37	12.32	-0.05	1.004	12.06	12.03	11.83	11.96	11.72
	C 2 H 2	1 Σ -u (π → π * ) 1 ∆ u (π → π * )	Val. Val.	12.28	7.37 7.74	7.05 7.46	-0.32 -0.29	1.026 1.025	7.28 7.62	7.24 7.56	7.26 7.59	7.15 7.48	7.10 7.44
	C 2 H 4	1 B 3u (π → 3s)	Ryd.	11.49	7.64	7.62	-0.03	1.004	7.35	7.34	7.29	7.42	7.39
		1 B 1u (π → π * )	Val.		8.18	8.03	-0.15	1.022	7.95	7.91	7.92	8.02	7.93
		1 B 1g (π → 3p)	Ryd.		8.29	8.26	-0.03	1.003	8.01	7.99	7.95	8.08	8.08
	CH 2 O	1 A 2 (n → π * )	Val.	12.00	5.03	4.68	-0.35	1.027	4.04	3.92	4.07	4.01	3.98
		1 B 2 (n → 3s)	Ryd.		7.87	7.85	-0.02	1.001	6.64	6.50	6.56	7.23	7.23
		1 B 2 (n → 3p)	Ryd.		8.76	8.72	-0.04	1.003	7.56	7.53	7.57	8.12	8.13
		1 A 1 (n → 3p)	Ryd.		8.85	8.84	-0.01	1.000	8.16	7.47	7.52	8.21	8.23
		1 A 2 (n → 3p)	Ryd.		8.87	8.85	-0.02	1.002	8.04	7.99	8.04	8.65	8.67
		1 B 1 (σ → π * )	Val.		10.18	9.77	-0.42	1.032	9.38	9.17	9.32	9.28	9.22
		1 A 1 (π → π * )	Val.		10.05	9.81	-0.24	1.026	9.08	9.46	9.54	9.67	9.43
	MAE				0.64	0.50			0.43	0.24	0.25	0.15	0.00
	MSE				0.64	0.48			0.14	0.02	0.03	0.14	0.00
	RMSE				0.70	0.58			0.55	0.33	0.33	0.20	0.00
	Max(+)				1.08	0.91			1.06	0.54	0.57	0.44	0.00
	Max(-)				0.20	-0.22			-1.77	-0.76	-0.71	-0.02	0.00

TABLE III .

 III Triplet excitation energies (in eV) for various molecules obtained with the aug-cc-pVTZ basis set computed at various levels of theory. Error (in eV) with respect to CC3 for singlet and triplet excitation energies of various molecules obtained with the aug-cc-pVDZ basis set computed within the static (red) and dynamic (blue) BSE formalism. R stands for Rydberg state. See TableIVfor raw data.

					BSE@G 0 W 0 @HF			Wave function-based methods
	Mol. H 2 O	State 3 B 1 (n → 3s)	Nature Ryd.	E GW g 13.58	Ω stat S 7.62	Ω 7.48 dyn S	∆Ω -0.14 dyn S	Z S 1.009	CIS(D) 7.25	ADC(2) 6.86	CC2 6.91	CCSD 7.20	TBE 6.92
		3 A 2 (n → 3p)	Ryd.		9.61	9.50	-0.11	1.007	9.24	8.72	8.77	9.20	8.91
		3 A 1 (n → 3s)	Ryd.		9.80	9.66	-0.14	1.008	9.54	9.15	9.20	9.49	9.30
	N 2	3 Σ + u (π → π * ) 3 Π g (n → π * )	Val. Val.	19.20	8.02 8.66	7.38 8.10	-0.64 -0.56	1.032 1.031	8.20 8.33	8.15 8.20	8.19 8.19	7.66 8.09	7.70 8.01
		3 ∆ u (π → π * )	Val.		9.04	8.48	-0.56	1.031	9.30	9.25	9.30	8.91	8.87
		3 Σ -u (π → π * )	Val.		10.11	9.66	-0.45	1.029	10.29	10.23	10.29	9.83	9.66
	CO	3 Π(n → π * )	Val.	16.46	6.80	6.25	-0.55	1.031	6.51	6.45	6.42	6.36	6.28
		3 Σ + (π → π * )	Val.		8.56	8.06	-0.50	1.025	8.63	8.54	8.72	8.34	8.45
		3 ∆(π → π * )	Val.		9.39	8.96	-0.43	1.024	9.44	9.33	9.56	9.23	9.27
		3 Σ -u (π → π * ) 3 Σ + u	Val. Ryd.		10.25 11.17	9.90 11.07	-0.35 -0.10	1.023 1.008	10.10 10.98	10.01 10.83	10.27 10.60	9.81 10.71	9.80 10.47
	C 2 H 2	3 Σ + u (π → π * ) 3 ∆ u (π → π * )	Val. Val.	12.28	5.83 6.64	5.32 6.23	-0.51 -0.41	1.031 1.028	5.79 6.62	5.75 6.57	5.76 6.60	5.45 6.41	5.53 6.40
		3 Σ -u (π → π * )	Val.		7.37	7.05	-0.32	1.026	7.31	7.27	7.29	7.12	7.08
	C 2 H 4	3 B 1u (π → π * )	Val.	11.49	4.95	4.49	-0.46	1.032	4.62	4.59	4.59	4.46	4.54
		3 B 3u (π → 3s)	Ryd.		7.46	7.42	-0.04	1.004	7.26	7.23	7.19	7.29	7.23
		3 B 1g (π → 3p)	Ryd.		8.23	8.19	-0.04	1.004	7.97	7.95	7.91	8.03	7.98
	CH 2 O	3 A 2 (n → π * )	Val.	12.00	4.28	3.87	-0.40	1.027	3.58	3.46	3.59	3.56	3.58
		3 A 1 (π → π * )	Val.		6.31	5.75	-0.56	1.033	6.27	6.20	6.30	5.97	6.06
		3 B 2 (n → 3s)	Ryd.		7.60	7.56	-0.05	1.002	6.66	6.39	6.44	7.08	7.06
	MAE				0.41	0.27			0.27	0.21	0.24	0.10	0.00
	MSE				0.41	0.06			0.23	0.10	0.14	0.05	0.00
	RMSE				0.45	0.33			0.31	0.27	0.30	0.13	0.00
	Max(+)				0.70	0.60			0.63	0.57	0.63	0.29	0.00
	Max(-)				0.11	-0.39			-0.40	-0.67	-0.62	-0.11	0.00
	���	��������
	���									�������
	���											��������������
	���				���������			�����������
	���
	���
	-���
	FIG. 2.

TABLE IV .

 IV Singlet and triplet excitation energies (in eV) for various molecules obtained with the aug-cc-pVDZ basis set computed at various levels of theory.

			BSE@G 0 W 0 @HF
	Molecule acrolein	State 1 A (n → π * ) Val. 11.67 4.62 4.28 -0.35 Nature E GW g Ω stat S Ω dyn S ∆Ω dyn S	Z S 1.030 3.77 CC3
		1 A (n → π * ) Val.	6.86 6.70 -0.16	1.023 6.67
		1 A (n → 3s) Ryd.	7.57 7.53 -0.04	1.004 6.99
		3 A (n → π * ) Val.	3.97 3.54 -0.43	1.031 3.47
		3 A (π → π * ) Val.	4.03 3.61 -0.42	1.032 3.95
	butadiene	1 B u (π → π * ) Val.	9.88 6.25 6.13 -0.12	1.019 6.25
		1 A g (π → π * ) Val.	6.88 6.86 -0.03	1.003 6.68
		3 B u (π → π * ) Val.	3.68 3.25 -0.43	1.032 3.36
		3 A g (π → π * ) Val.	5.51 5.01 -0.50	1.040 5.21
		3 B g (π → 3s) Ryd.	6.29 6.25 -0.04	1.005 6.20
	diacetylene	1 Σ -u (π → π * ) Val. 11.01 5.62 5.35 -0.28 1 ∆ u (π → π * ) Val. 5.87 5.63 -0.25	1.025 5.44 1.024 5.69
		3 Σ + u (π → π * ) Val. 3 ∆ u (π → π * ) Val.	4.30 3.82 -0.49 5.04 4.68 -0.36	1.031 4.06 1.027 4.86
	glyoxal	1 A u (n → π * ) Val. 10.90 3.46 3.14 -0.33	1.028 2.90
		1 B g (n → π * ) Val.	4.96 4.55 -0.41	1.034 4.30
		1 B u (n → 3p) Ryd.	7.90 7.86 -0.04	1.004 7.55
		3 A u (n → π * ) Val.	2.77 2.38 -0.39	1.028 2.49
		3 B g (n → π * ) Val.	4.23 3.75 -0.48	1.034 3.91
		3 B u (π → π * ) Val.	5.01 4.47 -0.55	1.034 5.20
	streptocyanine 1 B 2 (π → π * ) Val. 13.79 7.66 7.51 -0.15	1.019 7.14
	MAE		0.32 0.23	0.00
	MSE		0.30 0.00	0.00
	RMSE		0.38 0.29	0.00
	Max(+)		0.85 0.54	0.00
	Max(-)		-0.19 -0.73	0.00

E. E. Salpeter and H. A. Bethe, Phys. Rev.84, 1232 (1951).

G. Strinati, Riv. NuovoCimento 11, 1 (1988).

L. Hedin, Phys. Rev. 139, A796 (1965).

D. Golze, M. Dvorak, and P. Rinke, Front. Chem. 7, 377 (2019).

ACKNOWLEDGMENTS

The authors would like to thank Elisa Rebolini, Pina Romaniello, Arjan Berger, and Julien Toulouse for insightful discussions on dynamical kernels. PFL thanks the European

DATA AVAILABILITY

The data that support the findings of this study are available within the article.

Appendix A: Fourier transform of L 0 (1, 4; 1 , 3)

In this Appendix, we derive Eqs. (15) to (17). Combining the Fourier transform (with respect to t 1 ) of L 0 (1, 4; 1 , 3) [L 0 ](x 1 , 4; x 1 , 3 | ω 1 ) = -i dt 1 e iω 1 t 1 G(1, 3)G(4, 1 ), (A1) (where t 1 = t + 1 ) with the inverse Fourier transform of the Green's function, e.g.,

(where τ 13 = t 1t 3 ), we obtain

Applying the change of variable ω → ω + ω 1 /2, one gets

with τ 34 = t 3t 4 and t 34 = (t 3 + t 4 )/2. Finally, using the Lehman representation of the Green's functions [see Eq. ( 16)], and picking up the poles associated with the occupied (virtual) states in the upper (lower) half-plane for τ 34 > 0 (τ 34 < 0), one obtains, using the residue theorem,

with τ = τ 34 , and where pp and hh label the particle-particle and hole-hole channels (respectively) that are neglected here. 2 Projecting onto φ * a (x 1 )φ i (x 1 ) selects the first line of the righthand-side of Eq. (A5), yielding Eq. ( 17) with ω 1 → Ω s .

Appendix B: N|T [ ψ(6) ψ † (5)]|N, S in the electron-hole product basis

We now derive in more details Eq. ( 18 (B2)

Expanding now the field operators with creation/destruction operators in the orbital basis, i.e., -θ(-τ 65 ) N|â † q âp |N, S e -iε q τ 65 e -iΩ S t 6 ], (B6) leading to Eq. (18) with Ω S = E N S -E N 0 , t 6 = τ 65 /2 + t 65 , and t 5 = -τ 65 /2 + t 65 .


Sex Ratio Imbalances Among Children At Micro-Level: China And India Compared.

Christophe Z Guilmoto and S. Oliveau Guilmoto@ird.fr; sebastien.oliveau@univ-provence.fr

Sex ratio in Asia

Sex ratio imbalances found in Asia are primarily due to the increasing proportions of sons among children in China and India. 1 This gradual demographic masculinization over the last two decades results from the rapid progression in the number of sex-selective abortions aimed at ensuring the birth of male children. While other factors such as female infanticide or excess child mortality have played a sizeable role in the past in distorting sex ratios among the child population, the rise of pre-birth sex selection has entirely changed the situation. Sex selection through abortion -following detection of sex of the foetus-represents a technological breakthrough that has enabled millions of Asian couples to implement a new and more effective strategy of sex planning aimed at getting rid of unwanted female births. Figure 1 brings together child sex ratio estimates for China, India and the rest of the world (world without Asia). Sex ratio values (number of boys per hundred girls) around 107 have long been greater than normal in China, but the real increase dates from the 1980s. This change occurred at about the same period in India, although the magnitude of the rise in child ratio is much lower than in China. One reason for the slower increase recorded in India lies in the very uneven distribution of gender discrimination within the country. But as data for Punjab State in India indicates, regional values in India may be significantly higher than the Chinese average. In both Asian countries, the child sex ratio is probably still on the increase although lack of detailed annual series rules out any firm conclusion on most recent trends. During the last 40 years, the overall figure for the rest of the world has remained fairly stable at around 103, displaying almost no visible trend in sex ratio levels. 30 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 China India Punjab World without Asia 1 See [START_REF] Croll | Endangered Daughters: Discrimination and Development in Asia[END_REF]; Attané and Véron, (2005); [START_REF] Bhat | On the Trail of 'Missing' Indian Females[END_REF] 2/18

This paper offers a comparative study of sex ratio distribution in China and India. 2 We concentrate here on the measurement of sex ratio among the child population. Sex ratio imbalances among other, older age groups are more complex to interpret as they tend to be influenced by a host of other factors such as past trends in child sex ratio, sex-selective underregistration of young adults, adult mortality and sex-selective mobility behaviour [START_REF] Bhat | On the Trail of 'Missing' Indian Females[END_REF].

Mapping child sex ratio

One puzzling aspect of sex ratio imbalances among children relates to its social and geographical distribution in both countries. Whenever data are available for specific household sub-populations, it appears that sex ratio imbalances are more acute visible among affluent, urban groups while being often negligible among underprivileged groups. To some extent, there appears also to be broad similarities between fertility change and gender discrimination. However, this relationship is far from straightforward as examples abound of groups combining both low (or high) fertility and gender discrimination levels (see also [START_REF] Bhat | Fertility Decline and Gender Bias in Northern India[END_REF]. But lack of age and sex data for household populations -that would allow us to estimate child sex ratio for populations classified by socio-economic status-often limits our capacity to assess the correspondence between class membership and discriminatory behaviour towards female births.

Geographical variations are on the contrary easier to highlight thanks to the availability of age and sex pyramids for various administrative levels such provinces/states, urban agglomerations, and lower-level regional units. Census data from China and India indicate for example that gender discrimination levels may be both extremely pronounced in several regions of East China and North-west India whereas the rest of these countries often exhibit normal sex ratio levels below 105. Punjab (Figure 1) where child sex ratio is presently greater than 125 illustrate the variations that may be found within a single country (see also [START_REF] Arokiasamy | Regional patterns of sex bias and excess female child mortality in India[END_REF].

We start with detailed maps of both China and India based on results from the last census (see below for description of the data). The Chinese map (Figure 2) shows the great West-East division of the country in terms of child sex ratio. While its values tend to be usually normal in the western half of the country, the eastern half is characterized by high sex ratio levels. In East China, spatial patterns tend also to be more intricate. There are obviously several "hot spots" of masculinity in a few areas where values can at times cross the 150 mark. These areas are usually located away from the coast (except in South Guangdong and Hainan) as can be seen in Jiangxi or in Henan. The resulting picture remains fairly fragmented since hot spots may often be contiguous to areas where child sex ratio is at a normal level. We see several such instances of heterogeneity within Hubei, Anhui or Guangxi provinces. India's map (Figure 3) tends to be much easier to decipher than China's in spite of being based on a much larger number of units (see below for details). The vast majority of excess sex ratio values are shown to be clustering towards the Northwest. The highest values are recorded in the states of Punjab and Haryana, while other "core areas" may be found in Gujarat or Maharashtra. Interestingly enough, there appears to be a corridor of low child sex ratio separating North-western India from the contiguous state of Maharashtra where values are also high. This land strip closely corresponds to tribal settlements found at the borders of four States (viz. Rajasthan, Gujarat, Maharashtra and Madhya Pradesh) along the Narmada River basin. One may also distinguish in the map a few isolated hot spots such as in Tamil Nadu (Salem area) or in Orissa. The difference in spatial patterning observed in China and India may seem rather perplexing. China is to a large extent socially and culturally more homogeneous than India. Religious and linguistic differences tend for instance to be significantly more acute across regions in India than in East China which is overwhelmingly Han-populated (with standard Mandarin widely understood). Consequently, we would expect demographic behaviour to be more homogenously distributed in China. Mapping sex ratio at a fine level shows on the contrary a large amount of local disparities within Chinese provinces. We may therefore wonder whether the fragmented picture of China's sex ratio differentials is actually linked to local factors. Among possible such factors, local variations in family policy enforcement may explain why some localities would have very different demographic outcomes than their neighbours.

Another hypothesis relates to differentials in under-registration level caused by underreporting of girls either by local population or by government officials.

There may however be another explanations for local irregularities. This is best illustrated by India's map: in it, we see that Andhra Pradesh's distribution of sex ratio values is rather erratic in stern contrast to the rather regular geographical repartition observed elsewhere. Similar observations could be made for other states such West Bengal or Arunachal Pradesh. As we will show below, this irregular distribution is not caused by demographic behaviour, but rather by the insufficient demographic size of administrative units that have been mapped. This means that plotting values for small units as was done in the two previous maps may result in artificial statistical distortion when estimating local child sex ratio. Such local heterogeneity visible in some areas would then disappear when maps are based on larger units that tend to blur out local irregularity. Another consequence is that comparing China's sex ratio map with India's may not be advisable as long as administrative units for both countries are too dissimilar in both population size and area.

Measuring sex ratio

In this section, we will go back to our original datasets and describe them in greater detail. The last census was conducted respectively in 2000 in China and 2001 in India. No other data set provides an exhaustive sex and age count for both national populations. While China follows standard age classification (below 1 year, 1-4 years, 5-9 years etc.), India's data are not always available in the customary quinquennial age format. For some units, figures given by the Census relate to the populations aged "less than 7 years" or "7 years and more" as the latter age group is used to compute literacy rates. But we may add that because of the relatively high level of age heaping in India -also associated to sex of the children such as for exact age 5-, the 0-6 year age group tends to offer a more robust measurement of sex ratio than the usual 0-4 group. Another advantage relates to the larger size of the 0-6 age segment used to compute sex ratio (more on sex ratio variability below). China's age data are however on the whole of much better quality than their Indian counterparts. The difference in age group between countries may however have little impact in our comparison since adding ages 5 and 6 is unlikely to alter the overall sex ratio among children.

Child sex ratio itself happens to be a rather synthetic measure of gender discrimination as it incorporates the impact of its three main components: skewed sex ratio at birth, infanticide and sex differentials in infant and child mortality. At the same time, it is not known to be significantly affected by differentials in migration that may deeply alter sex ratio among older age groups. However, computation of sex ratio necessitates relatively large population sizes especially because it is a ratio (with both fluctuating denominator and numerator) and not a rate. A simple sensitivity analysis using the binomial distribution shows for instance that sex ratio fluctuates widely for smaller populations: for a population of 10,000 (and assuming a theoretical sex ratio of 105), the 95% confidence interval for the sex ratio is [99.3-111]. This militates against extending the geographical analysis to very small administrative units as random variations related to sample size may result in artificially high levels of local variability of sex ratio estimates.

Data from these censuses are available at various administrative levels in both countries. In China, the larger units are the 33 province-level units (provinces, municipalities and other regions), further divided into 333 prefectures while the finer units correspond to counties (Xiàn). There were 2391 counties in China at the time of the 2000 census out of which 2368 spatial units could be retained in this study. 3 In India, sex ratio data have long been examined on the basis of the 35 States and Union Territories or of the 593 districts existing during the 2001 census. This smaller number of districts in India, compared to China's counties, means that district units tend to be much larger in terms of population and are not directly comparable to China's counties.

Recently, data for so-called subdistrict units have been made available (CensusInfo 2005) for India. These units are identified in most States as "tahsils" although there are many other regional appellations such as "taluks" (in several south Indian states), "mandals" (Andhra Pradesh), "CD-blocks" in Bihar or West Bengal, "circles" in North-eastern India, etc. The main advantage of this subdistrict series that include the population below 7 is that they correspond to 5564 units, almost ten times more than districts. This grid has been used for the map of child sex ratio shown at the beginning of the paper. This administrative level combines however different regional administrative grids that tend to vary widely in area across Indian states. For instance, the 1125 mandals found in Andhra Pradesh are extremely small units. But the existence of such data (which have been further georeferenced by the authors using Census of India sources) will provide us the best dataset to confront local sex ratio variability in both China and India.

Geographical and administrative differences between China and India

Geographically, China and India present much dissimilarity. On the whole, the spatial distribution of the population tends to be highly skewed in China towards the Eastern part of the country where human density is extremely high. Large parts of China's western provinces (such as Xinjiang or Tibet) display on the contrary very low population density levels. They mostly correspond to barren, arid tracts where cultivable land represents less than 5% of the overall land surface. As a result, counties in Western China are often gigantic though possessing tiny populations. In 156 counties, the area exceeds 10,000 sq. km. Almost 400 counties have a child population that is below 10,000. Figures 4 and5 indicate that these counties are mostly located to the West of the country. Data shown on Table 1 show that Chinese counties are on average much larger in area than their Indian counterparts, but also that they are characterized by a sizeable standard deviation resulting from differences between the western and eastern parts of China. A further source of complication stems from the fact that these counties with large land area but limited population are also characterized by the lowest sex ratio observed in the country. Not only are minority groups often more numerous than ethnic Hans in these parts of the country, but fertility levels are on the whole higher than in the rest of the country. These areas are also predominantly rural and underdeveloped. These various features explain why child sex ratio in these low-density regions is at its lowest, but the combination of spatial and demographic features may also be seen as a potential source of trouble for spatial comparison. While not absent, variations in density levels across regions in India are far less pronounced than In China. Moreover, they do not correspond in any way to variations in child sex ratio levels and density heterogeneity is therefore unlikely to impact on variations in sex ratio as may be the case for China. However, the problem of small units is more acute for India where subdistricts often have limited population. As maps on Figures 6 and 7 illustrate, small and underpopulated subdistrict units tend to coincide in India, being clustered in Andhra Pradesh, Bihar or in North-East States such as Arunachal Pradesh or Nagaland. No less than 1655 subdistricts in India have a child population below 10,000 and they account for 30% of the overall number of subdistrict units. An even larger proportion of subdistricts have an area smaller than 500 sq. km. Why do these differences in administrative structures matter? The first reason relates to the above-mentioned variability of sex ratio. If computed on small populations (such as localities with less than 10,000 children), sex ratio may take values affected by sizeable random errors.

Another less obvious reason corresponds to the geographical processes that are captured in spatial analysis. In the geostatistical investigations that will follow, we will examine how far sex ratio observed in different units are affected by sex ratio values observed in other neighbouring administrative units. We will thus examine how far distance impacts on sex ratio variations. But distance across units is closely related to their average size: large regions are spatially more distant from each other than smaller ones lying closing to each other. If we want to examine the impact distance has on sex ratio variations within countries, we should therefore avoid using heterogeneous administrative grids combining very large and small units. This is especially true when the average area of administrative units tends to be correlated with the feature studied (here, child sex ratio). Moreover, to compare geographical distribution in both countries suggest using comparable units and India's are obviously too small in several regions.

The administrative grids in China and India are therefore too dissimilar to allow for direct comparison and some units appear also to be either too large or too small:

• Some units (mostly in China) that may be too large as well as underpopulated

• Some units (mostly in India) that may be too small as well as underpopulated

These problems call for different solutions. About China's large underpopulated countries, the only way was to remove from our study the provinces where they are found. We decided therefore to restrict our study to Central and Eastern China and remove the provinces of Xinjiang, Xizang (Tibet), Qinghai, Gansu and Nei Menggu (Inner Mongolia). 4 These provinces account for no more than 6% of China's overall population. Provinces retained in our analysis correspond roughly to the eighteen historical Chinese provinces (Qing Dynasty), to which have been added the northeastern provinces.

As for administrative units where population was deemed too small, we used spatial analysis tools to re-aggregate units using spatial proximity. The basic idea consists in merging counties that were close to each other and create a new spatial dataset made of slightly larger units. To do this, we have used Voronoi polygons to club together units according to their proximity. This method not only tends to reduce the number of small, underpopulated units where sex ratio computation may be fraught with estimation uncertainty, but it also produces new datasets for China and India that tend to be more similar in terms of average area.

After a few trials, it was decided to merge units that were less than 30 kilometres apart. This distance range has effectively reduced the number of underpopulated units from our original sample. Results shown on Table 1 above indicate that our new datasets for China and India are now of better quality and much more homogenous across countries. Because of spatial reaggregation procedure followed, the number of small units has drastically declined in both countries. This reflects on the number of units with child population below 10,000: they have almost disappeared in India while their proportion has greatly declined in China. We may however note that even after re-aggregation, there subsist a few units that are oversized or underpopulated in China and this is mostly caused by the presence of units located along the western border.

10/18

The number of units is now similar for both countries (1400 for Eastern China and 1512 for India) and the average area of these re-aggregated units is comparable (3,096 sq. km for China as against 2,142 for India). This procedure allows us now to compare more confidently the spatial distribution in both Asian countries using a new set of map based on re-aggregated units (Figures 8 and9). The new map of East and central China's sex ratio (Figure 8) is now easier to interpret as many local differentials have vanished. However, the patchy distribution of sex ratio values remains a major trait of China's map as no clear-cut geographical patterns has emerged. The new map drawn for India (Figure 9) has suppressed most local irregularities. This can be seen from Andhra Pradesh where sex ratio values are now uniformly below average except for a small pocket towards the South that was previously invisible. Similarly, the picture from States such as Bihar and West Bengal are now easier to interpret. The tribal corridor hinted at earlier that separates North-western India from Maharashtra remains clearly distinguishable.

Local autocorrelation

To go beyond the mere visual inspection of maps, we will resort to geostatistical tools to identify more systematically masculinity hot spots. Using Anselin's definition of local indicators of spatial association, we now submit our maps to the examination of local patterns. Each unit is now tested for local spatial autocorrelation [START_REF] Anselin | Local indicators of spatial association -LISA[END_REF][START_REF] Anselin | GeoDa 0.9 User's Guide[END_REF]. 5We start with a comparison between observed values and values obtained in neighbouring units. The procedure is also known as Moran's scatterplots. To do that, we'll use two different definitions of "neighbourhood". The first definition limits the notion of neighbours to that of immediately spatially adjacent units. The second, more comprehensive treatment of neighbourhood rests on third-order contiguity: units need to have less than 3 units between themselves to qualify as "neighbours". The resulting plots shown on Figures 10 and11 show local values (X-axis) to be strongly correlated to that of their neighbours (Y-axis). The correlation is as expected stronger for a restrictive definition of neighbourhood based only on adjacent units. Comparing now China to India, we see almost similar patterns of correlation at least using first-order contiguity. Moran's I (here computed following Anselin's formula), which corresponds to the fit of the scatterplots, is very high for both countries and at almost similar levels (.693 for India vs. .646 for China). In other words, in both China and India, there is a strong level of spatial dependency as far as gender discrimination is concerned: the correlation between adjacent units means that strictly local variations in sex ratio levels are uncommon. This of course calls to mind the possible role of diffusion processes to explain how local gender attitudes tend to be so homogeneous on this scale.

The difference between China and India appears more significant when we use a higher-order contiguity definition for neighbourhood. This implies that spatial autocorrelation tends to decrease faster in China and that distance has a stronger impact on heterogeneity between administrative units. However, Moran's values remain strong and significant even for distant neighbours.

Using scatterplots, we may now define hot spots as areas where child sex ratio values are the highest with neighbouring values also high ("high-high", in red on figure 12 & 13). In the same vein, "cold spots" can also be identified where sex ratio is below average in both individual units and their neighbours ("low-low", in blue). A more generous definition of neighbourhood -based on third-order contiguity-logically results in a larger number of hot spots. In India we distinguish the two core regions previously mentioned: North West India as well as inner Maharashtra. Other areas with high sex ratio tend to be too small to show up on the map of hot spots. A very large region extending from Kerala to the Southwest towards Northeast India is on the contrary characterized by homogenous, lower than average sex ratios. In China, the largest sex ratios belong to a large strip stretching from Hainan in the South to Henan in Yellow River Valley. This zone remains however slightly fragmented with lower than expected sex ratio (or spatial outliers) found in many interstices. Apart from parts lying along the Western border of our China map, there are several regions with homogeneous low sex ratio values such as the North-eastern states. Another visible such region is found in one of China's most modernized area, i.e. the Yangtze River delta in Shanghai and Zhejiang provinces.

Discussion

This paper has attempted to examine the spatial patterns of gender discrimination in China and India using geostatiscal tools that are becoming more common in demography [START_REF] Voss | Explorations in Spatial Demography[END_REF]). Spatial distributions of child sex ratio in both countries are rather dissimilar. India's patterns are oriented towards the Northwest whereas East China is characterized by a larger number of distinct hot spots. These differences may be related to the process through which gender discrimination has emerged in China and India. At the same time, it is also possible that such spatial variations between China and India could be partly artificial, resulting from the specific characteristics of each national administrative grid.

To examine the latter hypothesis, this paper has sought to harmonize the administrative grid in both countries. To do that, we have used for the first time India's subdistrict map that offers a much finer level of detail than district-based maps that are typically used. We have then conducted a re-aggregation procedure for both countries, after removing China's western provinces where population is too sparse. Results presented above are therefore less likely to be influenced by specific national administrative grids. Moreover, since the units used are of roughly comparable characteristics, the comparison of spatial autocorrelation across countries is now less risky.

What emerges is that local spatial autocorrelation at short distances is almost as strong in China as in India. In both countries, local sex ratio values tend to be closely correlated to values observed in neighbouring localities. This geographical clustering is a typical feature of sex ratio degradation in Asia and directly refers to the diffusion process of discriminatory behaviour. While it may be in theory conceivable that such spatial patterning proceeds from the spatial distribution of other characteristics, previous analysis has failed to identify any likely determinant that would explain away this peculiar spatial distribution.

In India's case, for instance, spatial autocorrelation has regularly been shown to be an irreducible component of regional variations in sex ratio levels [START_REF] Kishor | May God Give Sons to All": Gender and Child Mortality in India[END_REF][START_REF] Murthi | Mortality, Fertility, and Gender Bias in India: A District-Level Analysis[END_REF][START_REF] Malhotra | Fertility, Dimensions of Patriarchy, and Development in India[END_REF]). To take one example related to the religious factor, high sex ratio in Punjab is often related to the spatial distribution of the Sikh population among which female discrimination is intense. But by the same token, maps demonstrate that child sex ratio in neighbouring Haryana state, which is strongly Hindu-dominated, is equally as high.

Likewise, child sex ratio among Muslim communities, which as a whole records low sex ratio values in India, tends to be significantly higher in the Northeast than elsewhere. Another important feature of child sex ratio in India is that it turns out to be the most spatially autocorrelated variable [START_REF] Oliveau | Spatial autocorrelation and demography. Exploring India's demographic patterns[END_REF], which suggests that spatial patterns observed for child sex ratio are unlikely to be the sheer corollary of spatial patterns of some other factor.

Correlates of high sex ratio in China are often far less established than for India. Within Handominated areas, incriminated factors such as low fertility, family patterns, urbanization or socio-economic status seem to play but a minor role if any [START_REF] Lavely | Spatial Variation of Juvenile Sex Ratios In the 2000 Census of China[END_REF]. This leaves spatial clustering partly unexplained by social or economic factors and regression residuals remain spatially autocorrelated. Differential registration levels at the local level may be a hidden factor, but there is somewhat conflicting evidence on this aspect.6 Moreover, the rather strong spatial regularity of sex ratio differentials within East China in 1990 and 2000 suggests that local factors are not spurious. Moreover, strong spatial autocorrelation at the local level indicates that these other factors -such as variations in reporting levels-would exhibit spatial autocorrelation, a feature that may sound rather unlikely. Variations in the quality of registration tend indeed to be rather random and not permanent over time. Other unobserved factors may be ultimately to sub-ethnic value systems with Han China on which we have little testable information. There seems to be however no direct connection between Chinese dialect geography and sex ratio distribution.

These observations reinforce the hypothesis that spatial patterns closely reflect the mechanisms behind the propagation of discriminatory behaviour in both India and China during the last three decades. Several factors often thought to fuel discrimination such as availability of new technologies (scanning and abortion facilities) would tend to level off spatial differentials or to determine on the contrary a strong urban bias as technology diffusion proceeds in a top-down manner along the urban hierarchy. But no such urban bias can be detected from our maps. Fertility decline is also poorly correlated to local variations in sex ratio differentials and the correlation is at times presumed to be curvilinear [START_REF] Lavely | Spatial Variation of Juvenile Sex Ratios In the 2000 Census of China[END_REF][START_REF] Legge | Morality Policy and Unintended Consequences: China's "One-Child" Policy[END_REF], a feature most likely to disturb modelling attempts. Several instances of early declining regions such as Kerala in India or Shanghai in China do not indeed exhibit any substantial levels of discrimination towards girls.

The role of cultural macro-regions from where discriminatory practices seem to originate is a much more plausible factor behind the observed spatial distribution in India. Discriminatory attitudes and technology has diffused from core regions towards culturally homogeneous, adjacent regions over the years. New hot spots that have recently emerged such as in coastal Orissa are therefore likely to expand mechanically towards neighbouring regions. There is an unmistakable aspect of laissez-faire in the way demographic outcomes evolve in India, more often guided by changing level of demand than by supply factors or government intervention. But in China, these regions appear to be smaller in size and their impact on neighbouring areas appears more limited than in India. The somewhat fragmented pattern obtained in the latter country probably relates to the strong capacities of local and regional authorities to influence demographic outcomes. Local enforcement of regulations of family planning or sexselective abortion, not to mention local abilities to distort census data, do not however extend over distant zones and proceeds from entirely different mechanisms than unrestrained diffusion processes across social groups and localities as observed in India. Wide market mechanisms in India vs. local official regulations in China may therefore explain to a large extent the differences in spatial patterns in child sex ratio observed between both countries.

Figure 1 :

 1 Figure 1: Child sex ratio in China and India, 1950-2000

Figure 2 :

 2 Figure 2: Child sex ratio in China, counties, 2000.

Figure 3 :

 3 Figure 3: Child sex ratio in India, subdistricts, 2001.

Figure 5 :

 5 Figure 4: Chinese counties with area above 5000 sq. km

Figure 6 :

 6 Figure 6: Indian subdistricts with area below 500 sq. km.

Figure 7 :

 7 Figure 7: Indian subdistricts with child population below 10,000.

Figure 8 :

 8 Figure 8: Child sex ratio in China, re-aggregated counties, 2000.

Figure 9 :

 9 Figure 9: Child sex ratio in India, re-aggregated subdistricts, 2001.

Figure 10 :

 10 Figure 10: Moran scatterplots for China and India, first-order contiguity.

Figure 11 :

 11 Figure 11: Moran scatterplots for China and India, third-order contiguity.

Figure 12 :

 12 Figure 12: Child sex ratio hot spots in China, third-order contiguity.

Figure 13 :

 13 Figure 13: Child sex ratio hot spots in India, third-order contiguity.

Country Name of lower- level units Number Number of units used for this study Area (sq. km) Number of units with area Child population Number of units with less than 10000 children >10000 sq. km >5000 sq. km <500 sq. km

	China	Counties (Xiàn)	2,391	2,368	Av: 3,807	156	308	400	Av: 29,093	398
					Sd: 11,473				Sd: 25,914
	India	Subdistricts	5,564	5,454	Av: 579	2	14	3290	Av : 28,806	1,655
		(Tehsils, taluks, mandals, etc.)			Sd: 1,314				Sd: 36,741
	China	Re-aggregated		1400	Av 3,096	11	102	5	Av : 45,790	127
		counties			Sd: 2,902				Sd: 42,893
	India	Re-aggregated		1512	Av : 2,142	3	17	8	Av : 105,134	64
		subdistricts			Sd: 2,391				Sd: 117,224
	Notes: Av: average value; Sd: standard deviation.

Table 1 : Characteristics of lower-level administrative units, China (2000) and India (2001)

 1

For other comparisons between China and India, see Das[START_REF] Das Gupta | Why is son preference so persistent in East and South Asia?: A crosscountry study of China, India and the Republic of Korea[END_REF],[START_REF] Guilmoto | A spatial and statistical examination of child sex ratio in China and India[END_REF] and[START_REF] Guilmoto | The geography of deteriorating child sex ratio in China and India[END_REF].

Data for several lower-level units are missing for a variety of reasons in both India and China.

Taiwan is also removed for lack of data.

See also[START_REF] Bailey | Interactive Spatial Data Analysis[END_REF],[START_REF] Haining | Spatial Data Analysis: Theory and Practice[END_REF] and[START_REF] Goodchild | Spatially Integrated Social Science, coll. Spatial Information Systems[END_REF] on spatial autocorrelation measurement and interpretation.

See for instance[START_REF] Merli | Underreporting of Births and Infant deaths in Rural China: Evidence from Field Research in One County of Northern China[END_REF],[START_REF] Merli | Are Births Underreported in rural China? Manipulation of Statistical Records in response to China's Population Policies[END_REF],[START_REF] Li | Gender differences in child survival in contemporary rural China: a county study[END_REF],[START_REF] Wu | Determinants of High Sex Ratio among Newborns: A Cohort Study from Rural Anhui Province, China[END_REF] and[START_REF] Chu | Prenatal sex determination and sex-selective abortion in rural central China[END_REF].