Angel Alastuey

Riccardo Fantoni

Fourth moment sum rule for the charge correlations of a two-component classical plasma

Keywords: numbers: 05, 20, -y , 05, 20, Gg

may come

Fourth moment sum rule for the charge correlations of a two-component classical plasma

Angel Alastuey, Riccardo Fantoni

I. INTRODUCTION

Sum rules have been playing an important role in the study of charged systems for many years. In general, a sum rule provides a relation between microscopic correlations on the one hand, and macroscopic or universal quantities on the other hand. For charged systems, sum rules often express screening properties, so they shed light on the fundamental mechanisms at work. Furthermore, they also provide useful constraints for approximate theories. Sum rules have been derived for a large variety of charged systems, including classical, quantum and relativistic plasmas, while they concern both static and dynamic properties in the bulk or near interfaces. For instance, let us mention the last work [START_REF] Jancovici | Charge and current sum rules in quantum media coupled to radiation II[END_REF] by Bernard Jancovici, devoted to the study of the time-displaced charge correlations of a relativistic one-component plasma coupled to radiation. Other examples, including many contributions from Jancovici, can be found in two reviews [START_REF] Martin | Sum rules in charged fluids[END_REF][START_REF] Brydges | Coulomb systems at low density[END_REF].

One of the most well-known sum rules for classical ionic fluids was derived long ago by Stillinger and Lovett [START_REF] Stillinger | General restriction on the distribution of ions in electrolytes[END_REF], who shown that the second moment of equilibrium charge correlations is given by a simple universal expression, valid in any plasma phase and independent of the microscopic details of the considered models. That second-moment sum rule expresses the perfect screening of weak external charges. A few years later, Vieillefosse and Hansen [START_REF] Vieillefosse | Statistical mechanics of dense ionized matter. Hydrodynamic limit and transport coefficients of the classical one-component plasma[END_REF] derived another sum rule for the fourth moment of the charges correlations of the one-component plasma (OCP), where identical positively charged particles move in a rigid uniform neutralizing background. That fourth moment is expressed in terms of the isothermal compressibility. Soon after that work, there was an attempt [START_REF] Giaquinta | Hydrodynamic correlation functions for molten salts[END_REF] to extend such a fourth-moment sum rule to the two-component plasma (TCP) where both positive and negative charges are mobile. The corresponding expression for the fourth moment involves ill-defined thermodynamic quantities, so its validity remained quite doubtful. A more convincing approach for that TCP was introduced by van Beijeren and Felderhof [START_REF] Van Beijeren | Thermodynamics and correlation functions of plasmas and electrolyte solutions[END_REF].

However, the thermodynamical quantities involved in the expression of the fourth moment are defined through the application of suitable external potentials which are not explicited, while the derivation itself is rather tough. In fact, a similar expression was obtained later by Suttorp and van Wonderen [START_REF] Suttorp | Equilibrium properties of a multi-component ionic mixture: I. Sum rules for correlation functions[END_REF] for a multicomponent ionic mixture (MIM), where all mobile charges have the same sign and interact via the pure Coulomb potential, while a rigid uniform background of opposite charge ensures overall neutrality. Then, all involved thermodynamic quantities become well defined within the considered MIM.

The main goal of the present paper is to derive a fourth moment sum rule for a general TCP, namely to express such moment in terms of suitably defined thermodynamical quantities, similarly to the formulae derived for the OCP [START_REF] Vieillefosse | Statistical mechanics of dense ionized matter. Hydrodynamic limit and transport coefficients of the classical one-component plasma[END_REF] or the MIM [START_REF] Suttorp | Equilibrium properties of a multi-component ionic mixture: I. Sum rules for correlation functions[END_REF], [START_REF] Suttorp | Sum rules for correlation functions of ionic mixtures in arbitrary dimension d ≥ 2[END_REF]. Our strategy, inspired by Jancovici's style, consists in studying the response of the TCP to a weak external charge distribution with a plane wave structure. In the long wavelength limit, the induced local particle densities vary on macroscopic scales. This allows us to compute the response within some hydrostatic-like approach which involves local equilibrium states with arbitrary densities. As a crucial point, a proper definition of equilibrium homogeneous non-neutral states with arbitrary densities naturally emerges within the framework of density functional theory (DFT). Then the induced charge density is expressed in terms of well-behaved thermodynamical quantities of an auxiliary system, which is nothing but a TCP immersed in a rigid uniform neutralizing background. Comparing that expression with the general linear response formula, we obtain the required fourth moment sum rule for the charge correlations of the genuine unperturbed TCP.

According to the previous strategy, we first introduce in Section II the various systems which intervene in our analysis. Of course, we start by defining the TCP, where a shortrange regularization of the pure Coulomb interactions is essential for avoiding the classical collapse between oppositely charged particles. Two examples of such regularizations are provided, associated with either soft or hard spheres. After recalling that the TCP is always neutral and homogeneous in the absence of any external action on the particles, we show how the application of a suitable external potential produces homogeneous nonneutral states. This leads to the introduction of an auxiliary system, the TCP immersed in a charged background, for which equilibrium states are well defined for any set of particle densities.

The general framework of DFT is exposed in Section III, where we provide the fundamental DFT equation which relates particle densities to the applied external potentials.

The central object in that relation is the free energy, which is a functional of particle densities. For slow spatial variations, that free energy functional can be expanded in powers of the gradients of particle densities, where the local ingredients are equilibrium quantities of the above homogeneous auxiliary system. Let us mention that the idea of using density-gradient expansions was introduced a long time ago by van der Waals [START_REF] Van Der Waals | The thermodynamic theory of capillarity under the hypothesis of a continuous variation of density[END_REF] for studying capillarity.

In Section IV, within DFT, we compute the induced particle densities of the TCP submitted to a weak external charge distribution with a plane wave structure and wavenumber k. The resulting induced charge density exactly cancels the external charge distribution in the long wawelength limit k → 0, as expected from perfect screening arguments. Furthermore, its amplitude at the order k 2 included only depends on thermodynamical quantities of the auxiliary system. In other words the square-gradient corrections in the free energy functional, which intervene in the corresponding amplitudes of each induced particle density at this order, do not contribute anymore when forming the charge density thanks to cancellations. Then, by comparing this exact expression of the induced charge density obtained by DFT with the linear response formula, we obtain the required sum rule for the fourth moment of the charges correlations of the homogeneous neutral TCP. The corresponding thermodynamical expression of that fourth moment involves not only the isothermal compressibility of the TCP, but also partial compressibilities specific to the auxiliary system. We briefly discuss the content of previous approaches [START_REF] Giaquinta | Hydrodynamic correlation functions for molten salts[END_REF], [START_REF] Van Beijeren | Thermodynamics and correlation functions of plasmas and electrolyte solutions[END_REF], and we show how the known results for the OCP [START_REF] Vieillefosse | Statistical mechanics of dense ionized matter. Hydrodynamic limit and transport coefficients of the classical one-component plasma[END_REF] and the MIM [START_REF] Suttorp | Sum rules for correlation functions of ionic mixtures in arbitrary dimension d ≥ 2[END_REF] can be easily recovered within our general method.

It is worthy to check explicitly the fourth moment sum rule for specific models where exact calculations can be carried out for both microscopic and thermodynamical quantities. In Section V, we consider a model of charged soft spheres in the low density limit at fixed temperature. Within the Abe-Meeron resummations of the familiar Mayer diagrammatics for particle correlations, we first compute the lowest order terms in the density expansion of the fourth moment of charge correlations, namely the terms of order 1/ρ, 1/ρ 1/2 , ρ ln ρ and ρ 0 in the density ρ. Abe-Meeron resummed diagrammatics also provide the low density expansion of the thermodynamical quantities involved in the fourth moment sum rule : the corresponding expansion of the thermodynamical expression of the fourth moment exactly coincides with the previous purely microscopic calculation up to order ρ 0 included. That remarkable agreement holds for any values of the microscopic parameters defining the model.

In Section VI, we provide some additional comments about the derivation itself and its underlying assumptions, as well as extensions to three and more component systems. Beyond its own conceptual interest, we also discuss a possible use of the fourth moment sum rule for a better understanding of the conductor or dielectric nature of the critical point of the liquid-gas transition of an ionic fluid. It turns out that recent Monte Carlo simulations [START_REF] Das | When is a conductor not perfect? Sum rules fail under critical fluctuations[END_REF] strongly suggest that the fourth moment of charge correlations diverge when approaching the critical point, in a way close to that of the isothermal compressibility.

That observation was one of the motivations for the present work.

II. THE SYSTEMS OF INTEREST

A. Examples of two-component plasmas

We consider a two-component classical plasma (TCP) made with two species α = 1, 2 of mobile particles carrying positive or negative charges, let us say q 1 = Z 1 q > 0 and q 2 = -Z 2 q < 0 with Z 1 and Z 2 positive integers. The corresponding Hamiltonian for a total number of particles

N = N 1 + N 2 reads H N 1 ,N 2 = N i=1 p 2 i 2m α i + 1 2 i =j u α i α j (r i , r j ) (1)
where α i = 1, 2 is the species of particle i. The two-body potential u α i α j (r i , r j ) only depends on the relative distance r = |r ir j |, and it includes some short-range regularization of the Coulomb interaction, which is crucial for avoiding the classical collapse between oppositely charged particles. A first soft regularization is embedded in the simple expression

