Samuela ; Pasquali

Jerome K Percus

Mean Field and the Confined Single Homopolymer

Keywords: confined polymers, folding, long-range interactions PACS numbers: 05.20.Gg, 82.35.Lr, 87.15.A-I

come

in thermal equilibrium is a highly developed discipline. Along with specific physically motivated approximations have come tools of more general validity and utility. Our focus will be on a single polymer chain confined by external forces, and to minimize the needed information input, on single homopolymers. We will also attend in the main to idealized models in which the polymer is simply a chain of unit monomers of a few degrees of freedom, but will indicate how this restriction can be rewardingly removed, on the way to a realistic polymer representation.

We will aim at both analytic simplicity and reasonable suitability for the ultimately necessary computational procedures. For the former, we will work in a "monomer grand-ensemble" [7,8] in which the number of monomers per polymer is distributed, but will show that this need not be a drawback. For the latter, we will favor minimum principles to be able to better control computations. This does restrict the category of systems to be studied e.g. to purely attractive pair interactions, and will be generalized in later work, soon to be reported, as such restrictions are removed.

II. THE REFERENCE SYSTEM

The statistical mechanics of a free classical chain molecule, interpreted as a random walk is an old, old topic, although less attention has been paid to the effect of a structured external field, in part because the formal simplification afforded by translation invariance is then not available. In this case, the next neighbor pair energy then depends separately on the locations of each of the monomers in the pair, rather than on just their separation (in the point monomer model). We will tackle this problem by first (section II A) concealing the discreteness of the monomer population number, and then (section II B) concealing the inherent required inversion of the pair Boltzmann factor kernel.

A. Notation

Let us be a bit more explicit. We have in mind an ordered chain of N equivalent monomers, the j th being specified by its degrees of freedom r j . The order is maintained by a symmetric next neighbor interaction potential of Boltzmann factor w(r i , r i+1 ), depending of course on the inverse temperature β. Any two monomers can also interact via an interaction Boltzmann factor e(r i , r j ) = exp(-βφ(r i , r j )) and, crucially for the applications we have in mind, the polymer is constrained by an external potential u(r i ).

It is now convenient to imagine that the homopolymer in question is both in thermal equilibrium and in number equilibrium, i.e. that it is the result of monomer addition and absorption from a bath of non-interacting monomers (figure 1). The reaction equilibrium is analogous to that of the grand canonical ensemble for a fluid, but with a significant difference. Suppose that Q (m) is the monomer canonical partition function in its center of mass coordinate system, Q

N of the Nmonomer polymer (N ≥ 1, to recognize an object as a polymer). Then if the full system contains N ′ monomers, N of which are bound together, in a volume V , the full system partition function will be:

Q tot N ′ = N ≥1 (V Q (m) ) N ′ -N (N ′ -N )! Q (p) N . (1)
Since the monomers in the absence of a polymer would have a partition function:

Q (m) N ′ = (V Q (m) ) N ′ N ′ ! (2)
that attributed to the polymer will take on the form:

Ξ (p) [ζ 0 ] = lim N ′ →∞,N ′ /V fixed N ′ ! (V Q (m) ) N ′ Q tot N ′ (3) = N ≥1 ζ N 0 Q (p) N , where ζ 0 = (N ′ /V )/Q (m) . (4)
The obvious analogy with a fluid grand partition function (with the weight 1/N ! excised) will be very useful indeed.

The computation of Ξ (p) (ζ 0 ), which generates all thermodynamics and expectations in a thermal ensemble, is of course too general to be explicitely solvable, except in very special circumstances. Let us therefore start with the evaluation of Ξ (p) (ζ 0 ) for what may be regarded as the backbone of the polymer, that in which the arbitrary pair mutual interaction φ(r i , r j ) is set equal to zero. In this case, we have at once:

ζ N 0 Q (p) N = ζ(r 1 )w(r 1 , r 2 )ζ(r 2 ) . . . ζ(r N -1 )w(r N -1 , r N )ζ(r N )dr N (5)
where r) or regarding w(r, r ′ ) as the kernel of an integral operator w, ζ(r)δ(r -r ′ ) as that of the diagonal operator ζ, and the symbol |1 denoting the vector whose components are 1,