u αγ (r) = q α q γ r [1 -exp(-r/d αγ )] (2)
which remains finite at r = 0. The lengths d αγ control the exponential decay at large distances of the short-range part of the two-body potential.

A second regularization amounts to introduce hard cores, namely u αγ (r) = ∞ for r < σ αγ and u αγ (r) = q α q γ r for r > σ αγ .

The corresponding TCP of charged hard spheres is suitable for describing many ionic mixtures, where the hard-core interaction mimics the effective repulsion between the electronic clouds of two ions. If σ 11 and σ 22 can be viewed as the effective diameters of the ions, the characteristic crossed lengths σ 12 = σ 21 differ in general from the average (σ 11 + σ 22 )/2 which would arise if particles really were billiard balls. This so-called non-additivity can be understood by noticing that the σ αγ 's are the typical ranges of the repulsions between electronic clouds for which no pure geometrical considerations apply. The simplest version of that general asymmetric TCP is the celebrated Restrictive Primitive Model, which is fully symmetric with respect to the charges and the hard-core diameters, namely

|q 1 | = |q 2 | = q and σ 11 = σ 22 = σ 12 = σ.
Other short-range regularizations of the Coulomb interaction can be introduced. The corresponding most general TCP will be denoted S. The following derivations are valid for any S, beyond the above two examples.

B. The homogeneous neutral TCP

Let us first consider that S is enclosed in a box with volume Λ, while no external potential is applied to the particles. At equilibrium, all statistical ensembles should become equivalent in the thermodynamic limit which is also assumed to exist. Furthermore, in a fluid phase, the bulk is overall neutral, that is the homogeneous particle densities ρ 1 and ρ 2 far from the boundaries satisfy the local neutrality relation

q 1 ρ 1 + q 2 ρ 2 = 0 . (4)
Strictly speaking, these remarkable results have been only proved in the Debye regime, namely at sufficiently high temperatures and sufficiently low densities, for rather general regularized interactions and rational ratios q 2 /q 1 [START_REF] Brydges | Debye screening[END_REF]. Moreover, there exists a proof for charge symmetric systems, i.e. q 1 = -q 2 , for any values of the thermodynamic parameters [START_REF] Frohlich | Correlation inequalities and the thermodynamic limit for classical and quantum continuous systems[END_REF]. Let us also mention the beautiful proof for the quantum version with pure Coulomb interactions by Lieb and Lebowitz [START_REF] Lieb | The constitution of matter: existence of thermodynamics for systems composed of electrons and nuclei[END_REF]. According to all those rigorous results, it can be reasonably expected that both the existence of the thermodynamic limit and the local neutrality are valid for any classical TCP in the whole fluid phase.

Important features of the various statistical ensembles are associated with the neutrality relation (4). In the grand-canonical ensemble, the intensive thermodynamical parameters are the inverse temperature β and the chemical potentials µ α . It turns out that only the linear combination µ = (Z 2 µ 1 + Z 1 µ 2 )/(Z 1 + Z 2 ) is relevant and entirely determines the total particle density ρ = ρ 1 + ρ 2 . This can be readily understood within the following simple heuristic arguments. Let us introduce, for any arbitrary configuration, the total number of particles N = N 1 + N 2 and the corresponding total charge Q = Mq

with M = Z 1 N 1 -Z 2 N 2 .
According to the decomposition

µ 1 N 1 + µ 2 N 2 = µN + νM (5)
with ν = (µ 1 -µ 2 )/(Z 1 + Z 2 ), we see that µ controls the grand-canonical average < N > GC of the total particle number, while ν determines the grand-canonical average < Q > GC = q < M > GC of the net charge. In the thermodynamic limit (TL), Λ → ∞ with β and µ α fixed, the contributions of non-neutral configurations with Q proportional to the volume Λ become negligible, because the corresponding Boltzmann factors involve a positive self-electrostatic energy which diverges faster than Λ itself. Accordingly, the total charge density in the bulk

q 1 ρ 1 + q 2 ρ 2 = lim TL < Q > GC /Λ (6)
vanishes for any given ν, while the total particle density ρ = ρ 1 + ρ 2 is indeed entirely determined by µ and β.

In the canonical ensemble, the TL is defined by letting Λ → ∞ and N α → ∞, keeping β and N α /Λ fixed. All excess charges go to the the boundaries in the TL, and the remaining bulk is always neutral. The bulk thermodynamic quantities and bulk distribution functions computed within the canonical ensemble then become identical to their grand canonical counterparts. In particular the free-energy density in thermal units of this homogeneous neutral phase, which only depends on ρ and β = 1/(k B T ), can be computed through

f (ρ, β) = lim TL (β α µ α < N α > GC -ln Ξ)/Λ , ( 7
)
where Ξ is the grand-canonical partition function. This provides the familiar thermodynamic identity

f (ρ, β) = β(ρµ -P ) , (8)
with the pressure P = lim TL k B T Λ -1 ln Ξ. Since the pressure is also given by the thermodynamic relation

βP = ρ ∂f ∂ρ -f (ρ, β) , (9)
we infer

βµ = ∂f ∂ρ , (10)
which is analogous to the standard thermodynamical identity expressing the chemical potential for a one-component system with short-range interactions.

C. The homogeneous non-neutral TCP in an external potential

In order to obtain a non-neutral homogeneous state of S with arbitrary uniform densities, one must apply a non-vanishing external potential on the particles. Let us introduce the electrostatic potential ϕ B (r) created by an uniform charge density c B , and the corresponding external potentials U B α (r) = q α ϕ B (r) seen by the particles. At equilibrium, the total electrostatic field inside the bulk should identically vanish. According to that simple electrostatic argument, the induced particle densities should be homogeneous, while the resulting charge density q 1 ρ 1 + q 2 ρ 2 carried by the particles should cancel the external charge density c B . Interestingly, the above quite plausible scenario has been exactly demonstrated within a solvable model by Jancovici [START_REF] Jancovici | Exact results for the two-dimensional one-component plasma[END_REF]. He considered identical point particles in two dimensions with pure Coulomb interactions, which then take a logarithmic form. In addition the particles are submitted to a confining parabolic potential, associated with a fixed external uniform charge density. For a special value of the temperature, all equilibrium distribution functions can be exactly computed. The resulting particle density is indeed uniform in the bulk and such that the total charge density vanishes. Furthermore, all higher-order distribution functions in the bulk become translationally invariant in the TL.

D. The auxiliary system in a neutralizing rigid background

As suggested by the previous considerations, and for further purposes, it is convenient to introduce an auxiliary system S * , where now the mobile positive and negative charges of the TCP are immersed in an uniform rigid background with charge density c B . The corresponding Hamiltonian of S * reads

H * N 1 ,N 2 = N i=1 p 2 i 2m α i + 1 2 i =j u α i α j (r i , r j ) + N i=1 Λ dr q α i c B |r i -r| + 1 2 Λ 2 drdr ′ c 2 B |r ′ -r| , (11)
when the system is enclosed in a box with volume Λ. That system can be viewed as an extension of the well-known One-Component Plasma (OCP) made of identical charged particles immersed in a neutralizing rigid background. Now, there are two species which are immersed in the background, similarly to the case of a Binary Ionic Mixture (BIM).

However, notice that here we do need a short-range regularization of the Coulomb interaction in order to avoid the collapse between oppositely charged particles, while the BIM can be defined with pure 1/r Coulomb interactions because all mobile charges have the same sign.

Like the OCP or the BIM, the system S * should have a well-behaved thermodynamic limit, which is now taken with a fixed background charge density c B . Now, in the bulk region, which is again electrically neutral, the homogeneous particle densities satisfy the neutrality relation

q 1 ρ 1 + q 2 ρ 2 + c B = 0 . ( 12
)
The corresponding free-energy density f * of the homogeneous neutral system now depends on β, c B and one particle density. Equivalently, f * depends on β and on the two particle densities ρ 1 and ρ 2 . For any given set (ρ 1 , ρ 2 ), the charge background density is adjusted in order to satisfy the neutrality relation [START_REF] Brydges | Debye screening[END_REF]. This defines the function f * (ρ 1 , ρ 2 , β),

where now ρ 1 and ρ 2 are independent variables. That procedure is analogous to that which defines the free-energy density f OCP (ρ, β) of the OCP for any value of the particle density ρ where a suitable background charge density always ensure overall neutrality.

The homogeneous neutral TCP can be viewed as a particular realization of S * for densities (ρ 1 , ρ 2 ) satisfying the neutrality relation (4) in the absence of any background.

For such neutral sets, each particle density can be expressed in terms of the total particle number density as

ρ 1 = Z 2 Z 1 + Z 2 ρ and ρ 2 = Z 1 Z 1 + Z 2 ρ . (13)
The thermodynamic quantities of the homogeneous neutral TCP can then be inferred from their counterparts of S * for the neutral set of densities [START_REF] Frohlich | Correlation inequalities and the thermodynamic limit for classical and quantum continuous systems[END_REF]. For instance, the free-energy density of the homogeneous neutral TCP is given by

f (ρ, β) = f * (Z 2 ρ/(Z 1 + Z 2 ), Z 1 ρ/(Z 1 + Z 2 ), β) . (14)
For further purposes, it is useful to consider the isothermal compressibility defined by

χ T = -lim TL Λ -1 ∂Λ ∂P = ρ ∂P ∂ρ -1 = β ρ 2 ∂ 2 f ∂ρ 2 -1 , (15)
where all partial derivatives are taken at fixed β. According to identity ( 14), χ T can be recast as

χ T = β(Z 1 + Z 2 ) 2 ρ 2 [Z 2 2 ∂ 2 f * ∂ρ 2 1 + Z 2 1 ∂ 2 f * ∂ρ 2 2 + 2Z 1 Z 2 ∂ 2 f * ∂ρ 1 ∂ρ 2 ] -1 , (16)
where the second order partial derivatives of f * are evaluated at the neutral set [START_REF] Frohlich | Correlation inequalities and the thermodynamic limit for classical and quantum continuous systems[END_REF].

Eventually, all distribution functions of the homogeneous neutral TCP obviously reduce to those of S * for that set of densities.

III. DENSITY FUNCTIONAL THEORY

A. Grand-canonical description

Now we consider a general inhomogeneous state of S, where each particle of species α is submitted to an external potential U α (r). We define the inhomogeneous fugacity of each species by

z α (r) = exp[β(µ α -U α (r))] (2πλ 2 α ) 3/2 . ( 17
)
where λ α = (β 2 /m α ) 1/2 is the de Broglie thermal wavelength of species α. The classical grand-canonical partition function of S enclosed in a box with volume Λ reads

Ξ = ∞ N 1 ,N 2 =0 1 N 1 !N 2 ! N i=1 dr i z α i (r i ) exp(-βV N 1 ,N 2 ) , (18)
where V N 1 ,N 2 is the potential part of the Hamiltonian (1). The inhomogeneous particle density ρ α (r) can be expressed as a functional derivative of Ξ with respect to z α (r), namely

ρ α (r) = z α (r) δ ln Ξ δz α (r) (19)
while parameters Λ and β are kept fixed. The free-energy F in thermal units of S is given by the Legendre transformation

F = α Λ dr ρ α (r) β(µ α -U α (r)) -ln Ξ . (20)
The grand-partition function Ξ, as well as the free-energy F , can be considered as functionals of either z α (r) or ρ α (r). The functional derivative of F with respect to ρ α (r)

is readily computed as δF δρ α (r) = β(µ α -U α (r)) (21)
where we have used identity [START_REF] Bongiorno | Molecular theory of fluid interfaces[END_REF] as well as standard calculation rules for functional differentiation. The relation [START_REF] Hansen | Theory of simple liquids[END_REF] will play a key role in the following. The density profiles ρ α (r) for each given sets (µ α -U α (r)) can be determined from that relation if one knows the functional dependence of F with respect to the inhomogeneous densities : this is the strategy of density functional theories (DFT). However, the main difficulty of DFT is that the free-energy functional is not exactly known, except for hard rods in one dimension [START_REF] Percus | Equilibrium state of a classical fluid of hard rods in an external field[END_REF].

In general, approximate forms are used. Here, we will use exact asymptotic expansions for densities with infinitely slow spatial variations.

B. Homogeneous systems

As argued in the previous Section, homogeneous states of S with arbitrary densities (ρ 1 , ρ 2 ) are obtained by applying the external potentials

U B α (r) = q α ϕ B (r) = q α Λ dr ′ c B |r ′ -r| (22)
with the external charge density c B = -(q 1 ρ 1 + q 2 ρ 2 ). If we introduce the potential part

V * N 1
,N 2 of the Hamiltonian (11) for the auxiliary system S * with background charge density c B , we can rewrite

V N 1 ,N 2 + N i=1 U B α i (r i ) = V * N 1 ,N 2 -W B ( 23
)
where

W B = 1 2 Λ 2 dr dr ′ c 2 B |r ′ -r| (24)
is the self-electrostatic energy of the background charge density c B . Inserting relation [START_REF] Vieillefosse | Sum rules for the one-component plasma with additional shortrange forces[END_REF] into the general expression [START_REF] Yang | Molecular theory of surface tension[END_REF], we obtain for the grand-partition function of S submitted to the external potentials

U B α (r), Ξ{µ 1 -U B 1 (•), µ 2 -U B 2 (•)} = Ξ * {µ 1 , µ 2 } exp(βW B ) , (25)
where Ξ * is the grand-partition function S * for the same chemical potentials µ α and a background charge density c B , without any applied external potentials, i.e. U α (r) = 0.

In the derivation of identity ( 25), we have used that W B is a pure constant which does not depend on the particle degrees of freedom. This also implies that all the grandcanonical averages for S submitted to U B α (r) are identical to those for S * . In particular, both particle densities are identical, so the particle densities of S submitted to U B α (r) are indeed homogeneous and they are merely related to c B via the neutrality condition [START_REF] Brydges | Debye screening[END_REF] valid for S * . Moreover, all particle correlations of both systems are identical.

If we insert expression [START_REF] Mayer | Statistical Mechanics[END_REF] of Ξ into definition (20) of the free energy of S submitted to U B α (r), we obtain in the TL

F {ρ 1 , ρ 2 } - β 2 Λ 2 drdr ′ (q 1 ρ 1 + q 2 ρ 2 ) 2 |r ′ -r| ∼ F * {ρ 1 , ρ 2 } ( 26
)
while c B has also been replaced by -(q 1 ρ 1 + q 2 ρ 2 ) thanks to the neutrality condition [START_REF] Brydges | Debye screening[END_REF].

Thus, if we define, in general, the reduced free energy F red of S by subtracting to F the self-electrostatic energy of the charge distribution (q 1 ρ 1 (r) + q 2 ρ 2 (r)), namely

F red = F - β 2 Λ 2 drdr ′ (q 1 ρ 1 (r) + q 2 ρ 2 (r))(q 1 ρ 1 (r ′ ) + q 2 ρ 2 (r ′ )) |r ′ -r| , (27)
we find the remarkable identity

F red {ρ 1 , ρ 2 } ∼ F * {ρ 1 , ρ 2 } , (28)
which can be rewritten for the corresponding free-energy densities as

f red (ρ 1 , ρ 2 , β) = f * (ρ 1 , ρ 2 , β) . (29)
Notice that this subtraction from the free-energy functional of the self-electrostatic energy was first introduced by Hohenberg and Kohn [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] for studying the quantum electron gas.

C. Density functional expansions for almost homogeneous systems

For states of S with slow spatial variations of the particle densities, the corresponding reduced free-energy F red {ρ 1 (•), ρ 2 (•)} can be expanded in powers of the gradients of ρ 1 (r) and ρ 2 (r). The leading term in that systematic expansion is purely local and reduces to

dr f * (ρ 1 (r), ρ 2 (r), β) , (30)
where we have used identity [START_REF] Meeron | Plasma Physics[END_REF] for the reduced free-energy density of an homogeneous system. The first correction, the so-called square-gradient term, reads [18, 19]

1 12 α,γ dr M * αγ (ρ 1 (r), ρ 2 (r), β)∇ρ α (r) • ∇ρ γ (r) , (31)
where M * αγ (ρ 1 , ρ 2 , β) is the second moment of the short-range part of the direct twobody correlations, namely c SR αγ (r) = c αγ (r) + βq α q γ /r, for S * with homogeneous densities (ρ 1 , ρ 2 ). Similarly to the emergence of the free-energy density of S * in the purely local term [START_REF] Alastuey | Virial expansions for quantum plasmas: Diagrammatic resummations[END_REF], the direct correlations of S * arise in the square-gradient term because all the respective distribution functions of S * and S with the same homogeneous densities are identical as established above. Notice that for systems with short-range interactions, the square-gradient expansion of the free energy F involves second moments of the direct correlations themselves. Here, thanks to the subtraction [START_REF] Salpeter | On Mayer's theory of cluster expansions[END_REF] of the electrostatic selfenergy, the square-gradient expansion of the reduced free energy F red involves second moments of the short-range part c SR αγ (r), which do converge thanks to the large-distance behaviour c αγ (r) ∼ -βq α q γ /r when r → ∞.

The second correction to the purely local contribution (30) involves fourth-order spatial derivatives of the densities, and it has been explicitly computed in Ref. [START_REF] Alastuey | Contribution of the three-particle and higher order correlations to the computation of interface density profiles by density functional theories, for two-dimensional plasmas[END_REF]. The corresponding local ingredients are fourth moments of two-, three-and four-body direct correlations of S * with homogeneous densities. Higher-order corrections would exhibit similar structures with well-behaved local ingredients defined for the same system.

IV. LINEAR RESPONSE TO A WEAK SLOWLY-VARYING EXTERNAL CHARGE DISTRIBUTION

We start with S in the absence of any applied external potential, namely that realization of S is nothing but the homogeneous neutral TCP. In a second step, let us introduce an external charge distribution c ext (r) = q ext exp(ik • r), with an infinitesimal amplitude q ext . Our aim here is to determine the induced charge density in S by DFT and compare to its linear response expression. This will provide the required fourth moment sum rule for the charge correlations of the homogeneous neutral TCP. In a first step, within DFT, we compute the density responses δρ α (r) to the external potentials, U α (r) = q α ϕ ext (r) with ϕ ext (r) the electrostatic potential created by c ext (r), at leading order in q ext and in the limit of small wave-numbers k → 0.

A. Analysis within density functional theory

Since the applied external potential varies on an infinitely large scale length, the particle densities should also display infinitely slow spatial variations. Then, the free-energy functional can be replaced by its density-gradient expansion introduced above, namely

F {ρ 1 (•), ρ 2 (•)} = β 2 drdr ′ (q 1 ρ 1 (r) + q 2 ρ 2 (r))(q 1 ρ 1 (r ′ ) + q 2 ρ 2 (r ′ )) |r ′ -r| + dr f * (ρ 1 (r), ρ 2 (r), β) + 1 12 α,γ dr M * αγ (ρ 1 (r), ρ 2 (r), β)∇ρ α (r) • ∇ρ γ (r) + ... , (32)
where the terms left over do not contribute to the deviations δρ α (r) at the considered lowest orders in k, as shown further on. The fundamental equation ( 21) of DFT then becomes for each species,

∂f * ∂ρ 1 - 1 6 [M * 11 ∆ρ 1 (r) + M * 12 ∆ρ 2 (r)] - 1 12
∂M * 11 ∂ρ 1 (∇ρ 1 (r)) 2 + 2 ∂M * 11 ∂ρ 2 ∇ρ 1 (r) • ∇ρ 2 (r) + (2 ∂M * 12 ∂ρ 2 - ∂M * 22 ∂ρ 1 )(∇ρ 2 (r)) 2 + ... = βµ 1 -βq 1 ϕ tot (r) (33)
and

∂f * ∂ρ 2 - 1 6 [M * 22 ∆ρ 2 (r) + M * 12 ∆ρ 1 (r)] - 1 12
∂M * 22 ∂ρ 2 (∇ρ 2 (r)) 2 + 2 ∂M * 22 ∂ρ 1 ∇ρ 1 (r) • ∇ρ 2 (r) + (2 ∂M * 12 ∂ρ 1 - ∂M * 11 ∂ρ 2 )(∇ρ 1 (r)) 2 + ... = βµ 2 -βq 2 ϕ tot (r) , (34)
where ϕ tot (r) is the total electrostatic potential created by the charge distribution (q 1 ρ 1 (r) + q 2 ρ 2 (r) + c ext (r)). In Eqts. [START_REF] Martin | A new proof of the Stillinger-Lovett complete shielding condition[END_REF][START_REF] Kalinay | The sixth-moment sum rule for the pair correlations of the two-dimensional one-component plasma: Exact results[END_REF] all involved quantities of S * are evaluated for homogeneous densities identical to the local densities (ρ 1 (r), ρ 2 (r)) and inverse temperature β. Moreover, all terms left over involve at least fourth-order spatial derivatives of ρ 1 (r) and ρ 2 (r).

In order to compute the induced densities at lowest order in q ext , we can linearize Eqts. [START_REF] Martin | A new proof of the Stillinger-Lovett complete shielding condition[END_REF][START_REF] Kalinay | The sixth-moment sum rule for the pair correlations of the two-dimensional one-component plasma: Exact results[END_REF] with respect to δρ α (r). The third terms in the l.h.s. of those Eqts. do not contribute anymore since they are at least of order q 2 ext . The resulting deviations take the form of plane waves, like the forcing external charge c ext (r), namely

δρ 1 (r) = A 1 (k) exp(ik • r) and δρ 2 (r) = A 2 (k) exp(ik • r) (35)
where the amplitudes A α (k) are proportional to q ext . The total electrostatic potential

ϕ tot (r) satisfies Poisson equation ∆ϕ tot (r) = -4π [q 1 ρ 1 (r) + q 2 ρ 2 (r) + c ext (r)] = -4π [q 1 δρ 1 (r) + q 2 δρ 2 (r) + c ext (r)] ( 36
)
where the second equality follows from the overall neutrality of the unperturbed system S . In order to eliminate ϕ tot (r) in favor of the induced density deviations, it is then sufficient to take the Laplacian of the linearized versions of Eqts. [START_REF] Martin | A new proof of the Stillinger-Lovett complete shielding condition[END_REF][START_REF] Kalinay | The sixth-moment sum rule for the pair correlations of the two-dimensional one-component plasma: Exact results[END_REF]. This provides [START_REF] Aqua | Ionic criticality: an exactly soluble model[END_REF] with χ -1 αγ = ∂ 2 f * /∂ρ α ∂ρ γ and a αγ = M * αγ /6. Those reference quantities are evaluated for the set (ρ 1 , ρ 2 ) ensuring overall neutrality of the unperturbed system S. Notice that if the thermodynamic function χ -1 αγ is specific to the enlarged auxiliary system S * , the microscopic second moments a αγ entirely depend on the direct correlations of the genuine system S of interest.

4πβq 2 1 + χ -1 11 k 2 + a 11 k 4 + O(k 6 ) A 1 (k) + 4πβq 1 q 2 + χ -1 12 k 2 + a 12 k 4 + O(k 6 ) A 2 (k) = -4πβq 1 δq ext (37) 4πβq 1 q 2 + χ -1 21 k 2 + a 21 k 4 + O(k 6 ) A 1 (k) + 4πβq 2 2 + χ -1 22 k 2 + a 22 k 4 + O(k 6 ) A 2 (k) = -4πβq 2 δq ext
The linear Eqts. [START_REF] Caillol | A Monte Carlo study of the dielectric constant of the restricted primitive model of electrolytes on the vapor branch of the coexistence line[END_REF][START_REF] Aqua | Ionic criticality: an exactly soluble model[END_REF] are straightforwardly solved in terms of the determinant of the associated two by two matrix which reads

D(k) = 4πβ(q 2 2 χ -1 11 + q 2 1 χ -1 22 -2q 1 q 2 χ -1 12 )k 2 + χ -1 11 χ -1 22 -χ -2 12 + 4πβ(q 2 2 a 11 + q 2 1 a 22 -2q 1 q 2 a 12 ) k 4 + O(k 6 ) . (39)
The amplitudes A α (k) are then found to be

A 1 (k) = 4πβ D(k) (q 2 χ -1 12 -q 1 χ -1 22 )k 2 + (q 2 a 12 -q 1 a 22 )k 4 + O(k 6 ) δq ext (40)
and

A 2 (k) = 4πβ D(k) (q 1 χ -1 12 -q 2 χ -1 11 )k 2 + (q 1 a 12 -q 2 a 11 )k 4 + O(k 6 ) δq ext . (41)
Therefore the proportionality coefficient between a given amplitude and δq ext , behaves in the limit k → 0 as a constant, which depends only on the thermodynamic quantities χ -1 αγ , plus a term of order k 2 which depends on both χ -1 αγ and a αγ . Now, if we form the induced charge density

δc(r) = q 1 δρ 1 (r) + q 2 δρ 2 (r) = c ind (k) exp(ik • r) (42)
with the charge amplitude

c ind (k) = q 1 A 1 (k) + q 2 A 2 (k) , (43)
we find

c ind (k) = -δq ext [1 - (χ -1 11 χ -1 22 -χ -2 12 ) 4πβ(q 2 2 χ -1 11 + q 2 1 χ -1 22 -2q 1 q 2 χ -1 12 ) k 2 + O(k 4 )] . (44)
Remarkably, the proportionality coefficient between the induced and external charges goes to -1 when k → 0, in relation with perfect screening properties, as discussed further.

Furthermore the term of order k 2 in its small-k expansion now depends only on the thermodynamical functions χ -1 αγ , and no longer on the microscopic quantities a αγ .

B. The fourth moment sum rule

The resulting induced charge density, can be also determined within linear response theory, which provides

δc(r) = - 4πβ k 2 S(k)δq ext exp(ik • r) (45)
In the linear response formula (45), S(k) is the Fourier transform of the charge correlations of the unperturbed system S, i.e. the homogeneous neutral TCP,

S(k) = dr exp(ik • r) α,γ q α q γ ρ αγ (r) + α q 2 α ρ α δ(r) (46)
with ρ αγ (r) the two-body probability density for the spatial configuration where one particle of species α is fixed at the origin, while another particle of species γ is fixed at r.

The small-k expansion of the amplitude c ind (k) can be inferred from the linear response formula (45) by inserting the corresponding expansion of S(k), for the zeroth and second moments of S(r), the DFT expression ( 44) also provides a new sum rule for the fourth moment, namely

S(k) = I 0 + I 2 k 2 + I 4 k 4 + ... , (47)
I 4 = - ρ 2 (4π(q 1 -q 2 )) 2 β 3 (χ -1 11 χ -1 22 -χ -2 12 )χ T , (49)
which follows from the identification of the k 2 -terms. The compressibility χ T emerges in that sum rule, thanks to the identity ( 16) rewritten in terms of the charges q 1 = Z 1 q and q 2 = -Z 2 q.

C. Related sum rules for other models

Let us first consider the case of the OCP. A fourth moment sum rule for the corresponding charge correlations S OCP (r) was derived by Vieillefosse and Hansen [START_REF] Vieillefosse | Statistical mechanics of dense ionized matter. Hydrodynamic limit and transport coefficients of the classical one-component plasma[END_REF] through a macroscopic analysis of fluctuations. In their textbook [START_REF] Hansen | Theory of simple liquids[END_REF], Hansen and Mac Donald propose a simple derivation which is similar to ours. They compute the charge density induced by a weak external plane wave charge distribution within an hydrostatic approach, where the force associated with the local pressure gradient is balanced by the total electrostatic force created by both the external and induced charges. Notice that the corresponding equation can be merely obtained by taking the gradient of the fundamental DFT equation [START_REF] Martin | A new proof of the Stillinger-Lovett complete shielding condition[END_REF] restricted to a single species and where all non-local contributions, including that involving the second moment of the direct correlations, are omitted. Moreover, the corresponding f * can then be obviously replaced by f OCP . The fourth moment of S OCP (r) then reduces to [START_REF] Vieillefosse | Statistical mechanics of dense ionized matter. Hydrodynamic limit and transport coefficients of the classical one-component plasma[END_REF][START_REF] Hansen | Theory of simple liquids[END_REF],

I OCP 4 = - 1 (4πqρ) 2 βχ OCP T . ( 50
)
Notice that this expression has been recovered through manipulations of the BGY hierarchy, for pure Coulomb interactions [START_REF] Vieillefosse | Sum rules and perfect screening conditions for the one-component plasma[END_REF] and also including short-range interactions [START_REF] Vieillefosse | Sum rules for the one-component plasma with additional shortrange forces[END_REF].

The OCP result has been extended to a multicomponent ionic mixture (MIM) of all positive point charges immersed in a rigid neutralizing background [START_REF] Suttorp | Equilibrium properties of a multi-component ionic mixture: I. Sum rules for correlation functions[END_REF][START_REF] Suttorp | Sum rules for correlation functions of ionic mixtures in arbitrary dimension d ≥ 2[END_REF]. Interestingly, the derivation is intrinsic and does not rely on the response of the system to a weak external charge distribution. Like the analysis [START_REF] Vieillefosse | Sum rules and perfect screening conditions for the one-component plasma[END_REF] carried out for the OCP, it is based on suitable manipulations of the BGY hierarchy equations for the distribution functions of the infinite homogeneous neutral system. A priori the derivation is only valid for pure Coulomb interactions, without any short range regularization which is unnecessary here since all mobile charges repel together. It makes an explicit use of the remarkable homogeneity property of the resulting pair interactions. The fourth moment of charge correlations in real space is then given by formula (7.3) in Ref. [START_REF] Suttorp | Sum rules for correlation functions of ionic mixtures in arbitrary dimension d ≥ 2[END_REF], which reduces in three dimensions (d = 3) to dr r 4 S MIM (r) = -120 β α,γ q α q γ ∂ρ γ /∂µ α (51)

where we have used that q v = α q α ρ α , while µ α = β -1 ∂f MIM /ρ α . Each partial derivative ∂ρ γ /∂µ α is computed by fixing the inverse temperature β as well as all µ δ 's with δ = α.

Straightforward manipulations of the multi-variable functions ρ γ (β, {µ α }) and µ α (β, {ρ γ })

allow us to express all partial derivatives ∂ρ γ /∂µ α in terms of partial derivatives ∂µ α /∂ρ γ .

In the binary case, we find

∂ρ 1 ∂µ 1 = ∂µ 2 ∂ρ 2 ∂µ 1 ∂ρ 1 ∂µ 2 ∂ρ 2 - ∂µ 1 ∂ρ 2 ∂µ 2 ∂ρ 1 -1 ; ∂ρ 2 ∂µ 2 = ∂µ 1 ∂ρ 1 ∂µ 1 ∂ρ 1 ∂µ 2 ∂ρ 2 - ∂µ 1 ∂ρ 2 ∂µ 2 ∂ρ 1 -1 ∂ρ 1 ∂µ 2 = - ∂µ 1 ∂ρ 2 ∂µ 1 ∂ρ 1 ∂µ 2 ∂ρ 2 - ∂µ 1 ∂ρ 2 ∂µ 2 ∂ρ 1 -1 ; ∂ρ 2 ∂µ 1 = - ∂µ 2 ∂ρ 1 ∂µ 1 ∂ρ 1 ∂µ 2 ∂ρ 2 - ∂µ 1 ∂ρ 2 ∂µ 2 ∂ρ 1 -1 . (52)
Using the identity

I 4 = 1 120 dr r 4 S(r) (53)
and inserting relations (52) into formula (51), we find that the corresponding I BIM 4 exactly coincides with our general expression (49) specified to the BIM, where the free-energy density f * merely reduces to f BIM . Indeed, our derivation also applies to the BIM where q 1 and q 2 now have the same sign, while the auxiliary system S * becomes identical to the genuine BIM of interest with the background charge density c B = -q v =α q α ρ α .

D. About other approaches

To our knowledge, in the literature, there exist two attempts to derive a sum rule for the fourth moment of the charge correlations of the TCP. First, the hydrodynamic approach carried out in Ref. [START_REF] Giaquinta | Hydrodynamic correlation functions for molten salts[END_REF] provides an expression for the fourth moment, different from formula (49), which involves ill-defined thermodynamic quantities as well as particle masses. Its validity is then quite doubtful, in particular because classical equilibrium charge correlations do not depend on particle masses.

Second, van Beijeren and Felderhof [START_REF] Van Beijeren | Thermodynamics and correlation functions of plasmas and electrolyte solutions[END_REF] proceed to an intrinsic analysis of charge correlations within the grand-canonical ensemble, where they combine the Ornstein-Zernicke equations with DFT manipulations. In agreement with results previously derived by Mitchell et al. [START_REF] Mitchell | On the second-moment condition of Stillinger and Lovett[END_REF], who shown that the fourth moment cannot be expressed in terms of thermodynamic quantities of the sole TCP, they find that it is necessary to introduce non-neutral states of the TCP which can be realized through the application of a suitable external potential. However, they did not provide any scheme which determines that external potential. Thus their free-energy density f 0 , from which the thermodynamical chemical potentials are inferred through the usual identity written in formula (3.12) of

Ref. [START_REF] Van Beijeren | Thermodynamics and correlation functions of plasmas and electrolyte solutions[END_REF], remains a formal quantity, with no prescriptions for explicit calculations. This ambiguity might explain why their work is not always cited. According to our analysis, it can be easily clarified as follows. In fact, as shown in Section II, the external potential mentioned in Ref. [START_REF] Van Beijeren | Thermodynamics and correlation functions of plasmas and electrolyte solutions[END_REF] is nothing but our potential U B α (r) = q α ϕ B (r) where ϕ B (r) is the electrostatic potential created by an homogeneous background density. Therefore, f 0 is identical to our free-energy density f * of the TCP immersed in an uniform rigid background. Then, the relation between partial derivatives (52) allows us to exactly recast formula (6.26) of Ref. [START_REF] Van Beijeren | Thermodynamics and correlation functions of plasmas and electrolyte solutions[END_REF] as our expression (49), similarly to what occurs for the corresponding formula obtained for the BIM by Suttorp [START_REF] Suttorp | Sum rules for correlation functions of ionic mixtures in arbitrary dimension d ≥ 2[END_REF].

V. ASYMPTOTIC EXPANSIONS AT LOW DENSITIES

It is instructive to check the fourth moment sum rule for specific models and various ranges of thermodynamical parameters. Here, we consider the model of charged soft spheres with the pair interaction u αγ (r) given by formula [START_REF] Martin | Sum rules in charged fluids[END_REF]. First, we briefly describe how the pair correlations of S * can be represented by an infinite series of resummed Mayer graphs. Such resummed diagrammatics constitute a quite suitable framework for deriving low-density expansions of the quantities of interest. From the diagrammatic representation of charge correlations, we infer the low-density expansion of I 4 defined as the coefficient of the k 4 -term in the small-k expansion (47) of S(k). The diagrammatics for the pair correlations also give access to the free energy density f * through thermodynamical identities. The low-density expansion of the thermodynamical expression (49) is then computed, and it is shown to exactly match that of I 4 , as expected.

A. Exploiting the principle of topological reduction

Let ρ * αγ,T (r) = ρ * αγ (r) -ρ α ρ γ be the truncated pair distribution functions of S * , also called pair correlations, for an arbitrary set of densities (ρ 1 , ρ 2 ). As argued above, the distributions function of S * , which includes a background with charge density c B = -(q 1 ρ 1 + q 2 ρ 2 ), are identical to that of of a purely two-component system where the mobile particles are submitted to the external potential ( 22) U B α (r) created by the background. Therefore, pair correlations ρ * αγ,T (r) are represented by series of Mayer diagrams [START_REF] Mayer | Statistical Mechanics[END_REF] made with two root (white) points respectively fixed at the origin 0 and at r, and an arbitrary number of black points whose positions are integrated over. Each point carries a statistical weight

z α = exp[β(µ α -U B α )] (2πλ 2 α ) 3/2 , ( 54
)
while two points are connected by at most one Mayer bond

b M = exp(-βu αγ ) -1 . ( 55
)
Each diagram is simply connected, namely there exists at least one path connecting two arbitrary points.

The previous Mayer diagrams are difficult to handle because the fugacity weights (54) are inhomogeneous and depend on the positions of the points. A great simplification can be achieved by virtue of the principle of topological reduction, nicely exposed in Ref. [START_REF] Hansen | Theory of simple liquids[END_REF],

which consists in removing all articulation points. An articulation point is such that there exists at least one subdiagram attached to it and not connected to the rest of the diagram. In other words, the suppression of the articulation point leaves that subdiagram disconnected from the two root points. If one sums all those subdiagrams attached to a given articulation point, all articulation points are removed, while simultaneously all fugacity weights (54) are replaced by density weights ρ α [START_REF] Hansen | Theory of simple liquids[END_REF]. Furthermore, the topological structure of the diagrams is conserved through that reduction. Accordingly, the pair correlations ρ * αγ,T (r) are represented by Mayer diagrams made with the two root points fixed at 0 and r, and an arbitrary number of black points, where the point statistical weights are now the densities ρ α . Two point are still connected at most by one Mayer bond (55). Each diagram is again simply connected but is now free of any articulation point.

Thanks to the translational invariance of both density weights and Mayer bonds, the Mayer density diagrams reveal quite useful for explicit calculations as described further.

Notice that, remarkably, the background does not show in such diagrams, its effects being implicitly and entirely taken into account by the introduction of the homogeneous densities ρ α

B. Abe-Meeron resummations

Because of the long-range non-integrable decay of two-body interactions u αγ , every Mayer diagram diverges. All those divergencies can be removed via chain resummations, as first noticed by Mayer [START_REF] Mayer | The theory of ionic solutions[END_REF] and Salpeter [START_REF] Salpeter | On Mayer's theory of cluster expansions[END_REF], and then performed in a systematic way for the whole diagrammatical series by Abe [START_REF] Abe | Giant cluster expansion theory and its application to high temperature plasma[END_REF] and Meeron [START_REF] Meeron | Plasma Physics[END_REF]. A simplified presentation of that method can be found in Refs. [START_REF] Alastuey | Virial expansions for quantum plasmas: Diagrammatic resummations[END_REF] and [START_REF] Alastuey | Screened cluster expansions for partially ionized gases[END_REF]. It starts with the decomposition of each Mayer bond (55) as

b M = b T M -βq α q γ v C ( 56
)
with the truncated bond

b T M = exp(-βu αγ ) -1 + βq α q γ v C ( 57
)
and the Coulomb potential v C (r) = 1/r. After inserting the decomposition (56) into every Mayer diagram, one proceeds to systematic resummations of convolution chains of Coulomb bonds -βq α q γ v C . Thanks to remarkable combinatorial properties [START_REF] Alastuey | Screened cluster expansions for partially ionized gases[END_REF], all those resummations can be performed in terms of a single effective potential, which is nothing but the well-known Debye potential

φ D (r) = exp(-κ D r) r (58)
with the Debye inverse length κ D = ( α 4πβq 2 α ρ α ) 1/2 . The chain resummations give raise to two bonds [START_REF] Alastuey | Virial expansions for quantum plasmas: Diagrammatic resummations[END_REF], the Debye bond

b D = -βq α q γ φ D ( 59
)
and the short-range dressed bond

b R = exp(-β(u SR αγ + q α q γ φ D )) -1 + βq α q γ φ D , (60)
with the short-range part of pair interactions

u SR αγ = u αγ -q α q γ v C .
ρ * αγ,T (r) = ρ α ρ γ G 1 S G n i=1 α i dr i ρ α i b D b R G . ( 61
)
The so-called prototype graphs G are made with the two root points respectively fixed at includes the short-range regularization which also guarantees its integrability. We stress that representation (61) holds for any set (ρ 1 , ρ 2 ) of densities, and then appears to be quite useful for computing equilibrium quantities of S * . Moreover it is valid for any short-range regularization u SR αγ , including of course that describing soft or hard spheres.

C. Charge correlations

The Fourier transform (46) of the charge correlations of the homogeneous TCP can be recast as

S(k) = α,γ q α q γ ραγ (k) + α q 2 α ρ α (62)
where the Fourier transform ραγ (k) of pair correlations is given by the sum of the Fourier transforms of the contributions of all graphs G in the representation (61). Let us first consider the contribution of the simplest graph G D shown in Fig. 1, where the two root points are connected by a Debye bond. Its contribution to S(k) added to the constant term α q 2 α ρ α in the formula (62) provides the well-known Debye charge correlations

SD (k) = κ 2 D 4πβ k 2 k 2 + κ 2 D , (63)
which can be derived within a mean-field treatment of correlations, without any diagrammatic considerations. Now we stress that SD (k) saturates the first two moments sum rules for I 0 and I 2 , since SD (k) ∼ k 2 /(4πβ) when k → 0. Therefore all the remaining graphs in the representation (61) give no contributions to I 0 and I 2 . That remarkable property is related to the following reorganization of the series of graphs, which turns out to be also

K αγ,G d (r) = n i=1 α i dr i ρ α i b D b R G d (64)
and noticing that all four graphs G d , G Dd , G dD and G DdD have the same symmetry factor S G d , we can rewrite that total contribution as ρ α ρ γ /S G d times

K αγ,G d (r) - α 1 ρ α 1 dr 1 [βq α q α 1 φ D (r 1 )K α 1 γ,G d (|r -r 1 |) + K αα 1 ,G d (r 1 )βq α 1 q γ φ D (|r -r 1 |)] + α 1 ,α 2 ρ α 1 ρ α 2 dr 1 dr 2 βq α q α 1 φ D (r 1 )K α 1 α 2 ,G d (|r 2 -r 1 |)βq α 2 q γ φ D (|r -r 2 |) (65)
The corresponding contribution to S(k) in formula (62) can be readily computed by using the convolution theorem and φD (k) = 4π/(k 2 + κ 2 D ), with the result

k 4 (k 2 + κ 2 D ) 2 α,γ ρ α ρ γ q α q γ Kαγ,G d (k) (66)
divided by the symmetry factor S G d . Since that expression is at least of order k 4 when k → 0 for any G d , all the graphs beyond the Debye graph G D do not contribute neither to I 0 , nor to I 2 . Moreover, because of the prefactor of order k 4 , the resulting contribution to I 4 makes Kαγ,G d (0) appear. After adding the simple contribution of G D computed from the Debye formula (63), we eventually obtain the diagrammatic representation of I 4 ,

I 4 = - 1 4πβκ 2 D + 1 κ 4 D G d 1 S G d α,γ ρ α ρ γ q α q γ Kαγ,G d (0) . (67)
Representation ( 67) is well-suited for computing the low-density expansion of I 4 . Indeed, and as usual, because of the density weights carried by the black points, only a finite number of graphs contribute up to a given order. However, here we have to take care of the dependence on the density of the bonds b D and b R through the Debye wavenumber

κ D = ( α 4πβq 2 α ρ α ) 1/2 .
Consequently, the order of a contribution is not merely given by counting the number of black points on the one hand, while half-integer powers and logarithmic terms arise in the expansion on the other hand. We have computed the first three terms of that expansion, up to constant terms of order ρ 0 included. In the Appendix, we provide some technical details, as well as the complete list of graphs G d which contribute up to the considered order. The resulting expansion reads

I 4 = - 1 4πβκ 2 D + πβ 2 κ 5 D α,γ q 3 α q 3 γ ρ α ρ γ + 2πβ 3 3κ 4 D α,γ q 4 α q 4 γ ρ α ρ γ ln(8κ D d αγ ) + 4πβ κ 4 D α,γ q 2 α q 2 γ ρ α ρ γ d 2 αγ - 3πβ 2 κ 4 D α,γ q 3 α q 3 γ ρ α ρ γ d αγ + 1 κ 4 D α,γ q α q γ ρ α ρ γ dr exp(-βu αγ ) -1 + βu αγ -β 2 u 2 αγ /2 + β 3 u 3 αγ /6 + 11π 2 β 4 3κ 6 D α,γ,δ q 3 α q 3 γ q 4 δ ρ α ρ γ ρ δ - 52π 3 β 5 9κ 8 D α,γ,δ,η q 3 α q 3 γ q 3 δ q 3 η ρ α ρ γ ρ δ ρ η + o(ρ 0 ) . ( 68
)
Not surprisingly, the leading term in the expansion (68) is the purely Debye contribution, and it behaves as 1/ρ. The next correction is of order 1/ρ 1/2 and is independent of the short-range part of the interactions. Those short-range parts arise in further corrections of order ln ρ and ρ 0 . The last two terms of order ρ 0 are purely Coulomb contributions.

The expansion is valid for any set of densities (ρ 1 , ρ 2 ), including of course the neutral sets defining the neutral TCP. For other short-range interactions, like hard cores for instance, the structure of the low-density expansion of I 4 is identical to that (68) explicitly computed for charged soft spheres.

D. Free energy density

The excess free energy of S * for any set of homogeneous densities, can be obtained through the usual integration over the inverse temperature of the equilibrium average of the potential part of Hamiltonian [START_REF] Das | When is a conductor not perfect? Sum rules fail under critical fluctuations[END_REF]. The resulting free energy density in thermal units

f * (ρ 1 , ρ 2 , β) reduces to f * (ρ 1 , ρ 2 , β) = f id (ρ 1 , ρ 2 , β) + 1 2 α,γ β 0 dτ drρ * αγ,T (r)u αγ (r) + β 2 α,γ ρ α ρ γ dru SR αγ (r) . (69)
In that formula, the first term takes the familiar form

f id (ρ 1 , ρ 2 , β) = ρ 1 ln(ρ 1 (2πλ 2 1 ) 3/2 ) -1 + ρ 2 ln(ρ 2 (2πλ 2 2 ) 3/2 ) -1 (70)
which describes a mixture of ideal gases. The next two terms account for interactions. The second term involving pair correlations ρ * αγ,T (r) is obtained by adding and substracting u SR αγ to the purely Coulomb interactions in the particle-backround and background-background parts of the potential energy of Hamiltonian [START_REF] Das | When is a conductor not perfect? Sum rules fail under critical fluctuations[END_REF]. This also provides the third term which merely reduces to

β 2 α,γ ρ α ρ γ dru SR αγ (r) = -2πβ α,γ q α q γ ρ α ρ γ d 2 αγ . (71)
It is implicitly understood that pair correlations ρ * αγ,T (r) in formula (69) are evaluated at inverse temperature τ .

The contribution of pair correlations ρ * αγ,T (r) to f * follows by inserting its resummed diagrammatic representation (61) into the second term of formula (69). The contribution

of a given diagram G reads 1 2S G α,γ ρ α ρ γ β 0 dτ dr u αγ (r)K (τ ) αγ,G (r) (72)
where

K (τ )
αγ,G (r) is the integral (64) over black points with G in place of G d and bonds evaluated at inverse temperature τ . In the low-density limit, up to a given order in the density, only a finite number of contributions (72) needs to be retained. In the Appendix, we provide a few technical details of the calculations, as well as the list of graphs G which contribute to f * up to order ρ 2 included. Adding the simple ideal (70) and background (71) contributions, we eventually obtain the low-density expansion of f * ,

f * = ρ 1 ln(ρ 1 (2πλ 2 1 ) 3/2 ) -1 + ρ 2 ln(ρ 2 (2πλ 2 2 ) 3/2 ) -1 - κ 3 D 12π -2πβ α,γ q α q γ ρ α ρ γ d 2 αγ + 3πβ 2 2 α,γ q 2 α q 2 γ ρ α ρ γ d αγ - πβ 3 3 α,γ q 3 α q 3 γ ρ α ρ γ ln(8κ D d αγ ) - 1 2 α,γ ρ α ρ γ dr exp(-βu αγ ) -1 + βu αγ -β 2 u 2 αγ /2 + β 3 u 3 αγ /6 + o(ρ 2 ) . ( 73
)
The leading terms of order ρ ln ρ in the expansion (73) are ideal contributions. The next correction of order ρ 3/2 arises from pure Coulomb interactions, and is nothing but the well-known Debye term. Contributions from the short range part of the interactions appear in the terms of order ρ 2 ln ρ and ρ 2 . The terms left over are least of order ρ 5/2 ln ρ.

Expansion (73) is valid for any set of densities (ρ 1 , ρ 2 ), and it gives access to all the other thermodynamical functions of S * through suitable partial derivatives with respect to the independent parameters β, ρ 1 or ρ 2 defining an homogeneous equilibrium state of S * .

For other short-range regularizations, the low-density expansion of f * has the same structure as (73). However, notice that for hard core potentials, the ideal term in the decomposition (69) of the corresponding f * must be replaced by the free energy density of hard spheres, f HS (ρ 1 , ρ 2 , β). In the low-density limit, f HS (ρ 1 , ρ 2 , β) can be expanded in entire powers of ρ around the ideal term (70). In the resulting full expansion of f * , there are terms which depend only on the hard core diameters σ αγ and not on the particles charges.

E. Checking the sum rule at lowest orders

In order to check the sum rule (49), we first have to compute the low-density expansion of the partial compressibilities χ -1 αγ = ∂ 2 f * /∂ρ α ∂ρ γ . Using expansion (73) of f * , we find

χ -1 11 = 1 ρ 1 - πβ 2 q 4 1 κ D + 3πβ 2 q 4 1 d 11 -4πβq 2 1 d 2 11 - 2πβ 3 q 6 1 3 ln(8κ D d 11 ) - 8π 2 β 4 q 5 1 3κ 2 D α q 3 α ρ α + 8π 3 β 5 q 4 1 3κ 4 D αγ q 3 α q 3 γ ρ α ρ γ -dr exp(-βu 11 ) -1 + βu 11 -β 2 u 2 11 /2 + β 3 u 3 11 /6 + o(ρ 0 ) , (74)
χ -1 22 = 1 ρ 2 - πβ 2 q 4 2 κ D + 3πβ 2 q 4 2 d 22 -4πβq 2 2 d 2 22 - 2πβ 3 q 6 2 3 ln(8κ D d 22 ) - 8π 2 β 4 q 5 2 3κ 2 D α q 3 α ρ α + 8π 3 β 5 q 4 2 3κ 4 D αγ q 3 α q 3 γ ρ α ρ γ -dr exp(-βu 22 ) -1 + βu 22 -β 2 u 2 22 /2 + β 3 u 3 22 /6 + o(ρ 0 ) , (75)
and

χ -1 12 = χ -1 21 = - πβ 2 q 2 1 q 2 2 κ D + 3πβ 2 q 2 1 q 2 2 d 12 -4πβq 1 q 2 d 2 12 - 2πβ 3 q 3 1 q 3 2 3 ln(8κ D d 12 ) - πβ 3 q 3 1 q 3 2 3 - 4π 2 β 4 q 2 1 q 2 2 3κ 2 D α q 4 α ρ α + 8π 3 β 5 q 2 1 q 2 2 3κ 4 D αγ q 3 α q 3 γ ρ α ρ γ -dr exp(-βu 12 ) -1 + βu 12 -β 2 u 2 12 /2 + β 3 u 3 12 /6 + o(ρ 0 ) . ( 76
)
Notice that the leading contributions of order 1/ρ in both χ -1 11 and χ -1 22 arise from the ideal terms in f * , while the next correction of order 1/ρ 1/2 comes from the Debye term in expansion (73). The leading contribution of order 1/ρ 1/2 in χ -1

12 is also provided by that Debye correction. All terms which are left over in expansions (74), (75) and (76) are at least of order ρ 1/2 ln ρ.

According to the expression [START_REF] Percus | Equilibrium state of a classical fluid of hard rods in an external field[END_REF] of the isothermal compressibility, the thermodynam-quantities arise on the one hand, while perfect screening properties can be lost on the other hand, as mentioned below.

Our derivation also involves implicit assumptions about the existence of the thermodynamic limit, and of intrinsic bulk properties with bulk densities which become homogeneous far from the boundaries. Strictly speaking, to our knowledge, this has been only proved for the general three-dimensional TCP in the Debye regime [START_REF] Brydges | Debye screening[END_REF] and for its charge-symmetric version [START_REF] Frohlich | Correlation inequalities and the thermodynamic limit for classical and quantum continuous systems[END_REF]. In two dimensions, where the Coulomb potential takes a logarithmic form, both the neutrality and homogeneity of a TCP of point charges have been proved [START_REF] Leblé | Large deviations for the two-dimensional two-component plasma[END_REF]. Extensions of such results to all the systems introduced here would be quite valuable of course, and might constitute the first steps towards a complete proof of our sum rule. Meanwhile, physical arguments, in particular related to the beautiful proof for quantum Coulomb matter [START_REF] Lieb | The constitution of matter: existence of thermodynamics for systems composed of electrons and nuclei[END_REF], strongly suggest that the classical TCP, as well as its version immersed in a charged uniform background, do sustain a well-behaved TL.

Furthermore, there are strong evidences, arising either from specific models or mean-field approaches, that screening properties in the bulk can be disentangled from the reorganization of charges at the surface, so any boundary effects can be indeed a priori ignored.

In the absence of a complete mathematical proof, checking the sum rule within exact calculations for specific models or thermodynamical regimes is particularly valuable.

Here, such checking has been carried out for charged soft spheres in the low density regime, through the explicit calculation of the lowest order terms in density expansions of the quantities of interest. This illustrates the subtle interplay between short-range and screened Coulomb contributions which ultimately ensure the validity of the fourth moment sum rule at the considered orders. If there exists a simple reorganization of the Abe-Meeron diagrams which shows the validity of the second moment Stillinger Lovett sum rule at any order in the density expansion, a similar trick for the fourth moment sum rule, certainly more cumbersome, remains to be discovered.

If our derivation of the fourth moment sum rule is based on the response to external perturbations, more intrinsic derivations would be of great interest, both for enforcing its expected validity on the one hand, and for sheding light on the internal mechanisms at work on the other hand. For instance, the second moment Stillinger Lovett sum rule can be retrieved within suitable manipulations of the BGY hierarchy equations for the equilibrium distribution functions of the unperturbed homogeneous TCP, as shown by Gruber and Martin [START_REF] Martin | A new proof of the Stillinger-Lovett complete shielding condition[END_REF]. Moreover, the BGY hierarchy equations have been also used for deriving the fourth moment sum rule for the OCP [START_REF] Vieillefosse | Sum rules for the one-component plasma with additional shortrange forces[END_REF] and for the BIM with pure Coulomb interactions [START_REF] Suttorp | Sum rules for correlation functions of ionic mixtures in arbitrary dimension d ≥ 2[END_REF]. We are looking for extending such derivations to the TCP case, where the presence of short-range interactions requires further manipulations. Notice, that a full reorganization of Abe-Meeron diagrammatics as described above could be also seen as an intrinsic derivation. In the same spirit, let us mention that a sixth moment sum rule for the charge correlations of the two-dimensional OCP with logarithmic interactions was established through a full term by term analysis of the Abe-Meeron diagrammatics [START_REF] Kalinay | The sixth-moment sum rule for the pair correlations of the two-dimensional one-component plasma: Exact results[END_REF].

The fourth moment sum rule obviously extends to the TCP immersed in a charged background, namely the thermodynamical expression (49) of the fourth moment of charge correlations is not restricted to densities satisfying overall neutrality, but it is valid for any set of densities (ρ 1 , ρ 2 ). This is well illustrated by the low density calculations for charged soft spheres. Moreover, a similar DFT analysis combined with linear response theory can be carried out for an arbitrary number n of components. This would lead to formulae analogous to the thermodynamical expression (49), but with a more complicated structure arising from the inversion of a n×n matrix. In the case of the MIM, they should be equivalent to those derived by Suttorp [START_REF] Suttorp | Sum rules for correlation functions of ionic mixtures in arbitrary dimension d ≥ 2[END_REF]. Eventually, within our formalism, one can obtain sum rules for the zeroth and second moment of particle correlations, by comparing the DFT calculation of a given particle density to its linear response expression. Such sum rules are again equivalent to those obtained for the MIM [START_REF] Suttorp | Sum rules for correlation functions of ionic mixtures in arbitrary dimension d ≥ 2[END_REF]. Among the various possible applications of our new sum rule, we would like to emphasize its usefulness for a better understanding of the plausible lack of screening properties at the ionic critical point. The liquid-gas transition of a TCP has been widely studied the last twenty years. Let us mention for instance two recent works [START_REF] Das | Near critical electrolytes: Are the charge-charge sum rules obeyed?[END_REF], [START_REF] Caillol | Liquid-vapor transition and critical behaviour of the ultrasoft restricted primitive model of polyelectrolytes: A Monte Carlo study[END_REF]. Numerical simulations have convincingly shown that both liquid and gas phases display perfect screening properties, namely the second moment Stillinger-Lovett sum rule is satisfied. However, a first suspicion about the violation of that sum rule at the critical point was pointed out by Caillol [START_REF] Caillol | A Monte Carlo study of the dielectric constant of the restricted primitive model of electrolytes on the vapor branch of the coexistence line[END_REF]. Meanwhile, such violation was also observed for a solvable asymmetric mean-spherical model by Aqua and Fisher [START_REF] Aqua | Ionic criticality: an exactly soluble model[END_REF], which is expected to share common properties with an asymmetric TCP. More recently, and contrarily to various theoretical expectations, the violation of the Stillinger Lovett sum rule at the critical point was also observed for the fully symmetric RPM by Das, Kim and Fisher [START_REF] Das | When is a conductor not perfect? Sum rules fail under critical fluctuations[END_REF] : they provide strong numerical evidences by combining refined Monte Carlo simulations in the grand-canonical ensemble with finite-size scaling methods. Furthermore, they also show that the fourth moment of charge correlations diverges when approaching the critical point, in a way analogous to the isothermal compressibility. Clearly, our thermodynamical expression (49) of that fourth moment constitutes a promising tool for analyzing its behaviour near the critical point, as well as the underlying coupling between charge and mass fluctuations.

Eventually, let us conclude by a few comments regarding the two dimensional case. In two dimensions (2D), the Coulomb potential, defined as the solution of Poisson equation, takes the well-known logarithmic form. Thanks to the relatively weak singularity at the origin of the 2D Coulomb potential, the TCP of point charges is well behaved for coupling constants Γ < 2, namely at sufficiently high temperatures [START_REF] Deutsch | The partition function of a two-dimensional plasma[END_REF]. However, shortrange interactions need to be introduced for Γ ≥ 2 in order to avoid the collapse between positive and negative charges. The fourth moment sum rule derived here explicitly in three dimensions (3D), can be straightforwardly extended to the 2D case : this leads to the simple replacement of the factor 4π in formula (49) by the factor 2π, a direct consequence of the modification of Poisson equation when changing from 3D to 2D. At Γ = 2, Cornu and Jancovici [START_REF] Cornu | Two-dimensional Coulomb systems: A larger class of solvable models[END_REF] exploited a mapping with a field theory model valid for pure Coulomb interactions, which allowed them to derive analytical expressions for particle correlations of the 2D TCP. Such expressions are expected to become exact in the zero density limit for the well-behaved TCP including short-range interactions. Thus, they constitute a reliable starting point for further checking of the fourth moment sum rule at finite coupling. Similarly to its application to the study of ionic criticality in 3D, the fourth moment sum rule should also bring new insights for the celebrated Kosterlitz-Thouless transition [START_REF] Kosterlitz | Ordering, metastability and phase transitions in twodimensional systems[END_REF], [START_REF] Minnhagen | The two-dimensional Coulomb gas, vortex unbinding, and superfluidsuperconductive thin films[END_REF] : in the temperature-density plane, there appears a line of critical points separating a high-temperature conducting phase from a low-temperature dielectric phase [START_REF] Alastuey | Critical line near the zero-density critical point of the Kosterlitz-Thouless transition[END_REF]. The implications of the fourth moment sum rule should complete the results of a previous work [START_REF] Alastuey | Correlations in the Kosterlitz-Thouless phase of the twodimensional Coulomb gas[END_REF], where a plausible scenario for the large-distance decay of particle correlations in the dielectric phase was constructed in a way consistent with various sum rules. can be computed at lowest order in the density by merely replacing u SR αγ + q α q γ φ D by the bare pair potential u αγ since (exp(-βu αγ ) -1 + βu αγ -β 2 u 2 αγ /2 + β 3 u 3 αγ /6) is integrable in the whole space. The next density-dependent corrections to that leading contribution behave as ρ 1/2 and can thus be neglected in the considered calculation of I 4 up to order ρ 0 . The contributions of the other terms in the decomposition (81) are easily computed thanks to the simple analytic expressions of u SR αγ and φ D . Eventually, combining the above convolution and decomposition tricks, we obtain formula (68) for I 4 .

The five graphs in the series (61) for particle correlations which provide contributions (72) to the free-energy density f * are listed in Fig. 4. Each contribution follows from formula (72), so the value K given by expression (79), and of φD (k)4π/k 2 which reduces to 1 (2π) 3 dk exp(-ik • r)

16π 2 k 2 (k 2 + κ 2 D ) = 4π κ 2 D r (1 -exp(-κ D r)) . (83)
The final integrals over τ from 0 to β are then easily and explicitly performed for all terms which reduce to combinations of powers laws and logarithmic terms. It remains a term involving dr u αγ [exp(-τ u αγ ) -1 + τ u αγ -

τ 2 2 u 2 αγ ] , (84)

which only involves powers of k 2

 2 thanks to the expected exponential decay of charge correlations in real space. If we compare the resulting expansion of c ind (k) with the DFT result (44), we readily find I 0 = 0 and I 2 = 1 4πβ (48) which follow from respectively the absence of a 1/k 2 -term, and the identification of the constant terms. The vanishing of I 0 accounts for the perfect screening of internal charges. The universal value of I 2 , first demonstrated a long ago by Stillinger and Lovett [4], ensures the perfect screening of weak external charges. Beyond those well-known results

  The topological structure of the genuine Mayer diagrams remain unchanged, with bonds which can be either b D or b R , and with the additional rule excluding convolutions b D * b D in order to avoid double counting. Within the Abe-Meeron resummations, the genuine whole set of Mayer diagrams representing ρ * αγ,T (r) is then exactly transformed into

0FIG. 1 :

 1 FIG. 1: The Debye diagram in the resummed diagrammatic representation (61) of particle correlations. The two root (white points) are fixed at 0 and r respectively. The straight line represents a Debye bond b D (59).

quite useful for computing the fourth moment I 4 .FIG. 2 :

 42 FIG. 2: The three dressed Debye diagrams associated with a given diagram G d in the resummed diagrammatic representation (61) of particle correlations.

ACKOWLEDGMENTFIG. 3 : 3 (k 2 + κ 2 D ) 3 = 2π 2 κ 3 D( 1 +FIG. 4 :

 3323314 FIG. 3:The seven diagrams G d which contribute to I 4 up to order ρ 0 included in the diagrammatic series (67). The bubbles represent short-range dressed bonds b R (60).

  (r) of each graph is first computed with bonds b D and b R evaluated at temperature τ . After multiplication of K (τ ) αγ,G (r) by the pair potential u αγ (r), the further integrals over r in the whole space are readily computed by using decomposition (81) with τ in place of β, as well as the inverse Fourier transforms of [ φD (k)] 2

ical quantity in the right hand side of sum rule (49) can be rewritten as

The low-density expansion of that thermodynamical expression is straightforwardly computed by using the expansions (74), ( 75) and (76) of the χ -1 αγ 's. Its leading behaviour is immediately obtained by noticing that both χ -1 11 and χ -1 22 diverge faster than χ -1 12 in the zero-density limit, and the corresponding purely ideal behaviours χ -1 11 ∼ 1/ρ 1 and

which coincides with the leading term in expansion (68) of I 4 . The calculation of the next correction of order 1/ρ 1/2 remains simple, since it requires to retain only the first Debye corrections of order 1/ρ 1/2 to the ideal terms in both χ -1 11 and χ -1 22 , while χ -1 12 can be replaced by its leading Debye behaviour. The determination of the terms of order ρ 0 and ρ 0 ln ρ is still straightforward but more cumbersome. Eventually, we find that all those corrections to the ideal behaviour (78) of the thermodynamical quantity (77) exactly match the low-density expansion (68) of I 4 inferred from its microscopic definition. Thus, the fourth moment sum rule perfectly works, at least up to the considered order in the density.

VI. CONCLUDING COMMENTS AND PERSPECTIVES

In this paper, we have derived a new sum rule for the fourth moment of charge correlations of a TCP. Since the Stillinger-Lovett second moment sum rule naturally emerges as a by-product of our analysis, we believe that this new sum rule holds in any conducting phase, although all the steps of its derivation are not under a complete mathematical control at the moment. In particular, we expect that the free energy functional can be safely expanded around homogeneous states inside the conducting phase. We stress that all the partial compressibilities of the auxiliary system, namely second partial derivatives with respect to particle densities of the free energy density, have then to be well defined. Thus critical points must be dealt with some care, since singularities in the thermodynamical