ζ(r) = ζ 0 e -βu(r) = e βµ(
Ξ (p) [ζ] = 1|(ζ -1 -w) -1 |1 , (6)
subject of course to convergence of the series [START_REF] De Gennes | Scaling concepts in polymer physics[END_REF]. The corresponding "grand potential" is:

Ω (p) [ζ] = - 1 β ln 1|(ζ -1 -w) -1 |1 (7)
and as an immediate consequence the monomer density is given by

n(r) = -βζ(r) δ δζ(r) Ω (p) [ζ] (8) = 1|(ζ -1 -w) -1 |r r|(ζ -1 -w) -1 |1 ζ -1 (r)/ 1|(ζ -1 -w) -1 |1 (9)
The absence of the statistical weight 1/N ! in (4) suggests an unusually broad distribution of monomer number N about its mean, in the ensemble. In Appendix A, we examine the form of this distribution, its consequences, and the nature of required correction terms.

B. Minimum Principle

Expression (7) for the backbone grand potential is very concise, but it requires the inversion of an integral operator, ζ (-1) -w, and so is technically non-trivial, although computationally not a great impediment. More importantly, it makes unfeasible the direct use of (7) as potential building block for the full interacting system. Fortunately, (7) can be constructed, making use of a minimum principle in which no such inverses appear. Doing so depends upon the fact that if K is positive semidefinite, and a and ψ are arbitrary, then according to the Schwartz inequality we have

a|K -1 |a ψ|K|ψ = K -1/2 a|K -1/2 a K 1/2 ψ|K 1/2 ψ (10) ≥ K -1/2 a|K 1/2 ψ 2 = a|ψ 2 , (11)
βΩ (p) [ζ] = M in ψ ln ψ|ζ -1 -w|ψ -2 ln 1|ψ . (13)
The representation (13), it must be emphasized, is exact. There is no explicit inverse contained therein, and this indeed allows it to be used as a first step in analysis of the backbone with non-neighbor pair forces.

III. THE MEAN FIELD STRATEGY A. Exact Results

Having the backbone of the polymer in an external field potentially under control, we now approach our real objective: to find the effect of mutual interactions, φ 1 (r, r ′ ), between any two monomers at r and r ′ . The most common leading estimate in such situations, termed mean field, is one in which the mean value of the sum of pair interactions acting on a given unit is taken as an additional external field. This acts self-consistently on each unit, and so control of the noninteracting system then suffices. For classical fluids, Widom's insertion theorem [9] tells us that if the calculation is done instead for each member of the ensemble -fluctuations hence automatically included -and then averaged over the ensemble, it becomes exact in the sense that n(r) = e β(µ(r)-φ 1 (r,r ′ )n(r ′ )dr ′ ) (

where

n(r) = i δ(r -r i ) ( 15
)
is the microscopic, per configuration, particle density. Mean field thus amounts to replacing n(r ′ ) by n(r ′ ) = n(r ′ ). Carrying out corrections however leads to questions about the pair distribution, which can get quite involved. For our purposes, another exact representation is preferable, in which the interaction is replaced by a suitable average over an equivalent fluctuating external field. This is of very general validity, but it applies rigorously only to the additions of purely attractive (negative definite) pair forces, which certainly can distort local properties of a real system. In the case of a polymer backbone, which already includes next neighbor pair forces, this assumption is not lethal, and we can anticipate a reasonable representation of global structure. Suppose then that, quite generally, we set for convenience (φ is now positive definite as a continuous matrix) and imagine that φ is added to a reference system all of whose properties are known.

φ 1 (r, r ′ ) = -φ(r, r ′ ) ( 16
Ξ[µ, φ] = N . . . W (N ) 0 (r N , µ)e β 2 ′ i,j φ(r i ,r j ) dr N (17) = N . . . W (N ) 0 (r N , µ)e -β 2 i φ D (r i ) e β 2 i,j φ(r i ,r j ) dr N (18) = N . . . W (N ) 0 (r N , µ)e -β 2 i φ D (r i ) e β 2 n(r)φ(r,r ′ )n(r ′ )drdr ′ dr N (19)
where ′ omits the i = j contribution, φ D is the diagonal part of the matrix φ(r i , r j ).

The device of Kac, Siegert, Hubbard and Stratonovich [START_REF] Kac | Applied Probability[END_REF][START_REF] Siegert | [END_REF][12][13] is to represent the Gaussian in [START_REF] Poland | Theory of Helix-Coil Transitions in Biopolymers[END_REF] (in obvious notation) as a functional Laplace transform

e β 2 n•φn = e -β 2 v•φ -1 v e -βv•n Dv / e -β 2 v•φ -1 v Dv (20) = e -β i v(r i ) e -β 2 v•φ -1 v Dv / e -β 2 v•φ -1 v Dv (21) Since N . . . W (N ) 0 (r N , µ)dr N = Ξ 0 [µ],
eq.( 17) can thereby be rewritten as

Ξ[µ, φ] = Ξ 0 [µ - 1 2 φ D -v]e -β 2 v•φ -1 v Dv / e -β 2 v•φ -1 v Dv (22)
with a possible interpretation that the interaction -φ has been replaced by an ensemble average over a fluctuating external field v(r), serving as a sort of graviton shuttling back and forth between units, and in fact, particle configuration space has now been transformed to field space. The physical interpretation is comfortable, but more importatntly, it leads to a new method of evaluation as well, and applies to fluids, polymers,.... The kernel of ( 22) is a Boltzmann factor in field space, and so we may define a field average as

G[v] = G[v]Ξ 0 [µ - 1 2 φ D -v]e -β 2 v•φ -1 v Dv /Ξ[µ, φ] (23)
A very suggestive consequence of this notation follows from the observation that for the density

n(r), n(r) = 1 Ξ[µ, φ] δ δβµ(r) Ξ[µ, φ] (24)
= - δ δβµ(r) Ξ 0 [µ - 1 2 φ D -v]e -β 2 v•φ -1 v Dv /Ξ[µ, φ] e -β 2 v•φ -1 v Dv, (25)
or integrating by parts in v-space (assuming the absence of boundary terms) identifying -φn(r) as the "mean field" v(r) . Eq(27) is exact, but depends very much on the details of the system under study e.g. in (26). Exact calculations in function space are out of the question, and so it is the approximate evaluation of of v(r) that we must attend to.

n(r) = -φ -1 v(r)Ξ 0 [µ - 1 2 φ D -v]e -β 2 v•φ -1 v Dv /Ξ[µ, φ] e -β 2 v•φ -1 v Dv (26) = -φ -1 v(r) , (27)

B. A Mean Field Approximation

There are any number of ways of setting up approximation sequences for the evaluation of v(r) . We will adopt one which is fairly obvious in the context of ( 22) and ( 26), with the great advantage of a very pictorial leading order -still appearing as a minimum principle -and a physically transparent first order. It goes as follows. If the field Boltzmann factor of ( 23) is sharply peaked about a function v(r), than of course (27) becomes simply

n(r) = -φ -1 v(r), (28)
and the field v must now satisfy (Ω is again the grand potential

-1 β ln Ξ) δ δv(r) Ω 0 [µ - 1 2 φ D -v] + 1 2 v • φ -1 v = 0, (29)
or simply

n 0 (r|µ - 1 2 φ D -v) = -φ -1 v(r), (30)
the unsurprising result (comparing with (28)) that n(r) is equal to the "bare" density n 0 in the presence of v, to within an additional shift of 1 2 φ D . But in addition, we now have in leading order approximation the very explicit

Ω 1 (µ, φ) = Ω 0 [µ - 1 2 φ D -v] + 1 2 v • φ -1 v (31)
which, since (29) represents the minimization of the kernel of ( 24), takes the variational form

Ω 1 (µ, φ) = Min v Ω 0 [µ - 1 2 φ D -v] + 1 2 v • φ -1 v. (32)
Furthermore, consistency is established by noting that by virtue of (32), Ω 1 of (31) implies

n(r) = - δΩ 1 δµ(r) = -φ -1 v(r), (33)
reproducing (28).

To apply (32) to the single homopolymer under discussion, we need only insert (13), serving as Ω 0 , into (32). We then have where µ ′ (r) = µ(r) -1 2 φ D (r). Eq. ( 34) can be simplified by finding and eliminating the field ψ, and replacing v by -φn. To do so we observe that

βΩ 1 [µ, φ] = Min v,ψ β 2 v • φ -1 v -2 ln 1|ψ + ln ψ|e -β(µ ′ -v) -w|ψ (
0 = δΩ 1 δv = φ -1 v + e -β(µ ′ -v) ψ 2 / ψ|e -β(µ ′ -v) -w|ψ (35)
Integrating over the implicit r, with N = n(r)dr, then

N = ψ|e -β(µ ′ -v) |ψ ψ|e -β(µ ′ -v) -w|ψ , (36)
reducing (35) to

n = N e -β(µ ′ -v) ψ 2 / ψ|e -β(µ ′ -v) |ψ (37)
Since ( 37) is homogeneous of degree 0 in ψ, we are free to adopt the normalization

ψ|e -β(µ ′ -v) |ψ = N (38) converting (37) to n = e -β(µ ′ -v) ψ 2 (39)
(leading back to (38)), and consequently replace (34) by the simple

βΩ 1 [µ, φ] = Min n β 2 n • φn -2 ln 1|n 1/2 e β 2 (µ ′ +φ•n) + ln n 1/2 | I -e β 2 (µ ′ +φ•n) we β 2 (µ ′ +φ•n) |n 1/2 . ( 40
)
appearing as a µ-dependent density functional, reminiscent of "statistical models" [START_REF] Giacomin | Random Polymer Models[END_REF][START_REF] Flory | Statistical mechanics of chain molecules[END_REF][START_REF] Huggins | [END_REF], of polymers. Eq.( 40) is our main result, valid in the mean field level. The role of the interaction -φ is transparent: the external potential is augmented by the mean field :µ → µ ′ + φn, and the explicit energetic component -1 2 n • φn subtracted out. Eq.( 40) is of course a density functional representation with the density profile dependence on µ determined by dropping the implicit φ and making µ explicit, by

δΩ 1 [µ, n] /δn(r) = 0. ( 41
)
By definition, see Eq.( 40), we have ignored the fluctuations of the potential field v(r) about its mean v(r). Corrections must then correspond to including the fluctuations to some model order. This is analogous to including capillary waves in the description of a two-fluid interface, and they do make both energetic and entropic contributions. In Appendix B, we show briefly how to obtain the effect of potential field fluctuations to Gaussian order, but leave resulting numerical contributions to the future.

IV. NUMERICAL IMPLEMENTATION

The minimization required in (40) to obtain numerically the profile n(r) can be carried out in a number of ways. However, in many-monomer systems, which are of major interest, there is an important preliminary observation. In Appendix C, we analyze the system behavior as the monomer number increases at fixed external confinement potential. We show that, unsurprisingly considering the floppiness of the chains we are working with at this introductory stage, increasing mean monomer number serves only to increase the monomer density uniformly over the confinement volume. It is therefore only the normalized monomer density p(r) = n(r)/N that is relevant.

In terms of p the minimization functional (40) becomes

Ω 1 [µ, φ] = Min p N 2 2 p • φp -2 ln 1|p 1/2 e β 2 (µ ′ +N φψ) + ln p 1/2 | I -e β 2 (µ ′ +N φψ) we β 2 (µ ′ +N φψ) |p 1/2 . ( 42
)
Once we have determined p(r), the expression to compute the particle number and the true density are

n(r) = p(r) p(r) 1/2 [I -ζwζ]p(r ′ ) 1/2
(43)

N = p(r) 1/2 [I -ζwζ]p(r ′ ) 1/2 -1 (44)
For the numerical implementation, we considered a chain confined by hard walls, with harmonic φ is turned on the distribution becomes more peaked (figure 2).

We also looked at the effect of introducing an asymmetry in the external potential: a uniform , gravitational, potential along one of the axes. As expected the minimization of the functional converges to a probability density off-centered with respect to z and centered in the middle of the box with respect to x and y (figure 3).

Even if the theoretical derivation of eq.( 40) is valid only for strictly positive definite φ, in the algorithm we can invert the sign of the long range interaction and make it repulsive, instead of just attractive. We can therefore test our method beyond the known validity of the approximations used, taking a repulsive long-range interaction. The method still holds, provided we consider repulsive interactions which are not too strong. As expected a repulsive φ gives a less peaked density compared to the reference system (figure 4). The monomer interections we consider are both repulsive and attractive. Under this condition we could expect the polymer to undergo a coil-globule transition as the strength of the attractive interaction overcomes the repulsion [START_REF] Poland | Theory of Helix-Coil Transitions in Biopolymers[END_REF].

In our system, when no long-range interaction is present, the monomers interact harmonically, and the polymer configuration should be that of an extended coil, with the monomers sitting at their equilibrium separations. When the long-range interction φ is turned on the monomers start feeling a pure attraction, to all other monomers. There will never be a balance between attraction and repulsion because the repulsion acts only between nearest neighbors, while the attraction is to all other particles in the system. On energetic grounds alone, we would then expect our system to collaps to a globular conformation as soon as the long-range attraction is introduced, even when this is very weak. However, the increased entropy of an open structure softens the rigidity of this statement. Through direct Monte Carlo simulations of a model polymer inside a spherical cavity we indeed verified this behavior for an attractive φ. When repulsive long-range interactions are introduced instead, the polymer becomes somehow rigid and tends to wrap around in the cavity using all avaliable space to become as close as possible to a straight rod.

With the method developed in this paper, we have systematically looked at the effect of turning on the long range interaction, performing minimizations at regular φ increments. In Figure 5 we plot one of the coefficient of the exponent of the gaussian density (we have chosen B x , but results are independent of this choice) as a function of φ. We observe a highly non-linear dependence, indicating a compactification of the system even beyond the coil-globule transition. To observe a realistic coil-globule transition the long-range interaction itself should have a repulsive core and an attractive tail.

When looking at the details of the polymer configuration, proper quantities to measure are the end-to-end distance or the radius of gyration. At the present stage we measure the average density only. Other quantities are accessible with our formulation, but their theoretical treatment, though straight forward, is lengthy and will involve completely new numerical implementations. This is beyond the scope of this paper, where we want to establish the basic setting for the general formulation of confined polymer behavior in terms of a minimization principle and of mean field.

V. CONCLUDING REMARKS

The simplifying assumptions attend more directly to the physics, and these assumptions depend very much on the nature of the system to be studied. Taking these assumptions in order, we first For a Gaussian density like the one we have, a measure of B x gives an idication of how peaked the density distribution is. We observe a strong non-linear dependence, indicating a rapid compactification of the polymer as φ increases.

emebedded out system in a monomer grand ensemble. Since preprocessing assures that in practice one does not deal with the resulting extreme polydispersity, a fixed N ensemble is more relevant than fixed ζ. The corresponding inverse mapping has been attended to on numerous occasions (see e.g. [START_REF] Lebowitz | [END_REF]). The same formalism is indeed avaliable here (hinted at in eqs.(A15),(A16)).

Another restriction was to attractive long-range forces. We found however that the profile equation could indeed be pushed into the partially repulsive regime, although the validity of the minimum principle was in question; this is closely related to the functional Fourier transform for the repulsive component, likewise under uncertain control. An alternative approach lies in the use of the mean spherical model [19,20] and its extensions. This is the aim of ongoing research.

Of course, there is the implicit assumption that pair forces suffice, whereas the action of pairs on singlets is a frequent important observation, leading e.g. to dihedral angle dependence in chains.

Typically (see e.g. [21] for a very primitive example) one can simply create a multi-unit monomer to encompass only such forces, which then appear once more as pair forces.

Most blatantly, our approach has been restricted to homopolymers. Since the set of degrees of freedom of a monomer can include monomer type, this is no restriction at all if one is studying the effect on the full population of an assumed relative frequency of occurrence of next neightbors hetero-pairs. However, if we attribute a sense to the chain and the AB frequency differs from the BA frequency, the very convenient symmetry of the operator w is lost, and with it, the possibility

βF N -βΩ (p) = -ln Q N + ln Ξ (p) [ζ] (A13) ∼ = N ln ζ 0 -(N -1) ln λ 0 -ln(1 -λ 0 ). (A14)
It follows that for any (N -independent) parameter variation

β(δF N -δΩ (p) ) = 1 1 -λ 0 - N -1 λ 0 δλ (A15) which vanishes if ζ 0 in (A5) is chosen so that N = 1/(1 -λ 0 ) (A16)
We conclude that for this choice of λ 0 , expectations at fixed N and in the monomer number ensemble are in fact identical to leading order. This is not a surprise. Quite generally, for fluids, (see e.g [20]) inversion of the transformation from an N -ensemble to µ-ensemble follows the direct prescription above to leading order, with a correction term for expectations going as (δN ) 2 /N 2

. For fluids, this indeed goes as 1/N , but for a chain, using (19), as N 0 , and is therefore not a correction that vanishes as N → ∞. Another means of assessment is mandatory. In (12), we see that the relative correction to N has an amplitude ∼ (1 -λ 0 )/(1 -λ 1 ); the same is true for other expectations. Now suppose that, at some value of the system chemical potential µ, one has the pair λ 0 , λ 1 . Then if µ is increased by δµ, ζ (1/2) wζ (1/2) → exp(βδµ)ζ (1/2) wζ (1/2) , so that λ 0 → exp(βδµ)λ 0 , λ 1 → exp(βδµ)λ 1 . Then, as exp(βδµ)λ 0 is made arbitrarily close to 1, exp(βδµ)λ 1 will never go beyond λ 1 /λ 0 , and the ratio (1 -λ 0 )/(1 -λ 1 ) indeed goes as 1/N .

APPENDIX B: FLUCTUATIONS TO GAUSSIAN ORDER

We now take a first step towards including the effect of field fluctuations on the molecular density pattern. In this inital study, we have made one major approximation and several simplifying restrictions. The approximation is of course that of mean field, or selective neglect of fluctuations.

On the assumption that fluctuations are Gaussian to leading order (examples in which this is not the case are far from rare, see e.g. [23]), a correction sequence is in primciple routine: we expand ln Ξ 0 in ( 22) and ( 27) around v of (29). Using It readily follows (see e.g. [24]) that the density profile is given to the next order by simply averaging the density over the Gaussian field fluctuation:

δ 2 Ω 0 /δv(r)δv(r ′ ) = δn 0 (r)/δv(r ′ )| v, Ξ 0 [µ - 1 2 φ D -v] exp{- 1 2 β v • φ -1 v} = = exp -β Ω 0 [µ - 1 2 φ D -v] + 1 2 v • φ -1 v + 1 2 (v -v) • δn 0 δv | v + φ -1 (v -v) + . . . (
n(r) = n 0 [µ -1 2 φ D -v -∆]e -1 2 ∆ 2 M -1 (r,r) d∆ e -1 2 ∆ 2 M -1 (r,r) d∆ (B2)
where

M (r, r ′ ) = φ -1 (r, r ′ ) -n ′ 0 [µ -1 2 φ D -v -∆]
; here, ∆ is the field amplitude fluctuation at r.

APPENDIX C: LIMIT OF TIGHT CONFINEMENT

We have seen, in Sec. II B, that in the absence of non-neighbor interactions, the large N monomer distribution is determined by the "resonance state" ψ λ 0 satisfying

(ζ 1/2 wζ 1/2 )ψ λ 0 = λ 0 ψ λ 0 , (C1)
where

λ M AX = λ 0 ∼ 1. Explicitely, since r| ζ -1 -w -1 |1 = r|ζ 1/2 I -ζ 1/2 wζ 1/2 -1 ζ 1/2 |1 ∼ 1 1 -λ 0 r|ζ 1/2 ψ λ 0 ψ λ 0 |ζ 1/2 |1 (C2)
we have under these circumstances, from ( 9)

n(r) ∼ 1 1 -λ 0 ψ λ 0 (r) 2 , ( C3
)
or since ψ λ 0 (r) ∼ ψ 1 (r),

n(r) ∼ ψ(r) 2 1 -λ 0 = N ψ 1 (r) 2 (C4)
In other words, increasing the chain length increases the mean density uniformly, as if the floppy chain is simply winding around more under the same confinement. This uninteresting behavior is refined by two interactions that have been ignored in getting (C4). First is the stiffness of successive pair orientations, equivalent to the monomers having coupled orientational degrees of freedom, a topic that has been addressed to some extent in this format in the past and will be attended to more forcefully in the future. Second is the effect of non-next neighbor interactions, which we have here studied in a preliminary fashion. The mean field v(r) that has been enountered will also

Ω 1 (ζ, φ) = Min v,ψ [ 1 2 v • φ -1 v -2 ln (ρe -βv ) 1/2 |(ρ -1 e βv ) 1/2 ψ + ln (ρ -1 e βv ) 1/2 ψ|I -(ρe -βv ) 1/2 w(ρe -βv ) 1/2 |(ρ -1 e βv ) 1/2 ψ ] = Min v, ψ 1 2v • φ -1 -2 ln (ρe -bv ) 1/2 | ψ + ln ψ|I -(ρe -bv ) 1/2 w(ρe -βv ) 1/2 | ψ (C5)
Hence, according to (38), ψ = (ζ -1 e βv ) 1/2 ψ is normalized, ψ| ψ = 1 and from (39),

n = N ψ2 (C6)
Near resonance is now signaled by the approximate validity of

I -(ζe -βv ) 1/2 w(ζe -βv ) 1/2 ψ = λ ψ (C7)
with small λ > 0. Eq. ( 50) then tells us at once that using the exact

φ -1 v = -n, (C8)
we have

N = 1/ λ (C9)
There are two ways of making use of the approximation C7. Most directly, we substitute (C6-C9) into (C5), obtaining

βΩ 1 ∼ Min n 1 2 n • φn -2 ln (ζe βφn ) 1/2 |n 1/2 (C10)
This is obviously too sweeping an approximation: by using what is effectively a first order correction in the argument of a variational principle, the role of the next-neighbor w in βΩ 1 has vanished.

But we can pick up the next order by working instead at the "profile equation" level. It is only necessary to substitute (C6, C8, C9) directly into (C7) to rewrite the latter as

(ζe βv ) 1/2 w(ζe -βv ) 1/2 n 1/2 = 1 - 1 N n 1/2 , (C11)
and hence as From preliminary simulation studies we had noticed how in the absence of a repulsive core interaction between monomers we could not converge to a stationary density. Instead the density, and therefore the total particle number, either kept on growing or shrank to zero, depending on the value of the external chemical potential. This numerical instability occures as a result of the closeness to the resonance. In these simulations we noticed, on the other hand, that despite the fact that the density itself kept changing, the global shape of density profile reached a steady state.

e -β 2 (u-φn) we β 2 (u-φn) n 1/2 = 1 - 1 N e -βµ ′ 0 n 1/2 ≡ Λn 1/2 (C12) (µ ′ 0 is
We therefore decided to study the normalized density p(r) = n(r)/N . Once we had formulated the problem consistently in terms of p all the numerical instabilites desappeared and we were able to determine the corresponding stady state densities and total particle numbers. The minimization is done using a Metropolis Monte Carlo algorithm at fixed temperature, starting from a suitable ansatz for the functional form of p(r). This is one of many possible techniques, and is not optimized in any well-defined sense. We can freely impose the form of the next neighbor potential, and therefore w(r, r ′ ), the external confining potential u(r), and the long range interaction φ(r, r ′ ). The main numerical difficulty arises from the need of evaluating a functional involving 6-dimensional integrals nesting 3-dimensional integrals, such as p(r)e β 2 (µ ′ (r)+N φ(r,r ′′ )p(r ′′ )dr ′′ ) w(r, r ′ )e β 2 (µ ′ (r ′ )+N φ(r ′ ,r ′′ )p(r ′′ )dr ′′ ) p(r ′ )dr dr ′ (D1)

It turns out it is prohibitive to try to evaluate these terms directly by the use of a simple grid.

Already when using 10 discretization points simulations are too slow, while the precision is very poor. To achieve higher accuracy in the integration, we use a Legendre-Gauss quadrature method with either 6 or 8 points in each linear direction.

As a test system we have chosen a square well external potential with µ = -µ 0 inside; a harmonic n.n. potential, leading to w(r, r ′ ) = exp(-C h (|rr ′ | -ρ 1 ) 2 ); a long-range step potential φ = C lr Θ(x -ρ 2 ), where µ 0 , C h , C lr are energy constants, and ρ 1 and ρ 2 are characteristic lengths.

In the following we take the system size to be L = 1, C h = 5, and ρ 1 = 0. probability density we take the composition of three Gaussian in the three spatial directions, for a total of 6 parameters: B x ,B y ,B z ,X 0 ,Y 0 and Z 0 .

p(x, y, z) = A e -Bx(x-X 0 ) 2 e -By(y-Y 0 ) 2 e -Bz(z-Z 0 ) 2 (D2)

At each iteration of the algorithm the normalization of n is verified so that p(x, y, z)dx dy dz = 1, and the value of A is changed accordingly.

FIG. 1 :

 1 FIG. 1: Graphical sketch of the monomer bath. The chain is in grand canonical equilibrium with a bath of inert monomers.

Page 6 of 25 URL

 25

Page 8 of 25 URLFIG. 2 :

 252 FIG. 2: Left: Interaction between neighboring particles along the chain: harmonic potential only and harmonic potential plus attractive interaction. Right: Slices of the probability density profiles taken at y = z = L/2.

FIG. 3 :FIG. 4 :

 34 FIG.3: Probability density profile at z = L/2 in the presence of an external graviational field.

Page 10 of 25 URL

 25

Page 11 of 25 URLFIG. 5 :

 255 FIG. 5: The coefficient B x measured in terms of the strenght of the long-range attractive interaction φ.

  . of a specific long sequence of monomer types. The case of non-symmetric w has indeed been studied[22], but exercising the kind of control that we have in our current formulation remains a challenge. that if F N is the canonical Helmholtz free energy,

  have the effect of correlating pair orientations, and will of course alter the nature of the long chain resonant state. Let us see how this works. To keep the extrapolation from (C1-C4) transparent, us rewrite (34) (using ζ = e βµ ′ ) as,

  the global chemical potential). The operator in (C12) has a positive kernel, hence by Jentchke's extension of Perron-Frobenius[25], the positive eigenfunction is unique (up to a multiplicative constant) and Λ is real and maximal. This leads to what is in principle a simple iteration: start eg. with φ = 0, compute the eigenfunction n 1/2 and normalize it to get n(r)dr = N . At the k th stage of the M-fold iteration, replace φ by (k/M )φ and repeat the process using the current function n. An alternative strategy is to parameterize n and determine the parameters by Galerkin, i.e. integrate (C12) with weight function and solve the resulting algebraic equation.APPENDIX D: DETAILS ON NUMERICS

2 .

 2 As ansatz for the

Page 18 of 25 URL

 25 ://mc.manuscriptcentral.com/tandf/://mc.manuscriptcentral.com/tandf/://mc.manuscriptcentral.com/tandf/://mc.manuscriptcentral.com/tandf/://mc.manuscriptcentral.com/tandf/tmph Molecular Physics F o r P e e r R e v i e w O n l y Page 25 of 25 URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Page 4 of 25 URL: http://mc.manuscriptcentral.com/tandf/tmph

			Molecular Physics
	F o r
	P
	e
		e r
			R
			e
			v i e
				w
				O n l
					y
	so that a|K -1 |a ≥ a|ψ 2 / < ψ|K|ψ . Then indeed
	Max ψ	a|ψ 2 ψ|K|ψ	= a|K -1 |a	atKψ =	ψ|K|ψ a|ψ	a,	(12)

Appendix A: Number Distribution

We first observe from (9) that if N denotes the monomer number in a given configuration, than

where

) and λ 1 is the next largest eigenvalue; λ 1 < λ 0 in a confined system.

If N is large, then

is dominated by the "resonance" at λ 0 , so that if

with normalized ψ, then

There are now two consequences. On the one hand, we have

so that the distribution of N is given by

the very broad geometric distribution. It would appear that the N -ensemble must give a very poor representation of a given N . But on the other hand, we have from ( 5)


Maria Isabel

Garcia De Soria

Pablo Maynar

Dr Jean-Jacques Weis

María Isabel García De Soria

Emmanuel Trizac

(SI)Energy fluctuations in a randomly driven granular fluid

Keywords: numbers: 51.10.+y, 05.20.Dd, 82.20.Nk

come

I. INTRODUCTION

Recent years have witnessed ongoing interest for the microscopic and macroscopic properties of granular media. In such systems, a simple ingredient -energy dissipation resulting from collisions-has far reaching consequences [START_REF] Barrat | [END_REF], with rich phenomenology: non-Gaussian velocity distributions [2,3], non-equipartition of energy [4][5][6][7], or spontaneous symmetry breaking [8][9][10] to name but a few. Theoretically, one of the tools used to understand this body of phenomena is kinetic theory, which is extended naturally to these systems by introducing an "inelastic collision rule" in which the energy is not conserved. Although most of the work carried out until now has focused on the one-particle properties, and on the study of the corresponding Boltzmann equation, it has been shown that correlations are also important and, as a matter of fact, necessary to understand the behavior of the system when vortices or cluster are developed [11,12], or even in simpler situations where spatial homogeneity is enforced [13,14].

As a consequence of the dissipation in collisions, the total energy of an isolated granular system decays monotonically in time. Under certain conditions, the system reaches a homogeneous cooling state in which the time dependence of the one-particle distribution function is entirely embodied in the kinetic (so-called granular) temperature, which evolves with time as t -2 [2,15]. It is experimentally difficult to probe such a regime (see however [16,17]). Nevertheless, it is possible to maintain a granular system in the fast-flow regime by injecting energy in such a way that a stationary state is reached. In these states, the energy injected by the thermostat is compensated by the energy dissipated in collisions. Several mechanisms can be introduced in order to get a stationary state. If, for example, energy is injected by a moving boundary such as a vibrating piston, the system reaches an inhomogeneous stationary state [18]. In this work, we will focus on a granular gas which is heated uniformly, coupling the velocity of each particle to a white noise, this is the so-called "stochastic thermostat" [3,[19][20][21][22][23][24][25][26][27]. For this kind of forcing, which is relevant for some two-dimensional experimental configurations with a rough vibrating piston [28], the system reaches a homogeneous stationary state after a transient regime. The advantage of such a driving mechanism is that it lends itself to theoretical progress. In this context, the single particle distribution function has been characterised [3] and long range correlations predicted from a hydrodynamic treatment [21]. More recently, the fluctuations of the total energy have been analyzed [29] (see also [30] for a related numerical study in an inhomogeneous system). In Ref. [29], the second moment of the total energy fluctuations was evaluated by neglecting the correlations, which, by and large, could not explain the simulation results. The objective of the present work is to clarify and quantify the influence of the inelasticity induced correlations on the total energy fluctuations. The methods used bear some similarities with those reported in [31], where it was shown that for the unforced system, the contribution coming from the correlations is of the same order as that coming from the one-particle distribution function itself.

The paper is organized as follows. In section II, the equation for the two-particle distribution function is derived, taking due account of the thermostat while in section III, the results are particularized to the homogeneous stationary state, that will play the role of our reference state in subsequent analysis. There we also summarize the main results already known and pertaining to the one-particle distribution function. In section IV, we analyze the hydrodynamic equations for a homogeneous linear perturbation of the reference state, and obtain the corresponding modes and eigenvalues, which are finally used in section V to meet our objective and obtain an explicit expression for the variance of the total energy. We consider a gas of N hard disks (dimension d = 2) or spheres (d = 3) of mass m and diameter σ that collide inelastically with a coefficient of normal restitution α [START_REF] Barrat | [END_REF]. The system is heated uniformly by adding a random component to the velocity of each particle at equal times [21,22]. The driving is implemented in such a way that the time between random kicks is small compared to the mean free time. Then, between collisions, the velocities of the particles undergo a large number of kicks due to the thermostat. In addition, we will assume that the "jump moments" of the velocities of the particles verify

B ij,βγ ≡ lim ∆t→0 ∆v i,β ∆v j,γ ∆t = ξ 2 0 δ ij δ βγ + ξ 2 0 N (δ ij -1)δ βγ , (1)
i, j = 1, . . . , N and β, γ = 1, . . . , d where we have introduced ∆v i,β ≡ v i,β (t + ∆t) -v i,β (t), v i,β (t) being the β component of the velocity of particle i at time t. We have also introduced the strength of the noise, ξ 2 0 , and . . . , which denotes average over different realizations of the noise. The non-diagonal terms (corresponding to i = j and β = γ) are necessary in order to conserve the total momentum.

In the dilute limit, assuming molecular chaos, i.e. that no correlations exist between colliding particles, and that the sizes of the jumps due to the thermostat are small compared to the velocity scale on which the distribution varies, the equation for the single particle distribution function in our system is the Boltzmann-Fokker-Planck equation [3,22]

∂ ∂t f (x 1 , t) + L (0) (x 1 )f (x 1 , t) = J[f |f ] + ξ 2 0 2 ∂ ∂v 1 2 f (x 1 , t), (2)
where x i is a short-hand for position-momenta coordinates {r i , v i } and

L (0) (x 1 ) = v 1 • ∂ ∂r 1 . (3)
The inelastic collision operator J[f |f ] reads

J[f |f ] = σ d-1 dv 2 T0 (v 1 , v 2 )f 1 (r, v 1 , t)f 1 (r, v 2 , t), (4)
where

T0 (v 1 , v 2 ) = dσΘ(σ • g)(σ • g) α -2 b -1 σ -1 , (5)
with g = v 1 -v 2 the relative velocity, Θ the Heaviside step function, σ a unit vector joining the centers of the particles at contact and b -1 σ an operator replacing the velocities v 1 and v 2 appearing on its right by the precollisional values

v * 1 ≡ b -1 σ v 1 = v 1 - 1 + α 2α (g • σ)σ, (6)
v * 2 ≡ b -1 σ v 2 = v 2 + 1 + α 2α (g • σ)σ. (7)
The term

ξ 2 0 2 ∂ ∂v1 2 f (x 1 , t
) is a diffusive Fokker-Plank term, and is a signature of the external noise. As we shall study fluctuations, it is convenient to introduce the two-particle distribution function, f 2 (x 1 , x 2 , t). The quantity f 2 (x 1 , x 2 , t)dx 1 dx 2 is defined as the number of pairs of particles in which one lies inside the differential volume dx 1 centred in x 1 and likewise, with dx 2 , x 2 for the second particle. This definition is easily generalized to higher n-particle distribution functions, f n (x 1 , . . . , x n ). The evolution equation for f 2 is [31,32]

∂ ∂t + L (0) (x 1 ) + L (0) (x 2 ) f 2 (x 1 , x 2 , t) = δ(r 12 )σ d-1 T0 (v 1 , v 2 )f 2 (x 1 , x 2 , t) +σ d-1 dx 3 δ(r 13 ) T0 (v 1 , v 3 ) + δ(r 23 ) T0 (v 2 , v 3 ) f 3 (x 1 , x 2 , x 3 , t) + F T H , (8)
where we have introduced F T H that accounts for the external driving. The evolution equation (8) contains essentially three parts: the free streaming in the left-hand side, the two terms in the right hand side corresponding to collisions, 3 and the last term, F T H due to the thermostat. The collisional contribution is split in one part corresponding to collisions of particles with velocities v 1 and v 2 , and the other which involves collisions of particles with velocities v 1 or v 2 with a third particle with arbitrary velocity, v 3 . The collisional contribution is identical to the one that appears in the absence of forcing [31]. We concentrate now on the new term, F T H . Assuming that the sizes of the jumps due to the thermostat are small compared to the scale in which the distribution f 2 varies, we can expand F T H in the spirit of the Fokker-Planck description [START_REF] Van Kampen | Stochastic Proccesses in Physics and Chemistry[END_REF]

F T H ≃ 1 2 d β,γ=1 2 i,j=1 B ij,βγ ∂ ∂v i,β ∂ ∂v j,γ f 2 (x 1 , x 2 , t) = 1 2 ξ 2 0 ∂ 2 ∂v 2 1 + ∂ 2 ∂v 2 2 - 2 N ∂ ∂v 1 • ∂ ∂v 2 f 2 (x 1 , x 2 , t), (9)
where we have taken into account equation ( 1), and we have explicitly assumed that the jump moments B ij,βγ do not depend on the magnitude of the velocities of the particles.

Let us introduce the two-particle and three-particle correlation functions through the usual cluster expansion

f 2 (x 1 , x 2 , t) = f 1 (x 1 , t)f 1 (x 2 , t) + g 2 (x 1 , x 2 , t), (10)
and

f 3 (x 1 , x 2 , x 3 , t) = f 1 (x 1 , t)f 1 (x 2 , t)f 1 (x 3 , t) + g 2 (x 1 , x 2 , t)f 1 (x 3 , t) +g 2 (x 1 , x 3 , t)f 1 (x 2 , t) + g 2 (x 2 , x 3 , t)f 1 (x 1 , t) + g 3 (x 1 , x 2 , x 3 , t). ( 11
)
The equation for the correlation function g 2 (x 1 , x 2 , t) can be obtained following the same lines as in references [31,[START_REF] Maynar | [END_REF]. Neglecting the three-body correlations, g 3 , in Eq. ( 8), we obtain

∂ ∂t + L (0) (x 1 ) + L (0) (x 2 ) g 2 (x 1 , x 2 , t) = δ(r 12 )σ d-1 T0 (v 1 , v 2 )f 1 (x 1 , t)f 1 (x 2 , t) + [K(x 1 , t) + K(x 2 , t)] g 2 (x 1 , x 2 , t) - ξ 2 0 N ∂ ∂v 1 • ∂ ∂v 2 f 1 (x 1 , t)f 1 (x 2 , t), (12)
where we have introduced the linear operator K(x i , t) defined as

K(x i , t) = σ d-1 dx 3 δ(r i3 ) T0 (v i , v 3 )(1 + P i3 )f 1 (x 3 , t) + ξ 2 0 2 ∂ ∂v i 2 , (13)
and where the permutation operator P ab interchanges the labels of particles a and b in the quantities on which it acts.

As will become clear below, the 1/N term in equation ( 12) is crucial for the calculation of the energy fluctuations.

III. THE STATIONARY STATE

It has been shown numerically that, after a transient time, the system reaches a homogeneous stationary state [21] in which the energy input from the thermostat is compensated by the energy lost in collisions. In this section we will particularize the equations of the previous section to this state, summarizing the results that are already known about the one-particle distribution function and that are required for our theoretical analysis.

The Boltzmann-Fokker-Planck equation ( 2) for the distribution function, f H (v 1 ), in the stationary homogeneous state is

ξ 2 0 2 ∂ ∂v 1 2 f H (v 1 ) + J[f H |f H ] = 0. ( 14
)
It is convenient to introduce the scaled distribution function χ where n H is the homogeneous density, v H = 2TH m 1/2 is the thermal velocity defined from the granular temperature

H f H (v) = n H v d H χ H (c) , (15)
T H = 2 dn H dv 1 2 mv 2 f H (v), (16)
and c = v vH is the rescaled velocity. The distribution function has been studied in reference [3], where an approximate expression for χ H (c) was derived to second order in Sonine polynomials [START_REF] Landau | Physical Kinetics[END_REF] χ

H (c) = e -c 2 π d/2 1 + a 2 (α)S 2 d/2-1 (c 2 ) , (17)
with

S 2 d/2-1 (c 2 ) = 1 2 c 4 - 1 2 (d + 2)c 2 + 1 8 d(d + 2), (18)
and a 2 (α) a coefficient related to the kurtosis of the function χ H (c)

d d + 2 c 4 H c 2 2 H = 1 + a 2 (α). ( 19
)
An approximate expression for a 2 reads (see [23,[START_REF] Coppex | [END_REF] for a discussion on various possible approximations)

a 2 = 16(1 -α)(1 -2α 2 ) 73 + 56d -24αd -105α + 30(1 -α)α 2 . ( 20
)
The expression for the temperature in the first Sonine approximation is

T H = m dξ 2 0 √ π (1 -α 2 )Ω d n H σ d-1 (1 + O(a 2 )), (21)
where

Ω d = 2π d/2 /Γ(d/2) is the d-dimensional solid angle.
We now turn to the equation for the correlation function,

g 2,H (x 1 , x 2 ). It is convenient to introduce the rescaled correlation function g H through g 2,H (x 1 , x 2 ) = n H ℓ d v 2d H g H (l 12 , c 1 , c 2 ), (22)
where ℓ = (n H σ d-1 ) -1 is proportional to the mean free path and l = r/ℓ. In these units, the equation for the reduced function g H reads

Λ(c 1 ) + Λ(c 2 ) -c 12 • ∂ ∂l 12 g H (l 12 , c 1 , c 2 ) = -δ(l 12 ) T0 (c 1 , c 2 )χ H (c 1 )χ H (c 2 ) + ξ 2 0 n H ℓ d N ∂ ∂c 1 • ∂ ∂c 2 χ H (c 1 )χ H (c 2 ), (23)
where we have introduced the linearized Boltzmann-Fokker-Planck operator Λ(c)

Λ(c i )h(c i ) = dc 3 T0 (c i , c 3 )(1 + P i3 )χ H (c 3 )h(c i ) + ξ 2 0 2 ∂ ∂c i 2 h(c i ), (24)
with rescaled noise amplitude

ξ 2 0 = ξ 2 0 l v 3 H . ( 25
)
As can be seen in equation ( 23), the correlation function, g H , is determined by the properties of the linearized Boltzmann-Fokker-Planck operator, Λ, and by the one particle distribution function χ H . It is consequently important to study the spectral properties of Λ, in particular the upper (hydrodynamic) part of the spectrum, in order to understand the fluctuations of global quantities. In the case of a granular gas in the homogeneous cooling state [START_REF] Brey | Granular Gas Dynamics[END_REF], and for a system under ballistic annihilation dynamics [START_REF] García De Soria | [END_REF], it has been shown that it is possible to find the hydrodynamic eigenvalues and eigenfunctions of the linearized Boltzmann-Fokker-Planck operator. Once these quantities are known, it becomes possible to evaluate the fluctuations of the relevant global quantities in the so-called "hydrodynamic approximation" [31,[START_REF] Maynar | [END_REF]. In the remainder, we will see that we can evaluate the fluctuations of the total energy in an equivalent approximation, but without the knowledge of the eigenfunction associated to the energy. The only information needed is the form of the linearized hydrodynamic equations and, in particular, the eigenvalues. In this section we focus on the linearized hydrodynamic equations around the homogeneous stationary state. The objective is to consider the linearized equations around a homogeneous perturbation in order to extract information about the linear behavior of a small perturbation of the total energy.

The complete non-linear hydrodynamic equations for the granular system heated by the stochastic thermostat are [21,25]

∂ ∂t n = -∇ • (nu), (26)
∂ ∂t u = -u • ∇u - 1 mn ∇ j P ij , (27)
∂ ∂t T = -u • ∇T - 2 dn (∇ • q + P ij ∇ j u i ) -ζT + mξ 2 0 , (28)
where P ij is the pressure tensor, q is the heat flux and ζ is the cooling rate, which is also a functional of the distribution function

ζ = (1 -α 2 )mπ d-1 2 σ d-1 4dΓ d+3 2 nk B T dv 1 dv 2 |v 1 -v 2 | 3 f (r, v 1 , t)f (r, v 2 , t). ( 29
)
Considering a homogeneous state, the previous equations reduce to

∂ ∂t n = 0, ∂ ∂t u = 0, ∂ ∂t T = -ζT + mξ 2 0 . (30)
In the long time limit, the system is expected to approach a steady state with a constant temperature given by the equation

ζ H (f H )T H = mξ 2 0 . (31)
Substituting the explicit form of the one particle distribution function (15) in the equation above, we obtain the temperature given in equation ( 21).

Let us consider now a homogeneous state close to this homogeneous stationary state. We can write the hydrodynamic fields as n(t) = n H + δn, u(t) = δu and T (t) = T H + δT . We also define the dimensionless hydrodynamic fields

δρ(τ ) = δn n H , δw(τ ) = δu v H , δθ(τ ) = δT T H , (32)
where we have introduced the dimensionless time scale τ , proportional to the number of collisions per particle, defined as

τ = t 0 dt ′ v H ℓ = v H ℓ t. (33)
Assuming that the deviations are small, and taking into account equations ( 30)-( 31), we can write the linearized evolution equations for the dimensionless hydrodynamic fields in this new time scale

∂ ∂τ δρ = 0, ∂ ∂τ δw = 0, ∂ ∂τ δθ = -ζ 0 δρ - 3 2 ζ 0 δθ, (34)
where ζ 0 = lζH vH is a dimensionless coefficient that is a functional of the one-particle distribution function in the stationary state. Its expression in the first Sonine approximation is [3]

ζ 0 = (16 + 3a 2 )π d-1 2 (1 -α 2 ) 8 √ 2dΓ d 2 . ( 35
)
To obtain the equation for δθ we have assumed that the perturbed distribution function scales as , and χ H is the same scaled distribution function as for the reference stationary state. This assumption has already been used and tested numerically in [21]. Then, the cooling rate ζ for the state under scrutiny is proportional to T 1/2 (t) and we obtain the equation for the linearized energy written above. Equations [START_REF] Maynar | [END_REF] indicate that a perturbation in the total number of particles or total momentum does not decay, as a consequence of the fact that these variables are conserved, but a perturbation in the total energy will decay (exponentially in τ ) to the stationary value, as expected. Moreover, as the equation for the temperature can be rewritten in the following form

f (v, t) = n v(t) d χ H v v(t) , (36)
∂ ∂τ 2 3 δρ + δθ = - 3 2 ζ 0 2 3 δρ + δθ , (37)
we can identify the hydrodynamic eigenvalues λ = 0 and γ = -3 2 ζ 0 , λ being (d + 1)-fold degenerate. For the sake of clarity, it proves convenient to relabel these eigenvalues as

λ 1 = 0, λ 2 = 0, λ 3 = - 3 2 ζ 0 , (38)
where λ 2 is d-fold degenerate. The associated hydrodynamic modes, {y β } d+2 β=1 are

y 1 = δρ, y 2 = δw, y 3 = 2 3 δρ + δθ. (39)
The + sign in the last equation stems from the fact that an increased density leads to enhanced dissipation, and hence, a smaller temperature.

B. Enforcing consistency with the linearized Boltzmann-Fokker-Planck equation description

We now turn our attention to the problem of finding the linearized hydrodynamic equations for a homogeneous perturbation, directly from the Boltzmann-Fokker-Planck equation. Enforcing consistency with the macroscopic considerations of section IV A above, we will infer useful properties on the hydrodynamic part of the spectrum of Λ(c). We first introduce the scaled deviation of the distribution function

δχ(c, τ ) = v d H n H [f (v, t) -f H (v)]. ( 40
)
The evolution of the scaled distribution is governed by

∂ ∂τ δχ(c, τ ) = Λ(c)δχ(c, τ ), ( 41
)
where the operator Λ(c) is the linearized Boltzmann-Fokker-Planck operator defined in (24). Let us also introduce the scalar product

f (c)|g(c) ≡ dcχ -1 H (c)f * (c)g(c), ( 42
)
where f * is the complex conjugate of f . Interestingly, the hydrodynamic modes introduced in (39) can then be written as

y β = ξβ |δχ , β = 1, 2, 3, (43)
where

ξ1 (c) = χ H (c), ξ2 (c) = cχ H (c), ξ3 (c) = c 2 d - 1 6 χ H (c). ( 44
)
Taking the scalar product of the linearized Boltzmann-Fokker-Planck equation (41) with the functions ξβ , we obtain the linear equations (34) (in the hydrodynamic time scale, that is, if we wait long enough so that fast modes have vanished) only if the spectrum of Λ admits the three eigenvalues written in [START_REF] García De Soria | [END_REF], and the associated "hydrodynamic" eigenfunctions, {ξ β } d+2 β=1 , obey the orthogonality condition ξβ1 In Appendix A, it is shown that the null eigenvalue is (d + 1)-fold degenerate, and the corresponding eigenfunctions, ξ 1 and ξ 2 , are worked out. Moreover, as a consequence of particle and total momentum conservation in a collision, ξ1 and ξ2 are the corresponding left eigenfunctions. We were not able to demonstrate that λ 3 is an eigenvalue of Λ, but we have shown explicitly that

|ξ β2 = δ β1β2 , β 1 , β 2 = 1, 2, 3. (45)
ξ3 |ξ β = 0, for β = 1, 2. ( 46
)
In the following, we will assume that Λ actually admits this third eigenvalue, with an unknown eigenfunction ξ 3 . With the help of this assumption, we will see in the next section that it is possible to define a projector in the hydrodynamic subspace, which opens the way for evaluating the variance of the total energy fluctuations.

V. ENERGY FLUCTUATIONS

In this section, we study the fluctuations of the global energy for a system in the stationary state. As we are interested in global quantities, it is convenient to define a global correlation function φ

H φ H (c 1 , c 2 ) ≡ dr 12 g H (r 12 , c 1 , c 2 ). ( 47
)
The energy fluctuations can be written as a functional of this correlation function and the one-particle distribution function χ H , as [31,[START_REF] García De Soria | [END_REF] (δE

) 2 H = m 2 4 N v 4 H dcc 4 χ H (c) + dc 1 dc 2 c 2 1 c 2 2 φ H (c 1 , c 2 ) . (48)
In order to evaluate the integral over the correlation function φ H , we start from (23), integrating over the position variable. Assuming periodic boundary conditions, the spacial gradient terms disappear and we have the following equation for φ

H [Λ(c 1 ) + Λ(c 2 )] φ H (c 1 , c 2 ) = Γ(c 1 , c 2 ), (49)
where

Γ(c 1 , c 2 ) = -T0 (c 1 , c 2 )χ H (c 1 )χ H (c 2 ) + ξ 2 0 ∂ ∂c 1 • ∂ ∂c 2 χ H (c 1 )χ H (c 2 ). (50)
The solubility condition for equation ( 49) is that Γ does not have components in the subspace associated to the null eigenvalue. In our case, this subspace is generated by { ξ1 , ξ2 }. Due to the conservation of the number of particles and total momentum in a collision, and to the symmetry of the second term of Γ, we have

ξ1 (c 1 ) ξ2 (c 2 )|Γ(c 1 , c 2 ) = ξ2 (c 1 ) ξ1 (c 2 )|Γ(c 1 , c 2 ) = 0, (51) ξ1
(c 1 ) ξ1 (c 2 )|Γ(c 1 , c 2 ) = ξ2,i (c 1 ) ξ2,j (c 2 )|Γ(c 1 , c 2 ) = 0, (52)
for i = j. The case i = j is analyzed in the Appendix B, where it is shown that ξ2,i (c 1 ) ξ2,i (c 2 )|Γ(c 1 , c 2 ) = 0. In order to prove this property, the presence of the second term of the right hand side of Eq. ( 50) is essential. Hence, the solubility condition holds and the problem of finding φ H with equation ( 49) is well defined. Let us also define a projector P 12 in the hydrodynamic subspace as

P 12 h(c 1 , c 2 ) = 3 β1=1 3 β2=1 ξβ1 (c 1 ) ξβ2 (c 2 )|h(c 1 , c 2 ) ξ β1 (c 1 )ξ β2 (c 2 ), (53)
where {ξ β } 3 β=1 are the right hydrodynamic eigenfunctions of the linearized Boltzmann-Fokker-Planck operator, { ξβ } 3 β=1 the orthogonal set introduced in the previous section, equation (44), and we have generalized the scalar product by

f (c 1 , c 2 )|g(c 1 , c 2 ) = dc 1 dc 2 χ -1 H (c 1 )χ -1 H (c 2 )f * (c 1 , c 2 )g(c 1 , c 2 ). (54)
Note that P 12 is a projector even if the function ξ3 is not the true left eigenfunction of Λ(c) (remember that ξ1 and ξ2 are the actual left eigenfunctions associated to the null eigenvalue). The fact that the set of functions { ξβ } 45) is enough to guarantee that P 2 12 = P 12 . Using this projector, we define the "hydrodynamic part" of φ H to be the function

φ (h) H (c 1 , c 2 ) ≡ P 12 φ H (c 1 , c 2 ) = 3 β1=1 3 β2=1 a β1β2 ξ β1 (c 1 )ξ β2 (c 2 ). (55)
The coefficients a β1β2 are the quantities we need to evaluate. As they are essentially the first moments of the correlation function φ H , they are directly related to the integral we have to calculate in (48). It is tempting to treat ξ 3 as if it was an actual left eigenfunction of the linearized Boltzmann-Fokker-Planck operator, and we will in the following use the approximation

P 12 Λ(c i ) = P 12 Λ(c i )P 12 , (56)
which allows to find a closed equation for φ

H . This approximation has already been invoked in other systems such as the freely evolving granular gas [31], or the probabilistic ballistic annihilation model [START_REF] Maynar | [END_REF]. With the information available on the linearized Boltzmann-Fokker-Planck operator, it seems the best that can be done technically. Let us also remark that the approximation is exact in the elastic limit. Then, applying the projector P 12 to the equation ( 49) and taking into account the approximation (56), we obtain the following expressions for the coefficients a β1β2

a β1β2 = - ξβ1 (c 1 ) ξβ2 (c 2 )| T0 (c 1 , c 2 )χ H (c 1 )χ H (c 2 ) λ β1 + λ β2 , (57)
where it has been assumed that λ β1 + λ β2 = 0. The coefficients associated to the vanishing eigenvalue cannot be calculated by equation ( 49), but are fixed by the boundary conditions. In Appendix C, the coefficients a β1β2 are evaluated. The expression for φ

(h)
H is finally given by φ (h)

H = a 11 ξ 1 (c 1 )ξ 1 (c 2 ) + a 13 [ξ 1 (c 1 )ξ 3 (c 2 ) + ξ 3 (c 1 )ξ 1 (c 2 )] + a 33 ξ 3 (c 1 )ξ 3 (c 2 ), (58)
where a 11 = -1, a 13 = -1 3 and a 33 can be obtained as a functional of the one-particle distribution function

a 33 = ξ3 (c 1 ) ξ3 (c 2 )| T0 (c 1 , c 2 )χ H (c 1 )χ H (c 2 ) 3ζ 0 . ( 59
)
An approximate expression is derived in appendix C, and reads

Taking into account (48) and (58), the variance of the energy fluctuations can finally be written as

σ 2 E = N (δE) 2 H E 2 H = (a 2 + 1) d(d + 2) 4 + d 2 a 33 -d 2 5 36 . (61)
In reference [29], the value of σ 2 E has been measured by means of the Direct Monte Carlo simulation method (DSMC). In Fig. 1, we compare our theoretical prediction (solid line given by ( 61)) with the DMSC simulations results (symbols). The agreement is satisfactory for the whole range of inelasticities, at variance with the theoretical attempt put forward in [29], which neglected velocity correlations. In particular, we note the non-trivial result in the elastic limit lim α→1 -σ 2 E (α) = d/3 (i.e. 2/3 on Fig. 1, which is well obeyed), while in the free cooling regime, this quantity vanishes [31]. We emphasize that the elastic limit is singular: the behaviour for elastic systems with α = 1 is not approached by taking the quasi-elastic limit α → 1 -(we note that the divergence of the different moments of the velocity distribution as α → 1 is nevertheless indicative of the absence of a stationary state when α = 1). Such a singularity has already been reported in 1 dimension [39], but, to our knowledge, not for two dimensional granular systems. It is also interesting to note that the singular nature of the quasi-elastic limit appears at 2 body level through the energy fluctuations, while, as far as rescaled distribution functions are considered, the 1 body level of description is regular, with a well behaved velocity distribution approaching a Gaussian form [22].

VI. CONCLUSIONS

The problem of the fluctuations of the total energy of a granular (inelastic) gas maintained in a non-equilibrium stationary state by a random acceleration has been addressed. A numerical study of this quantity had been performed by means of Monte Carlo simulations and an argument assuming uncorrelated non-Gaussian individual distribution function had been proposed in [29], without success. The main goal of this work was therefore to take due account of velocity correlations in order to study these fluctuations.

To this end, the standard description at the single particle level is not sufficient, and the two-particle correlation function is needed. We have derived the evolution equation for such an object, and particularized the analysis to the homogeneous stationary state that is reached by the system in the long time limit. Our work shows that this equation is not a straightforward generalization of its counterpart arising in the context of the undriven granular gas (i.e. by only changing the linearized Boltzmann-Fokker-Planck operator into its driven form). A non trivial non diagonal term appears in the Fokker-Planck equation for the two particle distribution function (contribution ∂ v1 • ∂ v2 in Eq. ( 9)), as a consequence of the coupling between velocities due to momentum conservation.

We have seen that, for our purposes, the exact knowledge of the hydrodynamic eigenfunctions is not needed. The important point is that we can construct a set of functions { ξβ } 3 β=1 , which are linear combinations of 1, c and c 2 , that are orthogonal to the right eigenfunctions {ξ β } 3 β=1 of the linearized Boltzmann-Fokker-Planck operator Λ. This orthogonality property holds for the "real" left eigenfunctions, that in our case correspond to the null eigenvalue (i.e. density and velocity fields associated to conserved quantities). The function ξ3 is not a left eigenfunction of Λ but it can be proved to be orthogonal to ξ 1 and ξ 2 . In a subsequent step, the linear hydrodynamic equations around the reference state are derived and from that knowledge, the hydrodynamic eigenvalues are identified and the variance of energy fluctuations subsequently derived.

Finally, our prediction has been successfully tested against the numerical results obtained by the Direct Monte Carlo simulation method for all the range of values of the coefficient of normal restitution α. This provides strong support for the theory developed here and assesses in retrospect the validity of our assumptions.

  GRANULAR GAS : TWO BODY KINETIC DESCRIPTION

a 33 =

 33 -15 + 7d + 14d 2 -3(-9 + d(9 + 2d))α + 30(1 + d)α 2 -6(9 + d)α 3 ) 9d(-19 + 2d(-7 + 3α) + 3α(9 + 2(-1 + α)α)).

Figure 1 :

 1 Figure 1: Scaled second moment of the energy fluctuations σ 2E as a function of the restitution coefficient α. The solid line is the theoretical prediction and symbols are the two dimensional Monte Carlo simulation results of Ref.[29].

dc 1 3 -

 13 dc 2 χ H (c 1 )χ H (c 2 )ϑ(c 1 , c 2 ), 3α)c 12 (C • c 12 ) 2 , (C12)and C = (c 1 + c 2 )/2. The coefficient b(α) can be evaluated using the expression of χ H (c) in the first Sonine approximation, Eq. (17), which yieldsb(α) = (1 + d)(3 + d)(2d(a 2 + 16(-1 + α) + 15a 2 α) + 16(-1 + α)(-1 + 2α 2 ) + a 2 (7 + 3α(-13 + 10(-1 + α)α)))Γ[(1 + d)/2] 2 1/2 π d+1/2 128d 2 (-2 + (5 + d)/2)(-1 + (5 + d)/2)Γ[d/2]Γ[-2 + (5 + d)/2](C13) If we take into account the explicit form of ζ 0 and a 2 , given in Eqs. (20) and (35) respectively, we obtain after some algebraa 33 = -15 + 7d + 14d2 -3(-9 + d(9 + 2d))α + 30(1 + d)α 2 -6(9 + d)α 3 ) 9d(-19 + 2d(-7 + 3α) + 3α(9 + 2(-1 + α)α)) .

Acknowledgments

It is a pleasure to dedicate this work to Jean-Jacques Weis. We thank Paolo Visco for useful discussions, and for providing us with the Monte Carlo data of Fig. 1. We would like to thank the Agence Nationale de la Recherche for

Appendix A: Eigenvalue problem for Λ

We consider here the eigenvalue problem for the homogeneous linear Boltzmann-Fokker-Planck operator Λ, defined in (24) Λ(c)ξ β (c) = λ β ξ β (c).

(A1)

We are interested in the eigenfunctions and eigenvalues associated with linear hydrodynamics and, to perform the analysis, similar techniques as in [31,[START_REF] Brey | Granular Gas Dynamics[END_REF][START_REF] García De Soria | [END_REF] will be required. Consider first the function

When the linearized operator Λ acts on χ H , we have

Taking into account the equation for χ H , Ec. ( 14), we obtain the following relation

Now, let us considerer the function

Taking derivate in the equation obeyed by χ H (c -w) with respect to w, and subsequently evaluating the result for w = 0, we obtain

Finally, we will consider the function

From the equation obeyed by ψ 3 (λc 1 ), we can take derivate with respect to λ, and evaluate the result for λ = 1. We arrive at an equation for ψ 3 (c 1 ),

From equations (A4), (A6) and (A8), we can identify two eigenfunctions of Λ. Making use of (A4) and (A8), it appears that

Hence, from Eqs. (A6) and (A9) we can conclude that the null eigenvalue is (d + 1)-fold degenerate with the eigenfunctions

Appendix B: Evaluation of the coefficient a2,i2,i

In this appendix, we show ξ2,i (c 1 ) ξ2,i (c 2 )|Γ(c 1 , c 2 ) = 0. The integral corresponding to the second term of Γ is simply

The other term can be written as

which is the desired result.

Appendix C: Evaluation of the coefficients a ββ ′

In this appendix we evaluate the coefficients a ββ ′ . As the number of particles and the total momentum are conserved quantities in our system, we have

Enforcing the above constraints, we obtain

As a consequence, the values of some coefficients follow

Of course, the coefficients associated to λ β1 + λ β1 = 0 could also have been calculated directly by equation (57), obtaining the same results. The coefficient a 33 is evaluated using (57) and it can be written in terms of the oneparticle distribution function as


Dario Marrocchelli

Mathieu Salanne

Paul A Madden

Christian Simon

P Turq

The construction

Keywords: GeO2, molecular dynamics, infrared spectrum, polarization effects, SiO2

come

In the vitreous and liquid states, at ambient pressure, germania (GeO 2 ) (a close structural analog of silica (SiO 2 )) forms a tetrahedrally coordinated three-dimensional network [START_REF] Micoulaut | The structure of amorphous, crystalline and liquid GeO 2[END_REF] .

Because of its lower abundance, its usage in practical applications is much less widespread than for the SiO 2 . Still, GeO 2 is used in several fields, mainly related to optical technologies.

For example, a mixture of SiO 2 and GeO 2 allows precise control of refractive index in optical fibres and waveguides. It is of interest to develop simulation methods to allow a detailed examination of the local structure in such mixtures and to predict the infrared spectrum, since this determines the long wavelength limit for their use as optical fibres.

The essential similarity between the structures of glassy GeO 2 , SiO 2 and BeF 2 was made clear by the study of vibrational properties by Galeener and co-workers 2 who showed that the inelastic neutron, infrared and Raman spectra of the different materials were closely related. More recent infrared studies of GeO 2 have also been reported [START_REF] Kamitsos | Raman and infrared structural investigation of xRb 2 O.(1x)GeO 2 glasses[END_REF][START_REF] Teredesai | Infrared spectroscopy of germanium dioxide (GeO 2 ) glass at high pressure[END_REF] . Recently a full set of partial structure factors were determined in this system by using the method of isotopic substitution in neutron diffraction experiments [START_REF] Salmon | Glass fragility and atomic ordering on the intermediate and extended range[END_REF][START_REF] Salmon | Structure of glassy GeO 2[END_REF] . The results show that the tetrahedral network structure is based on corner sharing Ge(O 1/2 ) 4 tetrahedra with a Ge-O average distance of 1.73 Å and with a mean inter-tetrahedral Ge-O-Ge angle of 132 • . They also show that the topological and chemical ordering in the network display two characteristic length scales at distances greater than the nearest neighbour.

Several molecular dynamics (MD) simulations have been undertaken to study the structural and vibrational properties of the disordered phases of GeO 2 . Classical MD simulations were performed on glassy and liquid GeO 2 [START_REF] Micoulaut | Simulated structure and thermal properties of glassy and liquid germania[END_REF][START_REF] Peralta | Structural and vibrational properties of amorphous GeO 2 : a molecular dynamics study[END_REF] , using pairwise additive potentials with partial charges developed by Oeffner and Elliott (OE) for modeling the α-quartz and rutilelike phases of GeO 2 [START_REF] Oeffner | Interatomic potential for germanium dioxide empirically fitted to an ab initio energy surface[END_REF] . More recently, first-principles molecular dynamics (FPMD) of glassy GeO 2 have also been reported [START_REF] Giacomazzi | Medium-range structural properties of vitreous germania obtained through first-principles analysis of vibrational spectra[END_REF][START_REF] Giacomazzi | Vibrational spectra of vitreous germania from first-principles[END_REF] . In principle, the amount of empirical information needed to set up a first-principles calculation is minimal and it would normally be the method of choice to study the physico-chemical properties of condensed phase systems. However, it is computationally very expensive compared to classical molecular dynamics, which is a major drawback when dealing with glassy systems. Classical and FPMD simulations were also compared one with each other in a combined study 12 for temperatures above 2530 K. In this study apparent inconsistencies between the properties predicted by the classical simulations an original one, which was fitted from an ab initio energy surface, and a so-called rescaled one, which was developed from the previous one in order to give a better reproduction of the vibrational properties. For example, the partial charge of the germanium ion was shifted from 1.5 e in the original potential to 0.94174 e in the rescaled potential. The apparent inconsistencies arise because different classical potentials were being used. Our objective is a single potential which will give a good description of structure and dynamics and allow long simulations on systems of many atoms.

The structure of the simulated system can directly be compared with the neutron / Xray diffraction experiments. Both classical and FPMD simulation provide structure factors which are in rather good agreement with experiments. For the classical molecular dynamics, the most noticeable difference is a shift to high q-values of the first sharp diffraction peak (FSDP), which was observed with both versions of the OE potential [START_REF] Micoulaut | Simulated structure and thermal properties of glassy and liquid germania[END_REF][START_REF] Peralta | Structural and vibrational properties of amorphous GeO 2 : a molecular dynamics study[END_REF] . As this peak is associated to medium range distances [START_REF] Salmon | Glass fragility and atomic ordering on the intermediate and extended range[END_REF] , this means that the Ge-O-Ge bond angle is too large and the topology of the tetrahedral network is not correctly described. The FPMD simulations provide a diffraction pattern with improved agreement [START_REF] Giacomazzi | Vibrational spectra of vitreous germania from first-principles[END_REF] , but the small size / timescales explored lead to large statistical uncertainties, mainly in the region of this FSDP.

The vibrational properties can also be straightforwardly calculated. The vibrational density of states (VDOS) was estimated from FPMD and classical MD using the rescaled OE potential, with an overall reasonable agreement with the available experimental data extracted from inelastic neutron scattering experiments. The infrared absorption spectrum, which originates from the polarization fluctuations associated with motion of the atomic charges, has only been reported from the FPMD simulation. The agreement with experiment 2 was overall good, despite a small shift of the peaks towards lower frequencies.

The dynamical properties which involve structural rearrangements are out of reach of the FPMD. The only studies on GeO 2 employ classical MD, and were aimed at determining the diffusion coefficients in the liquid state [START_REF] Micoulaut | Simulated structure and thermal properties of glassy and liquid germania[END_REF]12 . Hawlitzky et al. 12 observed important differences between their values and those obtained by Micoulaut et al. [START_REF] Micoulaut | Simulated structure and thermal properties of glassy and liquid germania[END_REF] ; their hypothesis for rationalizing this discrepancy was inadequate equilibration in the latter. In fact, as noted above, it is due to the use of the two different OE potentials in these studies. Indeed, in simulations performed with the rescaled OE potential the fluidity of the system is enhanced by more than one order of magnitude compared to the results obtained with the original In this work we describe the parameterization of a new classical MD interaction potential for GeO 2 , by using some techniques we have developed for halides and oxides [START_REF] Madden | From first-principles to material properties[END_REF][START_REF] Jahn | Modeling Earth materials from crustal to lower mantle conditions: a transferable set of interaction potentials for the CMAS system[END_REF] , which is entirely based on first-principles electronic structure calculations. In order to provide an accurate, transferable description of the interactions, these potentials include dipole polarization effects and the ions carry full valence charges. The various structural, dynamical and vibrational properties are then compared with the corresponding experimental results.

In future work we will use such potentials to study mixtures and the effect of pressure on glassy GeO 2

II. INTERACTION POTENTIAL DEVELOPMENT

The interaction model used in this work is related to that used in previous work [START_REF] Castiglione | Polarization effects in the simulation of lead (II) fluoride[END_REF][START_REF] Castiglione | Fluoride ion disorder and clustering in superionic pbf 2[END_REF] , with a couple of minor differences. The model (known as DIPPIM) includes a pair potential, together with an account of the polarization effects which result from the induction of dipoles on the O 2-ions. The parameters for these potentials were obtained by force-matching them to first-principles reference data [START_REF] Madden | From first-principles to material properties[END_REF] . Such an approach was successfully applied in the case of other oxide materials [START_REF] Jahn | Modeling Earth materials from crustal to lower mantle conditions: a transferable set of interaction potentials for the CMAS system[END_REF] . Tangney and Scandolo 17 also used a similar interaction model to study silica, in which all the parameters including partial charges on the ions were fitted against first-principles data. Here we use formal ionic charges (Ge 4+ , O 2-) which should ensure a better transferability. In the next paragraphs we will give a brief description of the model used, the first-principles reference calculations and the force-matching procedure.

A. The DIPole-Polarizable Ion Model (DIPPIM)

The interatomic potential is constructed from four components: charge-charge, dispersion, overlap repulsion and polarization. The first three components are purely pairwise additive:

V qq = i≤j q i q j r ij (1)
4 where q i is the formal charge on ion i. The dispersion interactions include dipole-dipole and dipole-quadrupole terms

V disp = - i≤j [ f ij 6 (r ij )C ij 6 r 6 ij + f ij 8 (r ij )C ij 8 r 8 ij ]. (2)
Here C ij 6 and C ij 8 are the dipole-dipole and dipole-quadrupole dispersion coefficients, respectively, and the f ij n are the Tang-Tonnies dispersion damping function, which describe short-range corrections to the asymptotic dispersion term. The short range repulsive terms are approximately exponential in the region of physical interionic separations. The full expression used here for the short range repulsion is:

V rep = i≤j A ij e -a ij r ij r ij + i≤j B ij e -b ij r 2 ij , (3)
where the second term is a Gaussian which acts as a steep repulsive wall and accounts for the anion hard core; these extra terms are used in cases where the ions are strongly polarized to avoid instability problems at very small anion-cation separations [START_REF] Castiglione | Polarization effects in the simulation of lead (II) fluoride[END_REF] . The polarization part of the potential incorporates dipolar effects only. This reads:

V pol = i,j -q i µ j,α f ij 4 (r ij ) -q j µ i,α f ji 4 (r ij ) T (1) α (r ij ) - i,j µ i,α µ j,β T (2)
αβ (r ij ) + i 1 2α i | µ i | 2 . ( 4
)
Here α i is the polarizability of ion i, µ i are the dipoles and T [START_REF] Micoulaut | The structure of amorphous, crystalline and liquid GeO 2[END_REF] , T [START_REF] Galeener | Comparison of the neutron, raman, and infrared vibrational spectra of vitreous sio 2 , geo 2 and bef 2[END_REF] are the charge-dipole and dipole-dipole interaction tensors:

T (1) α (r) = -r α /r 3 T (2) αβ (r) = (3r α r β -r 2 δ αβ )/r 5 . (5)
The instantaneous values of these moments are obtained by minimization of this expression with respect to the dipoles of all ions at each MD timestep. This ensures that we regain the condition that the dipole induced by an electrical field E is αE and that the dipole values are mutually consistent. The short-range induction effects on the dipoles are taken into account by the Tang-Toennies damping functions: The parameters d ij determine the range at which the overlap of the charge densities affects the induced dipoles, the parameters c ij determine the strength of the ion response to this effect. It is important to notice that anion-anion damping terms have been taken into account, contrary to what was done in [START_REF] Castiglione | Polarization effects in the simulation of lead (II) fluoride[END_REF] . The addition of anion-anion damping terms was found to greatly improve the ability to match the first-principles data.

f ij n (r ij ) = 1 -c ij e -d ij r ij n k=0 (d ij r ij ) k k! . ( 6

B. The DFT reference calculations

The parameters in the interaction potential are determined by matching the dipoles and forces on the ions calculated from first-principles on condensed phase ionic configurations [START_REF] Madden | From first-principles to material properties[END_REF][START_REF] Jahn | Modeling Earth materials from crustal to lower mantle conditions: a transferable set of interaction potentials for the CMAS system[END_REF] . Starting from the empirical pair potential in [START_REF] Micoulaut | Simulated structure and thermal properties of glassy and liquid germania[END_REF] , we obtained atomic configurations for GeO 2 by running short MD simulations on small cells (150 ions); a total of two liquid configurations were obtained for GeO 2 ; for each of these, the Hellman-Feynman forces acting on individual ions of the simulation cell were calculated using planewave-DFT code CPMD [START_REF]The CPMD consortium[END_REF] . In all the calculations we used norm-conserving pseudopotentials and planewave energy cut-offs of 1360 eV ; all calculations were performed using the generalized gradient approximation (GGA) according to the Perdew, Burke and Ernzernhof (PBE) scheme. For the calculation of first-principles dipoles, the Kohn-Sham orbitals are localized via a Wannier transformation to construct maximally localized Wannier functions (MLWF) and the dipoles determined from the positions of the centres of the Wannier functions associated with each ion [START_REF] Aguado | Multipoles and interaction potentials in ionic materials from planewave-dft calculations[END_REF] .

C. Fitting procedure

The potential parameters are optimized by fitting the forces and dipoles, obtained with the DIPPIM potentials for the reference configurations, to the respective results from the DFT calculations; the 2 configurations provide a total of about 1800 data points to fit, comprising three Cartesian force components of each individual ion and three components for the dipoles. While most of the potential parameters are left as free parameters in the fits, there are some exceptions. The O melts from first-principles calculations 20 . One problem in DFT calculations is the uncontrolled representation of the dispersion interaction. Although dispersion only contributes a tiny fraction to the total energy, it has a considerable influence on phase transition pressures and on the material density and stress tensor. For this reason, we decided not to include the dispersion parameters in the fits as discussed in [START_REF] Madden | From first-principles to material properties[END_REF] and to add them afterwards. We used the parameters in 21 for the oxygen terms and rescaled these by the Ge polarizability for the anion-cation terms. The Gaussian parameters too were added after the fit and then the comparison with the first-principles data was run again to check that its quality remained unchanged.

In Figure 1 we report the agreement between the DFT calculated forces and dipoles and those predicted by the fitted potential, for a set of representative ions; the abscissa shows the x component of the forces and dipoles, while the ordinate is simply the ion number.

The quality of the representation is quite good and comparable with the one obtained in [START_REF] Wilson | The construction and application of a fully flexible computer simulation model for lithium oxide[END_REF] ;

this can be regarded as a very good result, considering the simplicity of the model (both ion shape deformation effects and quadrupoles are not taken into account in the present model [START_REF] Madden | From first-principles to material properties[END_REF][START_REF] Jahn | Modeling Earth materials from crustal to lower mantle conditions: a transferable set of interaction potentials for the CMAS system[END_REF] ). The parameters obtained for the DIPPIM potential are summarized in Table I. We performed MD simulations in the NV T ensemble through the use of the Nosé-Hoover chain thermostat method [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Martyna | Nosé-hoover chains: the canonical ensemble via continuous dynamics[END_REF] . The simulation cell contained 288 oxide and 144 germanium ions, and its volume was set to 6856 Å3 in order to match the experimental density at 300 K, which is 3.66 g.cm -3 . We used a time step of 1 fs to integrate the equations of motion and the minimization of the polarization energy was carried out with a conjugate gradient method.

The system was studied in both the liquid and amorphous phases.

For the amorphous phase, the system was first equilibrated for 2 ns at the temperature of 4500 K. We have then cooled down the system with a cooling rate of 3.33×10 11 K.s -1 , by rescaling the velocities and decreasing the target temperature of the thermostat in order to reduce the temperature by 50 K every 150000 MD steps. Finally, a 5 ns long simulation was undertaken at 300 K, from which we computed all the data presented in this article. This procedure was performed with two different potentials, the DIPPIM, and the original OE, which has extensively been used in the literature [START_REF] Oeffner | Interatomic potential for germanium dioxide empirically fitted to an ab initio energy surface[END_REF] . A simulation at 300 K was also performed with the rescaled OE potential (with one of the configurations obtained from simulations with the original OE potential as a starting configuration) for the study of the vibrational properties.

For the liquid, a series of long simulations with the DIPPIM potential at the temperatures of 3600 K, 3800 K, 4000 K, 4200 K, 4400 K, 4500 K, 4600 K, 4800 K and 5000 K was undertaken. In all the cases, an equilibration run was performed so that the slower species, Ge, moved on average, at least a distance of 5.5 Å; a subsequent run of the same length was made to accumulate enough statistics for the mean-squared displacements curves. At these temperatures, the runs were between 100,000 and 1,000,000 steps. Both the original and rescaled OE potentials were also used within the same simulation conditions at a temperature of 3600 K.

IV. STATIC STRUCTURE

Neutron diffraction experiments with the isotope substitution technique allowed Salmon to determine very precisely the full set of partial radial distribution functions (RDF) [START_REF] Salmon | Glass fragility and atomic ordering on the intermediate and extended range[END_REF] . These provide a good check of the validity of the potential.

The experimental RDFs are given in figure 2, together with the ones obtained from our simulations at 300 K involving the DIPPIM and the original OE potential; it can be readily observed that our DIPPIM potential gives a closer match with the experimental data than the OE potential, for which all the characteristic first-neighbour distances are slightly overestimated. Concerning the shape of the function, the peaks obtained for the Ge-Ge and Ge-O partials are too sharp but it has be be remembered that other factors, of the structure might not be sufficient for it to be capable of predicting the dynamical and vibrational properties.

V. DYNAMICS FIG. 3: Time dependence of the mean squared displacement curves for both the oxygen and germanium ions. The temperatures ranges between 3600 K and 5000 K.

In this section we will calculate diffusion coefficients for the liquid and compare them both with those obtained with the OE potentials and with the experimental values. The diffusion coefficient of species α can be obtained from the slope of the mean squared displacement at long times, i.e.

D α = lim t→∞ 1 6t r 2 α (t) , (7)
where

r 2 α (t) = 1 N α Nα i=1 |r i (t) -r i (0) 2 > . ( 8
)
The mean squared displacement curves were evaluated from a series of long simulations at temperatures between 3600 K and 5000 K; these are shown on a log-log scale on figure 3. No plateau is observed even at the lowest temperature, which shows that the system remains liquid in the range of the study.

The diffusion coefficients were extracted using equation 7; these are reported in figure 4 and compared with those obtained with the OE and rescaled OE potentials. repulsion were lowered by a factor of approximately 2.5; as a consequence, these ions can diffuse much more easily in a softer environment.

VI. INFRARED SPECTRUM

The infrared spectra of ionic melts with polyvalent cations exhibit discrete bands attributable to the vibrational motion of the local coordination complexes around the cations.

They originate from the polarization fluctuations associated with motion of the ionic charges.

The inclusion of polarization effects for the oxide ions in our model may influence the predicted spectrum in two ways 26 . First, the interactions of the oxide ion dipoles may alter the local structure of the network and the strength of the bonds, which may introduce a shift of the vibrational frequencies. Second, the induced dipoles will themselves be responsible for absorption, as they too contribute to the total polarization fluctuations. The absorption coefficient in the presence of these extra moments is calculated from the imaginary part of the total dielectric function [n(ν)α(ν) = 2πν ( (ν))], which can be determined following Caillol, Levesque and Weis 27,28 as where

(ν) -∞ = β 3 0 V M(0) 2 + 2πıν M • M ν (9) +2 M • J ν + ı 2πν J • J ν
J • J ν = ∞ 0 e 2πıνt J(t) • J(0) dt, (10)
J(t) is the charge current J(t) = N i=1 q i v i (t) and M(t) is the total system induced dipole moment, M(t) = N i=1 µ i (t). Figure 5 shows the imaginary part of the total dielectric function calculated for simulations carried out with the DIPPIM potential, compared to the one obtained from FPMD simulations [START_REF] Giacomazzi | Medium-range structural properties of vitreous germania obtained through first-principles analysis of vibrational spectra[END_REF][START_REF] Giacomazzi | Vibrational spectra of vitreous germania from first-principles[END_REF] and from several experiments 2,3 (the spectrum proposed by Kamitsos et al.

shows the absorption coefficient, here we transformed it into the dielectric function to facilitate comparisons). Three main bands are observed on the spectra, in good agreement with the other studies. The corresponding characteristic peak frequencies for the absorption spectra are summarized in table II. The agreement is very satisfactory; the frequencies obtained match very closely to the various experimental results (even though these show a significant scatter), while the FPMD simulation seems to systematically underestimate them. The relative intensities of the bands is also in good agreement with experiment, though the highest frequency band at around 894 cm -1 is not intense enough in our spectrum.

The importance of polarization effects in determining the relative intensities of the bands can be demonstrated by separating the various contributions to the absorption spectrum.

In the case of glassy silica 26 and beryllium fluoride [START_REF] Heaton | A first-principles description of liquid BeF 2 and its mixtures with LiF: 1. potential development and pure BeF 2[END_REF] , it was observed that the interference between the induced dipoles and permanent charge contributions to the total polarization, contained in the M • J ν cross term, is responsible for the changes of the relative intensities. at the same intensities. Note that the cross-term strongly reduces the intensity of the two low frequency bands relative to that which would be obtained from the charge fluctuations alone, whereas for the high-frequency band the cancellation is much weaker. For the high frequency band the relationship of the charge-charge and cross-terms is different to that found previously for SiO 2 and BeF 2 26,29 . In the latter cases, the charge-charge and crossterms had the same sign, so that the net band intensity was slightly larger than that predicted by the charge-charge term alone. This difference in the behaviour of the calculated spectra might arise from the inclusion of an anion-anion damping term in the polarization part of the interaction potential here, which was not the case for the previous studies.

FIG. 6: Decomposition of the total IR absorption spectrum in its three component parts.

We can now attempt to associate these IR frequencies with the vibrational modes of the system. The upper part of figure 7 shows the vibrational density of states (VDOS) of each ion. Some of these features can be associated to the vibrations of the GeO of each GeO 4 unit. For example, for the symmetric stretching motion

v A i = iα=1-4 v iα (11)
where i labels the Ge of a tetrahedral complex, and iα the four oxide anions at the vertices.

v iα is the projection of the relative velocity of iα along the i → iα bond, i.e.

v iα = (v iα -v i ) • (r iα -r i ). (12)
Similar expressions may be written down for the velocities of the other symmetry coordinates [START_REF] Pavlatou | The interpretation of vibrational spectra of ionic melts[END_REF] . The corresponding DOS is then obtained by a Fourier transform of the corresponding velocity autocorrelation function:

DOS A (ν) = ∞ 0 e 2πıνt v A i (t) • v A i (0) dt. ( 13
)
Were the network to vibrate as a collection of isolated tetrahedra, the spectra of each tetrahedral-VDOS would exhibit a single peak corresponding to the characteristic frequency of the corresponding normal mode of vibration. The four tetrahedral-VDOS are shown on the lower part of figure 7, it is clear that it is not the case here. Each spectra consists in the superposition of several bands, which shows the existence of some important coupling between the symmetry coordinates. It is therefore difficult to assign the IR features to the different modes. In particular, all of the three bands observed in the IR spectrum appear in the tetrahedral-VDOS of the bending E and F 2 modes. The other F 2 modes, which is associated to some kind of stretching of the tetrahedra, only exhibits the medium and high frequency modes. Finally, the A 1 stretching spectrum exhibits two bands at frequencies centered on the values of 510 and 970 cm -1 , which do not appear on the IR spectrum. Such a mode of this symmetry is indeed not expected to be IR active for an isolated tetrahedral molecule, since it corresponds to the symmetric breathing mode. Coming back to the total VDOS of each ion, the low frequency vibrations are not associated to any of the tetrahedral modes. In fact these are due to some slow collective bending of two tetrahedra linked by a bridging oxide anion [START_REF] Trachenko | Floppy modes in silica glass under pressure[END_REF] , so that they would not appear in the DOS of a single tetrahedron.

The vibrational properties of the system as obtained with the original and rescaled OE interaction potentials have also been computed. As these potentials do not include any polarization effect, it is only possible to compare the charge-charge current contributions to the IR absorption spectrum with the one obtained with the DIPPIM potential. The obtained spectra are given on figure 8. Unsurprisingly, the original OE potential seem to fail completely in describing the vibrational properties of GeO 2 ; this is the reason why the rescaled OE potential had to be developed. For that potential, the agreement is better, in particular for both the small and high frequency bands. Concerning the medium frequency band, observed at around 604 cm -1 for the DIPPIM potential, it does not appear for the rescaled OE potential. Instead there appears a shoulder to the small frequency band, for frequencies ranging from 400 to 500 cm -1 . As these bands were associated to the F 2 stretching 17 with the three different potentials (the two OE curves were rescaled with a factor q 2 O /4, where q O is the partial charge carried by the oxide ion in the corresponding potential.

VII. CONCLUSION

In conclusion, the obtained potential seem to reproduce all the studied properties (structural, dynamical and vibrational) to a high degree of precision. Since there was no reference at all to any experimental data in the parameterization of this potential, this represents a strong test of the model reliability. The study of pressurized germania is among our top priorities for future work. The potential is of a similar form to that introduced to describe silicates [START_REF] Jahn | Structure and dynamics in liquid alumina: simulations with an ab initio interaction potential[END_REF] ; since these potentials have proved to be transferable, we will be able, in the future, to focus on GeO 2 -SiO 2 mixtures as well. These mixtures, as already mentioned in the introduction, attract great interest as they are widely used in optical fibers and waveguides.

  are noted. In reality, two different parameter sets were proposed by OE:

) 5 Page

 5

FIG. 1 : 7 Page

 17 FIG. 1: Agreement between the DFT calculated dipoles (top panel) and forces (bottom panel) and those predicted by a fitted DIPPIM potential model for a set of 50 representative ions. In each panel the ordinate gives the ion number and the abscissa gives the DFT values of the x component of the forces and dipoles as points and the fitted values as the continuous line.

8 Page

 8

FIG. 2 : 10 Page 10 of 32 URL

 21032 FIG. 2: Radial distribution functions obtained from the 300 K simulations with the DIPPIM and original OE interaction potentials. These are compared to the experimental ones obtained from neutron diffraction experiments 5 .

Figure 4 shows 11 Page 11 of 32 URLFIG. 4 :Page 12 of 32 URL

 1132432 FIG.4: Diffusion coefficients obtained with a DIPPIM potential at different temperatures; data points obtained with both the original and rescaled OE potentials at 3600 K are also presented.

FIG. 5 :

 5 FIG. 5: Imaginary part of the dielectric function, calculated for simulations carried out with the DIPPIM potential, compared to FPMD simulation 10 and experimental results 2,3 .

13 Page 13 of 32 URL

 1332

Figure 6 14 Page

 614 Figure6shows the decomposition of the absorption spectrum in the case of glassy GeO 2 .The two main contributions are the charge fluctuation and the cross term; both show bands

16 Page 16 of 32 URLFIG. 7 :

 16327 FIG. 7:The vibrational densities of states, for the DIPPIM, of both type of ions (top panel), and their projections in terms of the A (symmetric stretch), F 2 (asymmetric stretch), E (bend), and F 2 (bend) modes of the Ge0 4 tetrahedra (bottom panel).

  the CMPC for his PhD funding. This work was carried out under the HPC-EUROPA++ project (project number: 211437), with the support of the European Community -Research Infrastructure Action of the FP7 "Coordination and support action"

Page 4 of 32 URL: http://mc.manuscriptcentral.com/tandf/tmph

	Molecular Physics
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

  2-polarizability, for instance, was fixed to α O 2-= 11 a.u., i.e.

	Molecular Physics
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

the value Salmon et al obtained experimentally in 5 . This value for the oxide polarizability is also compatible with the range of values obtained for this quantity in several magmatic 6 Page 6 of 32 URL: http://mc.manuscriptcentral.com/tandf/tmph

TABLE I :

 I Parameters in the repulsive, polarization and dispersive parts of the DIPPIM potential.All values are in atomic units. The oxide ion polarizability was set to 11 au.

		O-O O-Ge Ge-Ge
	A ij 17.78 40.89 2958
	a ij 3.893 1.186 5.96
	B ij 50000 50000	0
	b ij 1.1 1.75	/
	C ij 6	44	4	0
	C ij 8	853	50	0
	c ij 6	1.0	1.0	/
	c ij 8	1.0	1.0	/
	d ij 6	1.0	1.5	/
	d ij 8	1.0	1.5	/
	d ij pol 2.208 1.977	/
	c ij pol 1.770 1.709	/
	III. SIMULATION DETAILS

TABLE II :

 II Characteristic frequencies of the imaginary part of the dielectric function as obtained from the DIPPIM model, compared to previous FPMD simulations[START_REF] Giacomazzi | Medium-range structural properties of vitreous germania obtained through first-principles analysis of vibrational spectra[END_REF][START_REF] Giacomazzi | Vibrational spectra of vitreous germania from first-principles[END_REF] and to experimental results[START_REF] Galeener | Comparison of the neutron, raman, and infrared vibrational spectra of vitreous sio 2 , geo 2 and bef 2[END_REF][START_REF] Kamitsos | Raman and infrared structural investigation of xRb 2 O.(1x)GeO 2 glasses[END_REF][START_REF] Teredesai | Infrared spectroscopy of germanium dioxide (GeO 2 ) glass at high pressure[END_REF] . All the frequencies are given in cm -1 .

	DIPPIM 300 601 894
	FPMD 10 244 498 787
	Exp. 1 2 280 567 858
	Exp. 2 3 315 585 915
	Exp. 3 4	560 870

  4 tetrahedra by comparing them with the tetrahedral-VDOS as obtained from the correlation functions of the velocities of various local tetrahedral symmetry coordinates. The latter would correspond to the velocities of GeO 4 normal modes, if the tetrahedra were isolated (i.e. not mechanically linked in a network). For a tetrahedral molecule, vibrational normal coordinates ν 1 (A 1 , symmetric stretch), ν 2 (E, bend), ν 3 (F 2 , asymmetric stretch) and ν 4 (F 2 , bend) are expected[START_REF] Pavlatou | The interpretation of vibrational spectra of ionic melts[END_REF] . We may obtain velocities associated with the symmetry coordinates 15

	Page 15 of 32
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

URL: http://mc.manuscriptcentral.com/tandf/tmph

  and bending modes on the tetrahedral-VDOS, this means that this modes are not correctly depicted by the rescaled OE potential. FIG. 8: Comparison of the charge-charge current contributions to the infrared spectrum obtained

	Page 17 of 32	Molecular Physics
	F o
	r
	P
	e e r
		R e v i e w
		O n
		l y

URL: http://mc.manuscriptcentral.com/tandf/tmph

M. Hawlitzky, J. Horbach, S. Ispas, M. Krack, and K. Binder. Comparative classical and ab initio molecular dynamics study of molten and glassy germanium dioxide. J. Phys.:Condens. Matter, 20(28):285106, 2008.

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Acknowledgements

We thank Philip Salmon for providing his diffraction data. Calculations were carried out at the EaStCHEM computing facility (http://www.eastchem.ac.uk/rcf). DM thanks 18 Page 18 of 32 URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics