T Guérin

David S Dean

in heterogeneous periodic nonequilibrium systems

Keywords: numbers: 05.60.Cd, 02.50.Ey, 05.10.Gg, 05.40.a

I. INTRODUCTION

Fokker-Planck (FP) equations, and their associated stochastic differential equations (SDE) arise in a huge variety of physical systems [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF][START_REF]Stochastic differential equations[END_REF][START_REF] Gardiner | Handbook of stochastic methods for physics, chemistry and the natural sciences[END_REF]. In fluid mechanics and hydrology, such equations describe the transport of passive scalars advected by a fluid and locally dispersed by a microscopic molecular diffusivity. They also describe the diffusion of colloids and polymers in soft matter physics, as well as charge and energy transport in solid state physics. In many systems, damping due to the environment is very strong and inertial effects are thus negligible. In such systems, the evolution of the probability density function (pdf) p(x, t) of a tracer particle at position x and time t obeys the FP equation

∂ t p = ∂ xi ∂ xj [κ ij (x)p] -u i (x)p = -H x p, (1)
where the Einstein convention on indexes summation is used, and the indexes i, j run from 1 to the spatial dimension d. The time independent tensor κ ij (x) is the local diffusivity tensor and can, in general, vary in space in media of non-uniform composition, since diffusivity is depends on the local material properties of the system. The local diffusivity may also vary in space due to temperature inhomogeneities, or in the case of particle diffusion in viscous fluid flows close to solid boundaries [START_REF] Brenner | [END_REF]. The term u i (x) represents a drift term generated by forces such as buoyancy forces or electromagnetic forces acting on colloids, fluid flows, as well as thermodynamic effects, such as thermodiffusion in the presence of temperature gradients [5]. The presence of impenetrable obstacles can also be treated by introducing surfaces at which there is no normal flux,

n i (x){∂ xj [κ ij (x)p] -u i (x)p} = 0, x ∈ S, (2)
where n i represents the normal vector to the surface S of the obstacles. Equation [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF] describes the motion of tracer particles at small scales in a wide range of heterogeneous media. At a coarse-grained level, the motion is described by effective transport coefficients characterizing the average flux of particles and the relative spreading with time of initially close particles. From a theoretical point of view, the dispersion of tracer particles is defined in terms of two quantities, the mean displacement µ i (t) during a time t in the direction i

µ i (t) = X i (t) -X i (0) , (3)
(where X(t) denotes the position of a tracer particle and • • • denotes ensemble averaging), and the correlation of the dispersion in the directions i and j defined by σ ij (t) = [X i (t)-X i (0)][X j (t)-X j (0)] -µ i (t)µ j (t). ( 4)

The quantity σ ij (t) thus characterizes the dispersion of a cloud of particles about its mean position. Of particular interest is the large time behavior of the dispersion properties,

µ i (t) V i t, (t → ∞) (5) σ ij (t) 2 t (D ij + C ij /t), (t → ∞), (6)
where V i is the average particle drift, D ij is the late time effective diffusion tensor, and C ij gives the leading order corrections and provides information on how fast the system approaches the diffusive limit. We will see that the asymptotic form ( 6) is valid if the medium is spatially periodic. The quantities V i and D ij characterize the late time, large scale transport properties of the system and are routinely measured in experimental systems by single particle tracking methods [6] or by ensemble based measurements such as fluorescence recovery after photobleaching (FRAP) [7]. Effective transport coefficient are also important for estimating the spread of pollutants, filtration, mixing and chemical reaction times [8].

Transport and dispersion properties in non-equilibrium periodic systems have been investigated at length in various contexts. In a pioneering, classical work, Taylor [9] considered the dispersion of particles moving in viscous flows in cylindrical containers, and showed that the effective diffusion coefficient could be orders of magnitudes larger than the molecular diffusivity. Effective transport in incrompressible fluid flows were also investigated at length, such as in the case of diffusion in Rayleigh-Bénard convection cells [10][11][12], in frozen turbulent flows [13] or in porous media [14][15][16][17][18][19][20] using homogenization theory. However, these approaches have been restricted so far to incompressible flows and it is not clear how to generalize them. In another context, results were also derived, with the methods of statistical physics, for particles diffusing in periodic potentials when the diffusivity is constant [21][22][START_REF] Zwanzig | Proc. Natl. Acad. Sci. U S A[END_REF][START_REF] Gennes | [END_REF][25]; in here the potential barrier slows down the diffusion as the particle tends to become trapped in local mining. The non-equilibrium problem of transport in one-dimensional potential in the presence of external constant force (still with uniform microscopic diffusivity) was also considered [26][27][28][29]; this system is predicted to exhibit a huge increase of the effective diffusion near a critical tilting force. Results have also been obtained in the more general case where the noise amplitude is a periodic function of position [30][31][32], notably in the context of the modeling of transport in periodic channels within the Fick-Jacobs approximation (where the full Fokker-Planck equation is approximated by an effective one dimensional one). In many of the above cases, results only exist in one dimension, however their generalization to higher dimensions is implicit in the general results we will present here.

In this paper, we present a general formalism to compute the time-dependent dispersion tensor σ ij (t) (section III), the late time diffusion tensor D ij (section IV) and the corrections C ij describing the approach to the diffusive limit (section V). A summary of the results is given in section II for the reader who might not be interested in the technical aspects of the calculations. The results derived are valid for any system described by Eq. ( 1) for periodic κ ij (x) and u i (x); they encapsulate all known results for this problem and also generalize them to more complex situations. We will explicitly compare them to the expressions of the literature in various cases (section VI), and we will also give an explicit fully analytical solution for the one dimensional problem (section VII). A brief derivation of our main results has already appeared in Ref. [33], where they were used to analyze diffusion in a system with a spatially varying diffusivity in the presence of a uniform applied force. The present paper extends this recently developed approach, new features are the computation of the time-dependent dispersion properties (instead of only the late time diffusive limit) and that the incorporation of the effect of impenetrable obstacles with no-flux boundary conditions. Furthermore, the derivation presented here is different from that appearing in Ref. [33], as it is not based on the SDE formulation used there. Finally, the present paper contains explicit expressions for the first temporal corrections C ij to the late time diffusion tensor. We will show that the 1/t decay of the late time correction [Eq. (6)] occurs in all dimensions. This universality can be understood as, for periodic systems, the Fokker-Planck equation restricted to the basic cell Ω will in general have a gap and not a continuum (whose density would be dependent on the spatial dimension) of eigenvalues near zero.

II. NOTATIONS AND SUMMARY OF RESULTS

We assume that the space can be divided into individual cells Ω, and we that the diffusivity tensor κ ij (x) and the drift u i (x) are periodic in space, and thus satisfy

κ ij (x) = κ ij (x + k), u i (x) = u i (x + k), (7)
for all vectors k joining the centers of the cells composing the periodic structure. Denoting by X(t) the position of an individual particle at time t, we define the propagator p(x, t|y) of the stochastic process X(t) in infinite space; p(x, t|y) is thus the probability density to observe a particle at position x at time t starting from y at time 0. Mathematically, p is the solution of the equation in infinite space

(∂ t + H x )p(x, t|y) = 0 ; p(x, 0|y) = δ(x -y). (8)
We also introduce the process X(t) as the position of the particle X "modulo Ω", X is defined as the position of the particle in a reference cell Ω 0 obtained after an integer number of translations of X along the vectors joining the centers of the cells composing the periodic structure. We introduce the propagator P (x, t|y) for the process modulo Ω, P thus represents the probability density, starting at position y at time 0, to be at time t at position x modulo Ω. P is the solution of the FP equation ( 1) with periodic boundary conditions on the boundaries of Ω and initial value P (x, 0|y) = δ(x -y). There is an obvious relation between these propagators,

P (x, t|y) = k p(x + k, t|y), (9)
where k represents the vectors joining the center of the cell containing x to the centers of all the other cells of the periodic structure. The value of P at late times is denoted by P 0 and is the steady-state distribution of particles in one period (modulo translations along the lattice vectors). In contrast, the propagator in infinite space p does not tend to a steady state, but becomes Gaussian at long times, with a variance growing linearly as t, representing the spreading of the particle density. In what follows, we always assume that at t = 0 we are in the steady state, that is to say that the process X(t = 0) modulo Ω has the distribution P 0 . Before entering the details of the calculations, let us give briefly our main results. First, we will show that the temporal Laplace transform σij (s) = ∞ 0 σ ij (t)e -st dt of the dispersion tensor is given by σij

(s) = 2 s 2 Ω dx κ ij (x)P 0 (x) - 1 s 2 Ω
dxdy[ũ i (x)ũ * j (y) + ũj (x)ũ * i (y)] P (x, s|y)P 0 (y). (10) In this relatively compact formula, P (x, s|y) represents the Laplace transform of P (x, t|y) -P 0 (x). We have also defined an effective velocity field ũ that takes into account the presence of obstacles (if present)

ũi (x) = u i (x) -n j (x)κ ij (x)δ S (x), (11)
where δ S is a surface delta function; to be clear with the notations we precise that

Ω dxn j (x)κ ij (x)δ S (x)ψ(x) = SΩ dS j κ ij (x)ψ(x) (12)
for any function ψ, with the infinitesimal normal surface vector dS i and the unit normal vector n i oriented towards the interior of the obstacles. Finally, we have also introduced the drift field u * i which can be interpreted as the drift field for the particles upon time reversal, and is given by

u * i (x) = u i (x) -2 J i (x) P 0 (x) , (13)
where we have introduced the current in the steady state

J i J i (x) = u i (x)P 0 (x) -∂ x k [κ ik (x)P 0 (x)]. ( 14
)
The notation ũ * i (x) denotes that we have added surface drift terms which prevent the particles from entering the obstacles,

ũ * i (x) = u * i (x) -n j (x)κ ij (x)δ S (x). (15)
Equation (10) states clearly that the large scale dispersion properties can be obtained from the properties of the propagator (with periodic boundary conditions), which is used to average the drift fields with and without timereversal; this is the first result of the paper.

The second main result is an expression for the late time diffusion tensor D ij . Taking the large time limit of Eq. ( 10) leads to

D ij = Ω dx κ ij (x)P 0 (x) + 1 2 Ω dx [u i (x)f j (x) + u j (x)f i (x)] - 1 2 SΩ dS l (x) [κ il (x)f j (x) + κ jl (x)f i (x)], (16)
where the function f j is defined by

f j (x) = - Ω dy P 0 (y)G(x|y)ũ * j (y) (17)
and

G(x|y) = ∞ 0 dt[P (x, t|y) -P 0 (x)]. (18)
G is the pseudo-Green function [START_REF] Barton | Elements of Green's functions and propagation[END_REF] of the operator H. For computational purposes, it is useful to note that f i can be obtained by solving the non-homogeneous linear partial differential equation

H x f i (x) = -u * i (x)P 0 (x) + P 0 (x) Ω dy u * i (y)P 0 (y) -P 0 (x) SΩ dS l (y)P 0 (y)κ li (y) (19)
with periodic boundary conditions and also with the condition

n j {u j (x)f i (x) -∂ x k [κ kj (x)f i (x)]} = -n j κ ji (x)P 0 (x) (20)
at the surface of the obstacles, together with the orthogonality condition

Ω dx f i (x) = 0. ( 21
)
The above expressions are exact, general and can be used to compute the dispersion properties in any advectiondiffusion problem with periodic diffusivity tensor and drift fields, in any dimension, in the presence of obstacles with reflecting boundaries. We will show that our expression for D ij contains existing expressions for the effective diffusion tensor in porous media for particles transported by incrompressible fluids [14,[START_REF] Carbonell | [END_REF], diffusion in periodic potentials [36], and diffusion in non-equilibrium one dimensional systems [26,28,32]. One may also derive find an alternative expression for D ij by defining a function f * i as

f * i (y) = - Ω dx P 0 (y)G(x|y)ũ i (x), (22)
in terms of which the effective diffusion tensor reads

D ij = Ω dy κ ij (y)P 0 (y) + 1 2 Ω dy [u * i (y)f * j (y) + u * j (y)f * i (y)] - 1 2 SΩ dS l (y)[κ il (y)f * j (x) + κ jl (y)f * i (y)]. ( 23
)
This function f * i can also be obtained as the solution of the partial differential equations where we have denoted H * the transport operator for the time-reversed process,

H * y f * i (y) = -u i (y)P 0 (y) + P 0 (y) Ω dx [u i (x)P 0 (x)] -P 0 (y) SΩ dS l (x)P 0 (x)κ li (x). (24)
H * y (•) = ∂ yi [u * i (y)(•)] -∂ yi ∂ y k [κ ij (y)(•)]. (25)
The function f * i has periodic boundary conditions, and satisfies at the surface of the obstacles

n j (y){u * j (y)f * i (y) -∂ y k [κ kj (y)f * i (y)]} = -n j κ ji (y)P 0 (y), (26)
and also satisfies the orthogonality condition

Ω dy f * i (y) = 0. ( 27
)
Thus, we have two formulas for the effective diffusion tensor, which can be straightforwardly evaluated numerically by solving a set of elliptic partial differential equation on the domain Ω. Finally, we will show that the tensor C ij describing the approach to the diffusive limit can be obtained from

C ij = 1 2 Ω dx f i (x)f * j (x) + f j (x)f * i (x) P 0 (x) . ( 28
)
Although numerous expressions exist in the literature for D ij for particular cases, the correction tensor C ij has to our knowledge only been studied for diffusion in a periodic potential [37] or with periodic diffusivity [38], the above expressions enable its computation in a systematic way and at a computational cost equivalent to that required to calculate D ij .

Our results can be applied to study tracer particles diffusing in a 2D array of circular impenetrable obstacles (hard disks) [Fig. 1(a)] and subjected to an external force. For this system, the elements of the effective diffusion tensor match perfectly with the results of stochastic simulations (the algorithm is presented in appendix A).

III. DERIVATION OF THE KUBO FORMULA FOR THE TIME DEPENDENT DISPERSION

Here we describe the derivation of the Kubo formula for the dispersion tensor. We first briefly compute the average displacement in the direction i during a time t, defined in Eq. ( 3), which is given by

µ i (t) = R d dx Ω dy P 0 (y)p(x, t|y)(x i -y i ). ( 29
)
Note that the integral over the whole space R d denotes the integral over the volume not occupied by obstacles (where p does not vanish). Taking the time derivative of the above expression, and using Eq. ( 1), we obtain

∂ t µ i (t) = R d dx Ω dyP 0 (y){-∂ xj [u j (x)p(x, t|y)] + ∂ xj ∂ x k [κ kj (x)p(x, t|y)]}(x i -y i ). ( 30
)
Integrating by parts over x and applying the divergence theorem yields

∂ t µ i (t) = R d dx Ω dyP 0 (y){u i (x)p(x, t|y) -∂ x k [κ ik (x)p(x, t|y)]}, (31)
where we have taken into account the reflecting boundary conditions at the obstacles, so that no surface integral has appeared. Decomposing the integral over x in integrals over each unit cell, using the periodicity property (7) of u i and κ ij and the definition (9), we get

∂ t µ i (t) = Ω dx Ω dyP 0 (y){u i (x)P (x, t|y) -∂ x k [κ ik (x)P (x, t|y)]}. (32)
We recall that, since P 0 is the steady state particle distribution, it satisfies Ω dyP 0 (y)P (x, t|y) = P 0 (x). (33) Using this property, and again applying the divergence theorem, (32) becomes

∂ t µ i (t) = Ω dx u i (x)P 0 (x) - SΩ dS j (x)κ ij (x)P 0 (x)} = V i , (34)
where S Ω represents the surface of the obstacles inside the cell Ω, and the infinitesimal surface vector dS j is oriented towards the interior of the obstacles. Eq. [START_REF] Barton | Elements of Green's functions and propagation[END_REF] shows that, when the initial conditions are those of the steady state, the average drift V i is constant in time. This is in agreement with the classic result of Stratonovich [START_REF] Stratonovich | Radiotekh Elektron[END_REF].

Let us now derive a Kubo formula for the time dependent dispersion tensor σ ij [defined in Eq. ( 4)]. Define ψ ij to be

σ ij (t) = ψ ij (t) -µ i (t)µ j (t). ( 35
)
By definition, the expression for ψ ij is

ψ ij (t) = R d dx Ω dyP 0 (y)p(x, t|y)(x i -y i )(x j -y j ). ( 36
)
We denote by φ(s) = ∞ 0 dtφ(t)e -st the temporal Laplace transform of any function φ. The Laplace transform of the FP equation ( 1) is given by sp(x, s|y) = δ(x -y)

+ ∂ xi ∂ xj [κ ij (x)p(x, s|y)] -u i (x)p(x, s|y) . (37)
Using the above equality, the expression for

ψ ij becomes ψij (s) = 1 s R d dx Ω dy P 0 (y)(x i -y i )(x j -y j ) ∂ x k {∂ x l [κ kl (x)p(x, s|y)] -u k (x)p(x, s|y)} . ( 38
)
We remark that

∂ x k [(x i -y i )(x j -y j )] = δ ki (x j -y j ) + δ kj (x i -y i ), (39)
so that the integration by parts over x in Eq. ( 38) leads to an expression of the form

ψij (s) = Êij (s) + Êji (s) (40)
with

Êij (s) = 1 s R d dx Ω dy P 0 (y)(x j -y j ) {u i (x)p(x, s|y) -∂ x k [κ ki (x)p(x, s|y)]} . ( 41
)
Integrating by parts the term containing κ ik leads to

Êij (s) = 1 s R d dx Ω dy P 0 (y)(x j -y j )u i (x)p(x, s|y) - 1 s S dS l (x) Ω dyP 0 (y)κ il (x)(x j -y j )p(x, s|y) + 1 s R d dx Ω dyP 0 (y)p(x, s|y)κ ij (x). ( 42
)
Using the periodicity property of κ and the Laplace transform of the steady state property (33), we obtain

R d dx Ω dyP 0 (y)p(x, s|y)κ il (x) = s -1 Ω dx P 0 (x)κ il (x) = s -1 κ ij 0 , (43)
where ••• 0 represents the average over P 0 . Decomposing the integral appearing in (42) over x on all the individual cells of the periodic structure, and changing of variable x → x + k, y → y + k in each of them (where k is the lattice vector such that x + k is inside the cell Ω), and summing over all lattice vectors k again, we see that we can exchange the integration domains between x and y, leading to

Êij (s) = 1 s Ω dxu i (x) Bj (x, s) - 1 s SΩ dS l (x)κ il (x) Bj (x, s) + κ ij 0 s 2 , ( 44
)
where we have defined

Bj (x, s) = R d dy p(x, s|y)P 0 (y)(x j -y j ). (45)
Now comes the key point of our derivation. In order to calculate B j , we first consider another probability distribution q defined as q(y, t|x) = p(x, t|y)P 0 (y)

P 0 (x) . ( 46
)
Using Bayes' theorem, we see that q(y, t|x) is the probability density at the position y at time 0 given that x is the position at the later time t; q is thus the propagator of the tracer particles under time reversal. The evolution of q satisfies

∂ t q(y, t|x)P 0 (x) = -∂ xi [u i (x)P 0 (x)q(y, t|x)] + ∂ xi ∂ xj [κ ij (x)P 0 (x)q(y, t|x)]. (47)
Expanding all the derivatives and using the fact that HP 0 = 0, we find that q satisfies

∂ t q(y, t|x) =[u * i (x)∂ xi + κ ij (x)∂ xi ∂ xj ]q(y, t|x), (48)
where

u * i (x) = 2 ∂ xj [κ ij (x)P 0 (x)] P 0 (x) -u i (x), = u i (x) -2 J 0i (x) P 0 (x) , (49)
is the drift of the time reversed process and where we have introduced the current in the steady state J i

J i (x) = u i (x)P 0 (x) -∂ x k [κ ik (x)P 0 (x)]. ( 50
)
If we interpret (48) as a backward Fokker-Planck equation [START_REF] Gardiner | Handbook of stochastic methods for physics, chemistry and the natural sciences[END_REF], we deduce that q also satisfies the associated forward Fokker-Planck equation

∂ t q(y, t|x) = -∂ yi [u * i (y)q(y, t|x)] + ∂ yi ∂ yj [κ ij (y)q(y, t|x)]. ( 51
)
Therefore q is the propagator of fictive particles, moving in an effective drift field u * i instead of u i . Let us say a little bit more about the properties of q. From its definition (46), we immediately see that the initial condition for q is q(y, 0|x) = δ(x -y).

(

) 52
Given that no particles can flow in or out of the obstacles, the time reversed process must also have no-flux boundary conditions at the surface of obstacles. That is to say that

0 = n i {u * i (y)q(y, t|x) -∂ yj [κ ij (y)q(y, t|x)]} = n i P 0 (x) {u * i (y)P 0 (y) -∂ yj [κ ij (y)P 0 (y)]}p(x, t|y) -κ ij (y)P 0 (y)∂ yj p(x, t|y) , (53)
and as a consequence we see that n i κ ij (y)∂ yj p(x, t|y) = 0, which recovers an established boundary condition for the no-flux boundary condition in terms of the starting coordinate y [START_REF] Gardiner | Handbook of stochastic methods for physics, chemistry and the natural sciences[END_REF]. It is interesting to note that the steady state current of the time reversed process is given by -J i , i.e. exactly the opposite of the current of the original process. In addition, we see that for currentless steady states, where J i = 0, we have that u * i = u i and thus the original and time reversed processes are statistically identical in that p(x, t|y) = q(x, t|y). Interestingly, for advection of a particle with constant molecular diffusivity by an incompressible flow, one can easily show that P 0 = 1/|Ω| where |Ω| is the volume of the unit cell and that J i is non-zero and given by J i = u i /|Ω|, consequently u * i = -u i ,i.e. the time reversed process has the opposite flow field to the original process.

Eq. ( 45) may be rewritten in terms of q as Bj (x, s) =

R d dy q(y, s|x)P 0 (x)(x j -y j ). ( 54
)
Laplace transforming the Fokker-Planck equation (51) for q yields

s q(y, s|x) = δ(x -y) -∂ yi [u * i (y)q(y, s|x)] + ∂ yi ∂ yj [κ ij (y)q(y, s|x)]. ( 55
)
Inserting the above equality into Eq. (54) and integrating by parts over y, we obtain Bj (x, s) = 1 s R d dy P 0 (x){-u * j (y)q(y, s|x)

+ ∂ y k [κ kj (y)q(y, s|x)]}. ( 56
)
where we have used the no-flux condition (this time of the time reversed process) at the obstacles boundaries.

Integrating by parts the term containing κ kj , and then switching back to the propagator p instead of q and invoking the periodicity property, we obtain Bj (x, s) = -1 s Ω dy P 0 (y) P (x, s|y)u * j (y) + 1 s SΩ dS l (y)P 0 (y) P (x, s|y)κ lj (y). (57)

Inserting the above expression into Eqs. ( 44), (36) gives a formal expression for the tensor ψij . Using the fact that P (x, s|y) P 0 (x)/s for s → 0, one can check that, for large times, ψ ij µ i (t)µ j (t). It is now useful to introduce the difference between the propagator and its stationary value

P (x, t|y) = P (x, t|y) -P 0 (x), (58)
and thus has Laplace transform P (x, s|y) = P (x, s|y) -

1 s P 0 (x). ( 59
)
In terms of P , we obtain from Eqs. ( 35),( 36), (44), and (57) the following final expression for the temporal Laplace transform of the dispersion tensor

σ ij (t) σij (s) = 2 s 2 Ω dx κ ij (x)P 0 (x) + 1 s 2 Ω dx[u i (x)K j (x, s) + u j (x)K i (x, s)] - 1 s 2 SΩ dS l (x)[κ il (x)K j (x, s) + κ jl (x)K i (x, s)], (60)
with K j (x, s) = -Ω dy P 0 (y) P (x, s|y)u * j (y) + SΩ dS l (y)P 0 (y) P (x, s|y)κ lj (y). ( 61)

and we recall that u * i is the effective drift field after time reversal symmetry defined by Eq. (49). Equations ( 49),( 60) and (61) give a closed form expression for the Laplace transform of the dispersion tensor σ ij (t) at all times and is the key result of the paper.

IV. THE EFFECTIVE DIFFUSION TENSOR

Here we extract the late time effective diffusion tensor D ij . The Eqs. (60) and (61) are written in Laplace space and the large time limit is thus extracted from the small s behavior, for which

σij (s) 2 D ij s 2 (s → 0). ( 62
Taking s = 0 in the above expression, we find that P is finite at s = 0 (as P (x, t|y) → P 0 (x) as t → ∞). Therefore, we can define a function G such that

P (x, 0|y) = ∞ 0 dt [P (x, t|y) -P 0 (x)] ≡ G(x|y). ( 64
)
Applying the operator H to the above equality, it is clear that G obeys

H x G(x|y) = δ(x -y) -P 0 (x), (65)
and also from the conservation of probability G satisfies

Ω dx G(x|y) = 0. ( 66
)
The function G is therefore the pseudo-Green's function [START_REF] Barton | Elements of Green's functions and propagation[END_REF] of the operator H (that is, the inverse of H in the subspace orthogonal to the uniform function). Since P (x, s = 0|y) is finite, we see that the inverse Laplace transform at late times can be extracted from Eqs. ( 60) and (61) upon setting s = 0 in P . This immediately yields

D ij = Ω dx κ ij (x)P 0 (x) + 1 2 Ω dx [u i (x)f j (x) + u j (x)f i (x)] - 1 2 SΩ dS l (x)[κ il (x)f j (x) + κ jl (x)f i (x)], (67)
For computational purposes, it is useful to find the partial differential equation satisfied by f i . Acting with H on the above expression and using Eq. (65) leads to

H x f i (x) = -u * i (x)P 0 (x) + P 0 (x) Ω dy [u * i (y)P 0 (y)] -P 0 (x) SΩ dS l (y)P 0 (y)κ li (y), (69)
which is valid for x ∈ Ω (outside the obstacles). An equivalent equation for f i is

H x f i (x) = -u * i (x)P 0 (x) + P 0 (x) Ω dy {u * i (y)P 0 (y) -∂ y k [P 0 (y)κ kj (y)]}. ( 70
)
This equation for f i (x) (for x outside the obstacles) is supplemented by the following conditions. First, because of the conservation of probability, it follows from Eq. (66) that

Ω dx f i (x) = 0. (71)
Second, f i (x) is periodic on Ω. Third, at the surface obstacles, f i satisfies

n j {u j (x)f i (x) -∂ x k [κ kj (x)f i (x)]} = -n j κ ji (x)P 0 (x) (72)
It is actually not obvious to derive the boundary conditions (72) from Eq. ( 68). The validity of these boundary conditions is checked explicitly in Appendix B, where we demonstrate that Eq. ( 68) is the actual solution of the Eqs. (70,71,72), proving that our formulation is correct. One may alternatively express the diffusion tensor by defining a function f * i as

f * i (y) = - Ω dx P 0 (y)G(x|y)u i (x) + SΩ dS l (x)P 0 (y)G(x|y)κ li (x), ( 73
)
in terms of which the effective diffusion tensor reads

D ij = Ω dy κ ij (y)P 0 (y) + 1 2 Ω dy [u * i (y)f * j (y) + u * j (y)f * i (y)] - 1 2 SΩ dS l (y)[κ il (y)f * j (y) + κ jl (y)f * i (y)]. ( 74
)
Obviously,

f * i (y) = - Ω dx P 0 (x)G * (y|x)u i (x) + SΩ dS l (x)P 0 (x)G(y|x)κ li (x), ( 75
)
where G * is the pseudo-Green function associated for the operator H * of the time reversed process, with the drift field u * , that reads

H * y (•) = ∂ yi [u * i (y)(•)] -∂ yi ∂ y k [κ ij (y)(•)] (76)
Hence, f * i satisfies the same equations as f i if one replaces u i by u * i and vice-versa. As a consequence, f * i is the solution of

H * y f * i (y) = -u i (y)P 0 (y) + P 0 (y) Ω dx [u i (x)P 0 (x)] -P 0 (y) SΩ dS l (x)P 0 (x)κ li (x) (77)
with periodic boundary conditions and also with the condition

n j {u * j (y)f * i (y) -∂ y k [κ kj (y)f * i (y)]} = -n j κ ji (y)P 0 (y) (78)
at the surface of the obstacles, together with the orthogonality condition

Ω dy f * i (y) = 0. ( 79
)
Thus, we have two formulas for the computation of the effective diffusion tensor.

Having two formulations for the Kubo formula is a useful check on numerical solutions as the two resulting results for the diffusion tensor can be compared. Furthermore, we will now see that both the fields f i and f * i are needed to compute the leading order finite time correction to the time dependent diffusion tensor.

V. LATE TIME CORRECTIONS TO THE EFFECTIVE DIFFUSION TENSOR

In section IV we have extracted, from the full time dependent Kubo formula, the effective late time limit of the diffusion tensor D ij . While many examples of such formulas exist for special cases, little is known about how the time dependent diffusion tensor relaxes to its asymptotic limit. To our knowledge the only examples known are for diffusion in a periodic potential [37] and diffusion in a system with periodic diffusivity [38]. Both of these examples have steady states with zero current, here we will extend these results to the most general cases both with and without current. Here we extract the late time correction tensor C ij that describes the approach to the diffusive limit of the system.

Expanding Eq. ( 63) in powers of s, we obtain P (x, s|y) G(x|y) + s G (2) (x|y) + ...

with

G (2) (x|y) = - ∞ 0 dt t [P (x, t|y) -P 0 (x)]. ( 81
)
Since the motion of the tracer particles is Markovian (memoryless), we can write the equality

P (x, t|y) = Ω dz P (x, t -t 1 |z)P (z, t 1 |y), (82)
which holds for 0 < t 1 < t. Integrating this relation over t 1 leads to

t P (x, t|y) = t 0 dt 1 Ω dz P (x, t -t 1 |z)P (z, t 1 |y). ( 83
)
Substracting tP 0 (x) and taking the Laplace transform for small values of s, we obtain

G (2) (x|y) = - Ω dz G(x|z)G(z|y). ( 84
)
The correction tensor C ij is identified from the relation for small s

σij (s) 2 s 2 (D ij + s C ij + ...) (s → 0). ( 85
)
Inserting this relation into Eqs. ( 60) and (61), and using (80) and (84) leads to

C ij = - 1 2 Ω dz Ω dx[u i (x)f j (z) + u j (x)f i (z)]G(x|z) - SΩ dS l (x)[κ il (x)f j (z) + κ jl (x)f i (z)]G(x|z) , ( 86
)
where we have used the definition of f in Eq. ( 68) to perform the integration over y. Using the definition of f * i in Eq. ( 75) and the relation G * (y|x)P 0 (x) = G(x|y)P 0 (y), we finally obtain

C ij = 1 2 Ω dz f i (z)f * j (z) + f j (z)f * i (z) P 0 (z) . ( 87
)
This relation is a compact and explicit form for the correction tensor C ij and is the main result of this section. One can also derive the above result by decomposing P in terms of left and right eigenvector of the Fokker-Planck operator H. In doing this is is straightforward to see that the temporal corrections at the next order decay as exp(-λ 1 t), where λ 1 is the lowest positive eigenvalue of H which will be strictly positive given that the domain Ω is taken to be finite.

VI. SPECIAL CASES AND COMPARISON WITH EXISTING RESULTS

Here we examine the form the Kubo formula takes for a number of systems and compare them to existing results in the literature.

A. Flow in frozen incompressible velocity fields with isotropic constant diffusion tensor

The Fokker-Planck transport operator for a tracer advected by an incompressible velocity field v with constant isotropic diffusivity κ ij (x) = κ 0 δ ij is defined via

H x = -κ 0 ∂ xi ∂ xi + v i (x)∂ xi . (88)
Here the drift field is thus u i = v i . From the incompressibility condition ∂ xi v i = 0, we see that the steady state distribution on the unit cell Ω is uniform and thus given by P 0 (x) = 1/|Ω|, where |Ω| denotes the available volume of the unit cell. The steady state current is thus given by J i (x) = v i (x)/|Ω|. We also find that u * i = -u i . Applying Eq. (67) leads to

D ij = δ ij κ 0 + 1 2 Ω dx [f i (x)v j (x) + f j (x)v i (x)] - κ 0 2 SΩ dS i (x)f j (x) + SΩ dS j (x)f i (x) . (89)
The equation satisfied by f i is a simplified form of Eq. (70):

-κ 0 ∂ 2 xj f i + v j ∂ xj f i = 1 |Ω| v i (x) - 1 |Ω| 2 Ω dy v i (y). (90)
The boundary conditions on the surface of obstacles are [see Eq. ( 72)]

n j ∂ xj f i = n i /|Ω|, (91)
(where we have assumed that no fluid enters the obstacles, v i n i = 0). It is then obvious that the last four expressions recover those in the hydrodynamicshomogenization literature (see e.g. Refs. [14,[START_REF] Carbonell | [END_REF]), if we identify the function f i of this paper with -f i /|Ω| to match with the notations of Ref. [START_REF] Carbonell | [END_REF].

B. Systems with no steady state current

A large class of models studied in statistical mechanics, such as diffusion in a periodic potential or diffusion in a medium having periodic diffusivity, have no current in the steady state. In these cases the original and its timereversed processes have the same Fokker-Planck evolution operator, and consequently the same pseudo-Green's function. However in general we have that

G(x|y)P 0 (y) = P 0 (x)G * (y|x), (92)
where G * (y|x) denotes the pseudo-Green's function for H * . In the case where there is no current we know however that G * (y|x) = G(x|y) and thus we find that the operator

K(x, y) = G(x|y)P 0 (y) = P 0 (x)G(y|x) (93)
is symmetric. In the currentless case we also have

D ij = Ω dx κ ij (x)P 0 (x) - 1 2 Ω dxdy[ũ i (x)ũ j (y) + ũj (x)ũ i (y)]K(x, y), (94)
from Eq. ( 10). The operator K, being symmetric, must have real eigenvalues, however K has the same eigenvalues as G which must have eigenvalues with a positive real part if a steady state regime can be attained. Consequently under these assumptions K must be a positive operator and we must have that the change in the diagonal terms of the diffusion tensor due to the presence of drift or inclusions is negative, i.e.

D ii -κ ii 0 = - Ω dxdyũ i (x)ũ i (y)K(x, y) < 0, (95)
that is to say that diffusion in equilibrium systems is slowed down by the presence of the drift. This is physically obvious for say diffusion in periodic potentials where the diffusing particle becomes trapped in local minima of the potential. This slowing down due to drift also occurs for diffusion in a medium of spatially varying diffusivity where

H x = -∂ xi κ ij (x)∂ xj , (96)
here P 0 = 1/|Ω| and J i = 0. The drift field in this case is, in our notation, given by

u i (x) = ∂ xj κ ij (x), (97)
and in terms of the pseudo-Greens' function G we find

D ij = 1 |Ω| Ω dx κ ij (x) - 1 |Ω| Ω dxdy [∂ x k κ ik (x)][∂ x l κ il (y)]G(x|y), ( 98
)
which recovers the formula derived in Ref. [21]. We see that the first term in the above is the spatial average of the local diffusivity tensor while the second terms gives a negative contribution when i = j. Finally we remark that the case of diffusion with constant diffusivity κ 0 in a periodic potential Φ gives u i = -κ 0 β∂ xi Φ(x), where β is the inverse temperature. Here the steady state distribution P 0 is given by the Gibbs-Boltzmann distribution on the unit cell Ω

P 0 (x) = exp [-βΦ(x)] Z Ω , (99)
where Z Ω is the partition function normalizing the distribution, thus recovering the result of Ref. [21] for this particular case.

VII. GENERAL RESULTS IN ONE DIMENSION

We consider a generic one dimensional system with periodicity denoted by L. Here we can give explicit formula for all terms in the static Kubo formula, although the computation is surprisingly involved. In what follows we will use notation similar to Refs. [26,28] to facilitate comparison with their results. In one dimension the steady state current is constant in space and the steady state probability distribution is given by

P 0 (x) = JI + (x), (100)
where the term I + reads

I + (x) = exp (Γ(x)) κ(x) ∞ x dx exp (-Γ(x )) , (101) with Γ
(x) = x 0 dx u(x) κ(x) . ( 102
)
Due to the periodicity of u and κ the function Γ obeys the relation

Γ(x + L) = Γ(x) + Γ(L). ( 103
)
When Γ(L) = 0 the system clearly has a steady state equilibrium distribution with no current. In writing Eq. (101) we have assumed, without loss of generality, that Γ(L) > 0 so that the integral on the right hand side converges. The steady state current is then obtained from the condition of normalization of P 0 and is thus given by

J = 1 L 0 dx I + (x) , (104)
and consequently the effective drift is given by V = JL.

To compute the effective diffusion constant we use the representation for D given in Eq. ( 74) in its one dimensional version, that is to say

D = L 0 dx κ(x)P 0 (x) + L 0 dxf * (x)u * (x).
We now define f * (x) = g(x)P 0 (x) and write the diffusion constant as

D = L 0 dx κ(x)P 0 (x)+ L 0 dxg(x) J - d dx [κ(x)P 0 (x)] .
Defining g(x) = -x + h(x) we find that h obeys

κ d 2 h dx 2 + u dh dx = JL. (105)
A first integration of this equation for h yields

dh(x) dx = JLI -(x), (106)
where

I -(x) = exp (-Γ(x)) x -∞ dx exp (Γ(x )) κ(x ) , (107)
and where we have used the fact that dg/dx = -1+dh/dx is periodic and thus dh/dx must be periodic. Integrating again yields

h(x) = J 0 L x 0 dx I -(x ) + C ( 108
)
where C is an integration constant that will be determined from the orthogonality condition L 0 dx P 0 (x)g(x) = 0. It is not immediately obvious from Eq. (108) that g(x) = -x + h is periodic, however it is easy to see that

g(x + L) = g(x) -L + JL L 0 dx I -(x). ( 109
)
Now by integration by parts one finds

L 0 dx I -(x) = L 0 dx I + (x), (110)
and using this and Eq. (104) in Eq. (109) yields g(x) = g(x + L). Incidentally this shows that we may also write

J = 1 L 0 dx I -(x) . ( 111
)
To proceed it is useful to define the two new functions

A + (x) = x -∞ dx exp (Γ(x )) κ(x ) , (112)
A -(x) = ∞ x dx exp (-Γ(x )). ( 113
)
Further more one can show that

A + (x + L) = exp (Γ(L)) A + (x), ( 114
) A -(x + L) = exp (-Γ(L)) A -(x). ( 115
)
In terms of these two functions we can write I + and I - as

I + (x) = dA + (x) dx A -(x), (116)
I -(x) = - dA -(x) dx A + (x). ( 117
)
Computing the constant C in Eq. ( 108) using the orthogonality relation then yields

D = L 0 dx κ(x)P 0 (x) dh dx + J [-x + h(x) + LxP 0 (x) -Lh(x)P 0 (x)] , (118)
where we have carried out an integration by parts in the last term of Eq. (105). The function h is explicitly given by

h(x) = JL x 0 dx I -(x ), (119)
and using the expression Eq. ( 117) for I -it can be written as

h(x) = -JL A + (x)A -(x) -A + (0)A -(0) - x 0 dx P 0 (x ) J . (120)
We also, by integration by parts, note the identity

L 0 dx x 0 dx P 0 (x ) = L - L 0 dx xP 0 (x). ( 121
)
Using these relations in Eq. ( 118), after some algebra, we obtain the compact expression

D = L 2 L 0 dx κ(x)I ± (x) 2 I ∓ (x) L 0 dx I ± (x) 3 , (122)
where ± indicates that the index can be taken to be + or -. One can verify that Eq. ( 122) agrees with the results of Reimann et al. [26,28] for the case of diffusion in a tilted potential φ(x) = V (x) -f x, where V is periodic. It also agrees with the formulas presented in Ref. [31] where the local mobility is also varying. The result given here shows how the diffusion constant of any one dimensional system with periodic diffusivity and drift can be obtained.

VIII. CONCLUSIONS

The problem of computing effective transport coefficients occurs in many settings of fundamental physics and has many applications. Many results exist in the literature, for example in fluid mechanics using homogenization theory [10][11][12][13][14][15][16][17][18][19][20] as well as in statistical physics [21][22][START_REF] Zwanzig | Proc. Natl. Acad. Sci. U S A[END_REF][START_REF] Gennes | [END_REF][25] where both exact and approximate results exist for problems such as diffusion in periodic or random potentials and diffusion in media with locally varying diffusivity. In these latter studies most results have been obtained for systems with equilibrium (currentless) steady states, given for example by the Gibbs Boltzmann distribution for diffusion in a periodic potential.

More recently results have been derived for diffusion in systems with non zero currents [26][27][28][29][30][31][32] In such non equilibrium systems interesting new physics such as massively increased diffusivity due to an external applied force has been discovered. This enhancement occurs, in the weak noise limit, at the external field where all minima disappear from the overall potential to be replaced by points of inflection [26]. The difference between different trajectories at this critical point is enhanced and dispersion in thus increased. Using the Kubo formula, given here, one could investigate the effect of an externally applied field on diffusion in higher dimensional periodic potentials. The formula given here cannot in general be evaluated analytically in higher than one dimension, however their numerical evaluation is straightforward and the transport coefficients can be accurately determined using standard packages to solve partial differential equations. Of course many of the systems considered here can be studied via numerical simulations where one integrates the corresponding SDE, however the simulation approach suffers from the need to estimate the errors due to statistical fluctuations and in the case of obstacles the correct imposition of the no-flux boundary condition requires careful treatment and is far from being obvious [START_REF] Peters | [END_REF][41][42].

Many other problems can be tackled using the present approach, in Ref. [33] it was shown how the presence of an externally applied uniform field on a system with spatially varying diffusivity modifies the dispersion of a cloud of tracer particles. It was shown that the diffusion in the direction of the applied force could be hugely increased for systems in two or more dimensions. In addition it was shown that the components of the effective diffusion tensor can exhibit counter intuitive monotonic behavior. The results thus imply that one can control the dispersion of a diffusing cloud with an external field, and such an effect may well have useful applications.

The formula given here are extremely general and should be valuable for the interpretation of experiments. For instance, potential fields in which colloidal particles diffuse are often generated by laser [43][44][45][46], in which case absorption of laser light can also lead to temperature gradients. The variation of the local temperature will clearly influence local transport properties through the temperature dependence of the surrounding fluid viscosity as well as due to thermodynamic forces associated with temperature gradients, i.e the Soret effect [5].

Finally there has been much recent study of fluctuation dissipation relations (FDR) in non-equilibrium systems [47][48][49][50]. A notable example of this type of FDR is the Stokes-Einstein relation between the differential response of the mean velocity V i and the diffusion tensor. In our preliminary study [33], we showed how the presence of a steady state current leads to corrections to the usual Stokes-Einstein formula. However the general analysis of the correction term needs further study in order to fully understand the non-equilibrium physics it encodes, as well as to make contact with the literature on the subject [47][48][49][50]. In the literature in question results are most often given in terms of an integral over a time dependent violation factor, the results given here could be useful to understand these results in terms of the static, time independent quantities, used here.

Appendix A: Details on the numerical simulations Here we briefly describe the algorithm used to simulate the motion of a Brownian walker in a 2D array of circular disks, in the presence of a force F , leading to the results presented in Fig. 1(b). Let X be the position at time t of the walker, then we define e α the unit vector with direction X -X c , with X c the nearest disk center. e α makes an angle α with the x-direction. During the time step ∆t, the motion in the direction e α is modified by the increment

dX ⊥ = √ ∆tf 1 (r ⊥ / √ ∆t) + u f 2 (r ⊥ / √ ∆t) + h cos α∆t (A1)
Here h = βF D, r ⊥ is the distance to the nearest disk surface, and f 1 and f 2 are the functions calculated by Peters et al. in the case of reflecting boundaries [START_REF] Peters | [END_REF], and u = ± √ D∆t with probabilities 1/2. In the direction parallel to the nearest surface obstacles, the increment is dX = -h sin α∆t + u, where u is a stochastic variable with zero mean and variance u 2 = 2D∆t. With this algorithm, when the particle is close to the boundaries, r ⊥ ∼ √ ∆t, the terms in f 1 and f 2 dominate the increment in the perpendicular direction in Eq. (A1) if ∆t is small enough. These terms account for the fact that, since trajectories that cross the surface of the obstacles are forbidden, the averaging over authorized trajectories between t and t + ∆t gives rise to an effective drift oriented towards the exterior of the obstacles, described by the term f 1 , and also modifies the variance of these trajectories (term f 2 ), see Ref [START_REF] Peters | [END_REF] for details. If the particle is far from the obstacles, r ⊥ √ ∆t, then f 1 → 0 and f 2 → √ 2 and one recovers Brownian motion under force in the absence of obstacles. Although the algorithm of Ref [START_REF] Peters | [END_REF] was proposed only for one-dimensional geometries, we expect that it is also valid in the present 2D case, since for sufficiently small ∆t the curvature of the surface should play no role. For the results of Fig. 1(b), the time step was ∆t = 10 -5 L 2 /D and was checked to be sufficiently small so that convergence was reached for the calculation of D ij which were estimated over 500, 000 trajectories. The results of simulations match with the predictions of our (exact) Kubo formulas.

FFIG. 1 :

 1 FIG.1:(a) Examples of stochastic trajectories of two tracer particles diffusing in a heterogeneous medium (here a 2D array of circular obstacles) under a force F (arrow) that biases the motion to the right. (b) Elements of the dispersion tensor for the system described in (a), calculated from our theoretical formulas Eqs.(16,19) (full and dashed lines) vs results of stochastic simulations (symbols). Parameters: disk radius R = 0.25L with L the length of the side of the periodic cell. Results for Dij are expressed in units of the (constant) diffusion coefficient of the tracers outside the obstacles.

  st [P (x, t|y) -P 0 (x)].

  with f j (x) = -Ω dy P 0 (y)G(x|y)u * j (y) + SΩ dS l (y)P 0 (y)G(x|y)κ lj (y).

Appendix B: Boundary condition (72) for the function fi

In this appendix, we show that the formulation of the problem under the form of Eqs. (70,71,72) is correct by checking that Eq. ( 68) is the actual solution of these equations. We consider the quantity

where the adjoint operator reads

The pseudo-Green function G can be shown to satisfy the adjoint equation

The following condition also holds

Using the above relations and Eqs. ( 69) and (71) we deduce that

On the other hand, using the explicit expressions of H [Eq.( 1)] and integrating by parts two times the term containing H in Eq. (B1), we obtain

The last term of this equation vanishes, since the relation dS j κ kj (y)P (x, t|y) = 0 is known for reflecting boundaries [START_REF] Gardiner | Handbook of stochastic methods for physics, chemistry and the natural sciences[END_REF]. Then, using the boundary conditions (72), we obtain

which, when compared to Eq. (B5), is exactly Eq. ( 68), thereby proving that the formulation of Eqs. (70,72,71) is correct.


Alain Burgisser
email: burgisse@cnrs-orleans.frbscaille@cnrs-orleans.fr

Bruno Scaillet

Redox evolution of a degassing magma rising to the surface

Volatiles carried by magmas, either dissolved or exsolved, exert a fundamental role on a variety of geological phenomena such as magma dynamics 1-5 or the composition of Earth's atmosphere 6 . In particular, the redox state of volcanic gases emanating at the Earth surface is widely believed to mirror that of their magma source, and thought to have exerted a first order control on the secular evolution of atmospheric oxygen 6-7 . However, fO 2 estimates retrieved from lava or gas chemistry may differ by as much as one log unit 8-10 , the reason of such a difference remaining obscure. Here, by using a coupled chemical-physical model of conduit flow, we show that the redox state evolution of an ascending magma, and thus of its coexisting gas phase, is strongly dependent on both the composition and amount of gas in the reservoir.

Magmas with no sulphur show a systematic fO 2 increase during ascent, by as much as

2 log units. Magmas with sulphur show also a change of redox state during ascent, but the direction of change depends on the initial fO 2 in the reservoir. Our calculations closely reproduce the H 2 S/SO 2 ratios of volcanic gases observed at convergent settings, yet the difference between fO 2 in the reservoir and that at the exit of the volcanic conduit may be as much as 1.5 log unit. This fO 2 evolution provides one explanation for the apparent conflict documented between lava and gas based fO 2 estimates. We thus conclude that the redox state of erupted magmas is not necessarily a good proxy of the redox state of the gases they emit. Our findings may require re-evaluation of current models aimed at quantifying the role of magmatic volatiles on geological processes.

Recent studies have investigated the role of redox equilibria during volatile exsolution but the solution of the numerical problem requires fixing either iron or sulphur redox states [START_REF] Moretti | On the oxidation state and volatile behavior in multicomponent gas-melt equilibria[END_REF] : as a result, the effect of decompression on oxygen fugacity (fO 2 ) cannot be evaluated. In the present work, we relax such an assumption by taking advantage of the fact that the solubility laws of key volatile species other than H 2 O and CO 2 have been recently determined, in particular that of H 2 , SO 2 , and H 2 S [START_REF] Clemente | The solubility of sulphur in hydrous rhyolitic melts[END_REF][START_REF] Gaillard | Rate of hydrogen-iron redox exchange in silicate melts and glasses[END_REF] . In this work we considered gas phases in the system H-O-S, with 6 species (H 2 O, H 2 , SO 2 , H 2 S, S 2 , and O 2 ) in which the following redox equilibria occur:

H 2 + ½ O 2 = H 2 O
(1)

H 2 S + ½ O 2 = H 2 O + ½ S 2 (2) O 2 + ½ S 2 = SO 2 (3)
Standard thermodynamic considerations [START_REF] Holloway | Thermodynamic modelling of geological materials: minerals, fluids and melts[END_REF][START_REF] Scaillet | Role of fO 2 on fluid saturation in basalts[END_REF] show that once total pressure, P, temperature, T, and two additional intensive parameters such as water and hydrogen fugacities (fH 2 O and fH 2 ) are known, the fugacities of all remaining species (fH 2 S, fO 2 , fS 2 , fSO 2 ) are fixed and the gas phase composition in the H-O-S system is fully determined. Each species i in the gas has a mole fraction m i :

Σm i = 1
The gas phase is modelled as an ideal mixture of non-ideal gases, a valid approximation in the pressure range considered here. Departure from ideal behaviour of end-member species is accounted for by the fugacity coefficient γ i which is fixed by P and T [START_REF] Shi | Thermodynamic modeling of the C-H-O-S fluid system[END_REF] . The total weight fraction of each species (w Ti ) is the sum of its exsolved part (w gi ) and its dissolved part:

w Ti = w gi + a i (γ i m i P) bi where a i and b i are solubility constants experimentally determined (Table 1). We use a homogeneous, one-dimensional conduit flow model [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF] to simulate magma ascent under closed system and equilibrium conditions. Magma rises in a cylindrical conduit at constant mass flux and volatile exsolution affects the flow through changes in buoyancy and viscosity (cf. Supplementary Information). Calculations are performed by first fixing P, T, fH 2 O fH 2 , and the amount of gas in the reservoir. The model seeks first the equilibrium distribution of each volatile species between gas and melt prior to ascent. Then, at each next lower pressure, mass conservation requires that total amounts of each element (O, H, and S) remain constant. Using this constraint, the model calculates the equilibrium distribution of volatile species, which in turn affects ascent dynamics.

Redox equilibrium during magma ascent involving dissolved iron can be written such as:

Fe 2 O 3 = 2 FeO + ½ O 2 (4)
The importance of such a reaction will be dictated by the initial amount of dissolved iron oxides. It is currently not possible to include iron as a species into our thermodynamic model because the activities of FeO and Fe 2 O 3 are not known. We thus simulated the redox effect of iron by using an empirical model [START_REF] Kress | The compressibility of silicate liquids containing Fe2O3 and the effect of composition, temperature, oxygen fugacity and pressure on their redox states[END_REF] . Runs under typical storage conditions of arc rhyolites (≤1 wt% total iron and fO 2 between NNO-1 and NNO+1, fO 2 being referenced to the solid buffer Ni-NiO, such that NNO+1 means an fO 2 one order of magnitude higher than NNO) show that reaction (4) partly buffers changes in fO 2 when little gas is present in the reservoir (<0.2 unit log without changing the redox trend during ascent, see SI). In contrast, in iron-rich liquids such as basalts, the buffering capacity of iron species will be higher. Thus our results primarily apply to magmas in which the residual melt is rhyolitic, as commonly observed in arc settings. We did not consider the role of iron in crystals because the role of solid buffers as a first order mechanism controlling redox state during magma ascent can be ruled out on kinetic grounds. The lattice diffusion processes that this mechanism requires are exceedingly slow compared to those in gas or liquid phases. Similarly, for simplicity we have not investigated the role of CO 2 since detailed petrologic studies of silicic magmas in arcs show them to have little or no CO 2 under pre-eruptive conditions [START_REF] Wallace | Volatiles in subduction zone magmas: concentrations and fluxes based on melt inclusion and volcanic gas data[END_REF] . In general, however, introduction of CO 2 will lower fH 2 O and thus drives our calculated values toward lower fO 2 via equilibrium (1). Our assumption of equilibrium implies that the model might not capture the chemical evolution of the gas within rapidly decompressed magmas such as in Plinian eruptions, in which the contrasted diffusive kinetics of volatile species may inhibit attainment of equilibrium. Both H 2 O and H 2 are, however, fast diffusing species relative to CO 2 and S-bearing ones [START_REF] Wallace | Volatiles in subduction zone magmas: concentrations and fluxes based on melt inclusion and volcanic gas data[END_REF][START_REF] Watson | Diffusion in volatile-bearing magmas[END_REF] . Thus, if physical fractionation of volatile species arises during ascent, the gas phase composition will be driven toward the system H-O (C and S species remain in the melt, such as documented for S for the 1991 Pinatubo eruption [START_REF] Westrich | Magmatic gas source for the stratospheric SO 2 cloud from the June 15, 1991 eruption of Mount Pinatubo[END_REF] ), which will exert a dominant control on the redox state of escaping gases. Our model represents thus a fundamental end-member case toward which magmas, notably those intermediate to silicic in arc settings, tend to evolve.

We have explored the following range of starting conditions typical of silicic arc magmas [START_REF] Scaillet | Experimental constraints on volatile abundances in arc magmas and their implications fro degassing processes[END_REF][START_REF] Scaillet | Petrological and volcanological constraints on volcanic sulfur emissions to the atmosphere[END_REF] : an initial pressure from 2000 to 3000 bars, fO 2 from NNO+2 to NNO-0.5, bulk iron contents (FeO*) up to 3 wt%, bulk water contents up to 10.4 wt% and bulk sulphur contents up to 3 wt%, the two latter parameters being adjusted by varying the amount of excess gas in the reservoir (up to 5 wt% [START_REF] Wallace | Quantification of pre-eruptive exsolved gas contents in silicic magmas[END_REF] ). The conduit radius was fixed at 5 m for all runs. Temperature has been fixed to 825°C, i.e. typical of rhyolite magmas [START_REF] Scaillet | Experimental constraints on volatile abundances in arc magmas and their implications fro degassing processes[END_REF] , and melt density to 2140 kg/m 3 , though different choices will not affect the trends observed.

Runs are constrained to reach atmospheric pressure at the vent, which yields initial ascent rates between 0.7 and 12 m/s. Because degassing occurs in equilibrium, changes in chemistry as a function of pressure occur regardless of ascent speed. In all cases, the simulations are carried down to atmospheric pressure though it can be anticipated that the last increments of fO 2 change we compute might not be reached owing either to the inhibiting effect of viscosity at low water content [START_REF] Gardner | Bubble growth in highly viscous silicate melts during continuous decompression from high pressure[END_REF] or to gas loss once a permeability threshold is reached such that the system becomes open to gas [START_REF] Eichelberger | Non explosive silicic volcanism[END_REF] .

We first consider the S-and Fe-free case. All simulations performed showed essentially identical behaviour in terms of redox state evolution, i.e. the fO 2 of the magma increases during ascent, the magnitude of increase being more or less dampened by the amount of free gas present in the reservoir. A representative example is shown on Figure 1, corresponding to a magma initially stored at 2000 bars, at an fO 2 of NNO-0.2, with a melt H 2 O content of 5.46 wt% and for excess gas contents of 0.1, 1, and 5 wt%. The simulations show that the magnitude of fO 2 change increases when the amount of excess gas decreases, and can reach more than two log units for a gas-poor magma (0.1 wt%).

We now explore the case of sulphur-bearing rhyolite magma with 1 wt.% of total iron and stored at 2000 bar, illustrating our calculations with three different initial fO 2 but similar initial dissolved sulphur contents (ca. 200-250 ppm). A magma starting at an fO 2 of NNO-0.5 displays a continuous increase in its fO 2 as it ascends, except in the last few hundred bars where a reversal in fO 2 toward reduction occurs (Fig. 2). The magnitude of change is strongly dependent on the amount of gas initially present in the reservoir. At low gas content (0.1 wt%), the fO 2 rises by 0.7 log unit relative to starting conditions. With 5 wt% gas, the magma has a redox state nearly constant up to a pressure of 100 bars. When the starting fO 2 in the reservoir is at NNO (Fig. 2), the magma undergoes oxidation only for gas-poor conditions (0.1 wt%). Higher amounts of gas in the reservoir impart a reducing trend in the fO 2 evolution during ascent, the final fO 2 differing by almost one log unit from the initial value for an initial gas content of 5 wt%. When the initial fO 2 in the reservoir is at NNO+1.5, the magma undergoes a significant reduction during ascent regardless of its initial gas content (Fig. 2). In this case, the drop in fO 2 may exceed 1.5 log units at near atmospheric conditions for gas rich conditions.

The change in fO 2 during decompression is accompanied by dramatic changes in gas phase composition (Fig. 3). Our calculated H 2 S/SO 2 ratios fall in the range 0.1-10, which is comparable to that of volcanic gases measured at convergent settings for silicic to intermediate magmas [START_REF] Symonds | Volcanic-gas studies: methods, results, and applications[END_REF] (Fig. 3). Clearly, a variety of H 2 S/SO 2 ratio can be produced from magmas having common initial redox states but different amounts of gases. Conversely, a given H 2 S/SO 2 ratio may be produced from a wide range of starting redox conditions. For instance a H 2 S/SO 2 ratio of about one can be produced from a magma initially stored at NNO+1.5 with 5 wt% gas, or from a magma initially at NNO-0.5 with 0.1 wt% gas (Fig. 3a). Thus, anticipating the H 2 S/SO 2 ratio of gases emanating from a given reservoir would require knowing not only the reservoir fO 2 but also the amount of free gas present in the reservoir and the depth at which gas and melt are physically separated. Although such a rich behaviour precludes a simple explanation for each trend calculated, a sensitivity analysis of our model suggests that water exsolution plays an important role in oxidizing the system, and that, on the other hand, the complex pressure dependence of the redox equilibria ( 1)-( 3) contribute to the reducing trends. Our findings have thus obvious implications for the use of volcanic gases as a monitoring tool of volcanic activity. They also illustrate the important role of sulphur. Magmas poor in sulphur or in which reactions involving sulphur are kinetically inhibited are likely to undergo a significant increase in their redox state during ascent.

The above results show that the redox state that magma records at depth does not necessarily mirror that of its escaping gases, in particular when they are released from levels shallower than the main reservoir. Thus, from a broader perspective, our findings have implications for our understanding on how past volcanic activity may have impacted Earth's atmosphere. Current models of the evolution of atmospheric oxygen implicitly assume that the redox state of magmatic rocks can be taken as equal to that of their outgassed products [START_REF] Holland | Volcanic gases, black smokers and the great oxidation event[END_REF][START_REF] Kasting | Mantle redox evolution and the oxidation state of the atmosphere[END_REF] . Our calculations show that, for silicic magmas, this assumption holds true only under a restricted set of conditions (for instance a magma starting at NNO-0.5 with 1 wt% gas, Fig. 2). Oxidized silicic magmas are particularly prone to redox change during ascent. A recent study [START_REF] Roeder | Morphology and composition of spinel in Pu'u'O'o lava (1996-1998), Kilauea volcano, Hawaii[END_REF] has stressed the difference in fO 2 retrieved from volcanic gas [START_REF] Moretti | On the oxidation state and volatile behavior in multicomponent gas-melt equilibria[END_REF] and volcanic glass [START_REF] Clemente | The solubility of sulphur in hydrous rhyolitic melts[END_REF] at Kilauea volcano: the basaltic glass, which is fully degassed, records an fO 2 lower by 1.2 log unit than the gases. Although our model is calibrated on Fe-poor liquids, its predictions are qualitatively in accord with such an observation which suggests that even for basaltic magmas in non-arc settings, redox change during degassing may occur [START_REF] Mathez | Influence of degassing on oxidation states of basaltic magmas[END_REF] . The corollary is that the iron redox state of a magma may significantly differ from its source, in contrast to conventional wisdom [START_REF] Carmichael | The redox states of basic and silicic magmas: a reflexion of their source regions[END_REF] . Altogether, this suggests that equating the redox state and composition of present day volcanic gases to those emitted in the geologic past [START_REF] Holland | Volcanic gases, black smokers and the great oxidation event[END_REF] may not be a correct assumption.

Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests and contributed equally to this work.

Correspondence and requests for materials should be addressed to A.B. (burgisse@cnrsorleans.fr) of fO 2 , each with three different gas contents (same initial conditions as in Fig. 2).

Natural range observed on active volcanoes in convergent settings [START_REF] Symonds | Volcanic-gas studies: methods, results, and applications[END_REF] , and range calculated by phase equilibra experiments 21 are also shown. B) Evolution of the gas composition for a magma oxidized at depth (∆NNO+1.5, 0.1 wt% gas). C) Evolution of the gas composition for a magma reduced at depth (∆NNO-0.5, 0.1 wt% gas).

Contents of O 2 are too low (<10 ppm) to be displayed

Supplementary Information

A coupled chemical-physical model of magma ascent allowed us to establish the relationship between the ascent of magma and its redox state. The chemistry formulation is based on that of Clemente et al. (2004), whereas the physical part is based on Burgisser and Gardner (2005). The gas phase is composed of n species, each with a molar fraction m i :

1 1 = ∑ = n i i m (1)
Conversion between molar fraction and weight fraction is:

∑ = j j j i i i M m M m x ( 2
)
where M i are molecular weights of each species. The total weight fraction of each species (w Ti ) is the sum of its exsolved part and its dissolved part:

( ) i b i i i gT Ti f a x w w + = (3)
where w gT is total gas weight fraction, a i and b i are solubility constants determined experimentally (Table 1), and f i is species fugacity. We assume the gas phase is an ideal mixture of non-ideal gases, which yields the following expression for fugacities:

P m f i i i γ = ( 4
)
where P is total pressure, and coefficients γ i are calculated at each pressure step following Shi and Saxena (1992). Using ( 2) and ( 4), (3) can be written as:

( i b i i i j j i i gT Ti P m a M m M m w w γ + = ) ∑ (5)

System H-O

The gas phase is composed of 3 species (H 2 O, H 2 , and O 2 ) that are in equilibrium at all times:

( )( ) ( ) 2 / 1 2 1 2 2 1 - - = fO fH O fH K (6)
where K 1 is an equilibrium constant calculated according to Robie et al. (1979).

Replacing the fugacities by their expressions in (4) yields:

( )( ) ( ) 2 / 1 2 2 1 2 2 2 2 1 - - = P m m m K O O H H O H O H γ γ γ (7)
Initially, w gT and fH 2 are set at a given pressure and temperature at depth. The molar composition of the gas phase can be determined by using (4) to calculate m H2 , (1) and (7) to calculate m O2 , and (1) to calculate m H2O . When decompression occurs, mass conservation requires that total amounts of oxygen and hydrogen of the system remain constant. Total weight percents of atomic oxygen (w TO ) and atomic hydrogen (w TH ) can be determined from initial gas composition by using:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + =
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = 2 2 2 2 2 2 H TH O H O TH H TH M W M w M w (8b)
and replacing the total amount of each species w Ti in ( 8) by their expressions in ( 5):

( ) ∑ ∑ + + = j j O gT O H b O H O H O H j j O H gT O TO M m m w M P m a M m m w M w O H 2 2 2 2 2 2 2 2 γ (9a) ( ) ( ) 2 2 2 2 2 2 2 2 2 2 2 2 2 H b H H H j j H gT O H b O H O H O H j j O H gT H TH M P m a M m m w M P m a M m m w M w H O H γ γ + + + = ∑ ∑ (9b)
When total pressure decreases, the molar fractions (m H2O , m H2 , and m O2 ) and the total amount of volatiles exsolved (w gT ) are constrained by mass conservation, chemical equilibrium, and solubility laws embedded in the four equations ( 1), ( 7), (9a), and (9b). Using ( 1) and ( 7), m H2O and m H2 can be expressed as a function of m O2 :

( )

1 2 2 2 1 2 2 2 1 1 - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - = P m K m m O O H O H O O H γ γ γ O H O H m m m 2 2 2 1 - - = (10)
Combining (9a) and (9b) as to eliminate w gT and using ( 10) yields an expression for m O2 that we solved with an iterative procedure to find its roots:

= + - - ) 2 ( ) 1 ( 2 2 2 2 O O H O TO O H TH m m M w m M w ( ) ( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - - + - 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 ) 1 ( 2 2 O O H O O H b O H O H O H O H b H H H m m m M P m a m M P m a O H H γ γ (11)
Because m H2O is on the order of 10 -1 and m O2 is usually <10 -10 , it is essential to carry all calculations with sufficient numerical precision. We used a precision of 20 digits, which ensured that errors on m O2 values were <0.1% when m O2 >10 -17 .

System H-O-S

The gas phase is composed of 6 species (H 2 O, H 2 , O 2 , S 2 , H 2 S and SO 2 ) that are in equilibrium at all times:

( )( ) ( ) 2 / 1 2 1 2 2 1 - - = fO fH O fH K ( )( ) ( ) 1 2 2 / 1 2 2 2 - - = fO fS fSO K ( )( ) ( ) ( ) 1 2 2 / 1 2 2 / 1 2 2 3 - - = O fH fS fO S fH K ( 12
)
where K 1 , K 2 , and K 3 are equilibrium constants (Ohomoto and Kerrick, 1977; Robie et al., 1979). The total weight percents of atomic oxygen (w TO ), atomic hydrogen (w TH ), and atomic sulfur (w TS ) are given by: Like the H-O system, the resolution is carried out in two steps. First, w gT , fH 2 O and fH 2 are set at depth so that the composition of the gas phase and the total amounts of each element could be determined. Second, each time pressure decreases, mass conservation (1) and ( 13), chemical equilibrium (12), and solubility laws (5) are used jointly to solve for the seven unknowns (molar fractions m i and the total gas weight fraction w gT ).

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + = 2 2 2 2 2 2 2 2 SO TSO O TO O H O TH O TO M w M w M w M w ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + = S H S TH H TH O H O TH H TH M w M w M w M w 2 2 2 2 2 2 2 2 2 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + =

Effects of Iron

When there is iron dissolved in the melt, the following redox reaction occurs:

3 2 2 2 1 O Fe O FeO = + (14)
It is currently not possible to include the redox effects of iron into our thermodynamic model because the activities of FeO and Fe 2 O 3 dissolved into a rhyolitic melt are not known. It is, however, possible to simulate such effects by using an empirical model that links, at constant temperature, the molar ratio of FeO and Fe 2 O 3 to the oxygen fugacity (Kress and Carmichael, 1991):

( ) T P h T P f m d c T b fO a m m KC KC k k k KC KC KC KC FeO O Fe 2 2 3 2 ln ln + + + + + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∑ (15)
where k is one of (FeO*, Al 2 O 3 , CaO, Na 2 O, and K 2 O), m k is the total molar fraction of these species, and a KC , b KC , c KC , d KC k , f KC , h KC are constants given in Kress and Carmichael (1991). We calculated the redox effect of iron in the H-O-S system by using (15) and keeping track of the total amounts of Fe (w TFe ), oxygen in the FeO and Fe 2 O 3 form (w TO(Fe) ), and oxygen in the O 2 , H 2 O, and SO 2 form (w TO ).

Because the relation (14) contains the total molar fraction of several oxides and our mass balance (13) is in weight fraction, conversions between the two units are needed, and the total composition of the melt should be known. Thus, we fixed the total weight fraction of 9 major oxides of a typical rhyolite (in wt%: SiO 2 =77.44, Al 2 O 3 =12.88, TiO 2 =0.07, CaO=0.6, Na 2 O=4.14, K 2 O=4.76, MgO=0.03, MnO=0.05, P 2 O 5 =0.03). These initial values add up to 100%, and need to be corrected to take in account the amounts of volatiles and iron. For each run, initial quantities of FeO* (w TFeO* ), w gT , fH 2 O and fH 2 are set at depth. Using the same procedure as described in "H-O-S system", total amounts of sulphur, hydrogen, and oxygen in the O 2 , H 2 O, and SO 2 form can be calculated. The initial amount of oxygen in the FeO and Fe 2 O 3 form is not known and should be determined using the empirical equation (14). Thus, the initial quantities of the 9 major oxides are normalized so that the 9 oxides, plus the volatiles species (w TO , w TS , w TH ), plus the total iron (w TFeO* ) add up to 100%. Then the total molar fraction of FeO*, Al 2 O 3 , CaO, Na 2 O, and K 2 O are calculated using:

∑ = j j Tj i Ti i M w M w m 1 (16)
The total amounts of iron and oxygen fixed by the iron can be expressed as:

) ( + + = (20)
The mass balance can now be rigorously carried out by adding the 9 oxides, the volatiles species (w TO , w TS , and w TH ), the total iron (w TFe ) and the oxygen fixed by the iron (w TO(Fe) ). The total molar fractions m k can be updated using (16) and are used in (15) for the rest of the simulation. Because the system is closed, all the quantities involved in the mass balance, except w TO and w TO(Fe) , are considered constant. The reaction (14) imposes an exchange of oxygen between that locked in the Fe-bearing species and that involved in the H-O-S redox reactions. The sum w TO + w TO(Fe) = w TO(Total) , however, remains constant. Thus, using (20), the total amount of "free" oxygen w TO can be related to fO 2 :

KC KC M m M m w w j j j O Fe Total TO TO 2 1 3 1 ) ( + + - = ∑ (21)
The numerical resolution of (21) needs to be coupled with that of the H-O-S system. For each pressure step, solving for the redox equilibrium of the volatile species (equations (1), ( 5), ( 12) and ( 13)) yields an fO 2 . This fugacity is used in (21) to obtain a new amount of "free" oxygen. The redox equilibrium of the volatile species is solved again with the updated w TO , and the process is conducted iteratively until desired precision (<10 -5 % on w TO ) is reached. An under-relaxation factor of 0.1 is applied to the updated w TO to ensure convergence in 20 to 30 iterations.

Runs with 1 wt% of total iron under typical storage conditions of arc rhyolites (fO 2 between NNO-0.5 and NNO+1.5) show that the buffering capacity of iron is limited to less than 0.2 unit log compared to the S-O-H system when little gas is present in the reservoir, and becomes insignificant when more gas is present in the reservoir (Fig. SI-1). Under such storage conditions, iron does not change the redox evolution of the magma with pressure. Similarly, the H 2 S/SO 2 ratio changes by less than 0.2 log units when little gas is present in the reservoir and is unaffected by iron when the reservoir hosts a large amount of gas. Stippled lines are runs with 1 wt.% total iron (FeO*) and solid lines are runs with no iron. Upper graph runs start with 0.1 wt.% gas at depth and lower graph runs start with 5 wt% gas at depth.

Physical part of the model

The conduit flow model is homogenous and one-dimensional. From mass and momentum conservation, the pressure evolution with depth (dP/dz) is given by: where v is magma velocity, ρ is magma density, r is conduit radius, and g is gravity.

The viscosity µ is calculated with the relationship by Hess and Dingwell (1996) for the melt and correcting for the effect of bubbles using the relationship by Dobran (1992). We note that using instead the relationship by Jaupart and Allègre (1991) for the effect of bubbles on viscosity did not change the trends caused by variations in volatile chemistry. Magma bulk density is:

l gT g gT w w ρ ρ ρ - + = 1 1 ( 23
)
where ρ g is gas density from perfect gas law, and ρ l is melt density. Replacing the derivative of ( 23) into ( 22) and rearranging gives: Above fragmentation (gas volume fraction >0.75), the viscosity is that of a dusty gas (Dobran, 1992), and degassing continues. Equation ( 24) is solved using a fourth-order Runge-Kutta algorithm with adjustable distance step. At each pressure step, the gas composition is calculated for the system considered (H-O or H-O-S). The terms dM/dz and dw gT /dz in (24) are evaluated implicitly, which can be justified by the fact that they are always several orders of magnitude smaller than inertia and viscous terms. Boundary condition is either atmospheric pressure at the vent, or initial velocity at depth. In the former case, a run starts with an initial guess for the velocity at depth, and the pressure is solved for until vent is reached. Initial velocity is then changed until upper boundary condition is satisfied. Because degassing occurs in equilibrium, changes in chemistry as a function of pressure are identical for the two cases.

Figure 1 :

 1 Figure 1: Fundamental relationship between magma ascent and magma redox state, for a

Figure 2 :

 2 Figure 2: Fundamental relationship between magma ascent and magma redox state, for a

Figure 3 :

 3 Figure 3: Evolution of the composition of an H-O-S gas phase during ascent of a rhyolite

  FIGURE 1

  ., Scaillet, B. Redox evolution of a degassing magma rising to the surface

=

  and the conversion(16) for w TFe and w TO(Fe) yields the total amount of oxygen fixed by the iron:

Figure SI- 1 :

 1 Figure SI-1: Effect of dissolved iron on the redox evolution of an ascending rhyolite.Stippled lines are runs with 1 wt.% total iron (FeO*) and solid lines are runs with no iron. Upper graph runs start with 0.1 wt.% gas at depth and lower graph runs start with 5 wt% gas at depth.

  M is the average molar mass of the gas phase:

Table 1 : Solubility constants

 1

	species	a i	b i
	H 2 O	1.063 10 -3	0.5399
	H 2	3.400 10 -7	1.2800
	SO 2	1.632 10 -10	1.3789
	H 2 S	8.239 10 -6	0.5145
	O 2	0	-
	S 2	0	-
	The a i and b i parameters were determined by fitting experimental
	solubility data of corresponding species to an empirical equation of
	the form w

i = a i f i bi . We have used data from ref .

12

for sulphur bearing species, ref. 13 for hydrogen and ref. 29 for water.

Burgisser, A., Scaillet, B. Redox evolution of a degassing magma rising to the surface

Acknowledgements

We thank M. Rutherford, P. Wallace, and an anonymous reviewer for critical comments that helped us to improve our model. A.B. acknowledges support from the Swiss National Science Foundation.

Supplementary Information is linked to the online version of the paper at www.nature.com/nature

Authors Information


Claude Le Bris

Pierre Rouchon
email: pierre.rouchon@mines-paristech.fr

Julien Roussel

Adaptive low

Keywords: numbers: 03.65.Yz, 02.60.Cb, 42.50.Pq, 42.50.Md, 02.70.Ss

and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Lindblad equations are notoriously challenging to simulate numerically. The two categories of approaches are deterministic approaches, on the one hand, and Monte-Carlo approaches [1], on the other hand. Both categories have their pros and cons. In the former category, the simulation is extremely effective when possible, but a major difficulty lies in the high dimensionality of the ambient space which drastically limits the applicability. In the latter category, dimensionality is not an issue, but the intrinsic noise of stochastic simulations affects the quality of the numerical results. The twofold purpose of this article is to present recent advances in either of the two categories of approaches.

In a previous work [2], the first two authors have presented a possible deterministic approach for the simulation of the Lindblad equation, actually borrowed from similar ideas introduced in [3,4] in the context of quantum filtering. The approach consists in approximating the evolution of the n × n density matrix ρ solution to the differential Lindblad equation using a reduced dynamics on the set of density matrices of some fixed rank m n. This reduced dynamics is obtained by taking the orthogonal projection of d dt ρ onto the tangent space to this set of rank-m matrices . The clear limitation of the approach lies in the fact that many practical problems are not reducible to a low-rank approximation, and further that, even when it is the case, the intrinsic dimensionality of the reduced dynamics is not necessarily known beforehand and may vary in time. So the question of adjusting on-the-fly the dimensionality m of the low rank dynamics immediately arises. As our first contribution in the present article, we describe below an adaptive low-rank simulation, the purpose of which is to significantly extend the applicability of the approach introduced earlier in [2]. The questions we examine in this work enjoy some similarity with questions arising in computational quantum chemistry, typically for multiconfiguration time-dependent Hartree and Hartree-Fock equations [START_REF] Lubich | From quantum to classical molecular dynamics: reduced models and numerical analysis[END_REF]6].

For problems definitely not amenable to deterministic simulation because of their prohibitively high dimensionality (which is indeed the case for many practically relevant problems), stochastic approaches are in order, see [1,[7][8][START_REF] Haroche | Exploring the Quantum: Atoms, Cavities and Photons[END_REF]. Although the low-rank dynamics no longer adequately represents the system, a reduced model, simulated deterministically, can however serve as a useful tool for the stochastic simulation of the high dimensional system. We employ the reduced dynamics as a control variate within a variance reduction method applied to the full, high dimensional stochastic system. Our second contribution is to demonstrate the efficiency of such a variance reduction method.

ADAPTIVE LOW-RANK APPROXIMATION

We consider throughout this article a Lindblad equation with, for simplicity (and this is by no means a limitation of our methods), a single decoherence operator L,

d dt ρ = -i[H, ρ] -1 2 (L † Lρ + ρL † L) + LρL † , (1)
where ρ is a n × n non-negative Hermitian matrix with Tr (ρ) = 1, H is a n × n Hermitian matrix and L is a n × n matrix. The reduced dynamics derived in [2] approximates, for n large, the above dynamics on the set of non-negative Hermitian matrices of rank m, m being an integer presumably much smaller than n. To make the approximation explicit, one introduces a system of two coupled differential equations for U and σ corresponding to the generic decomposition ρ LR = U σU † , where σ is a m × m strictly positive Hermitian matrix, U a n × m matrix with U † U = I m , and I m denotes the m × m identity matrix. That system reads as:

d dt U = -iHU + (I n -U U † ) -1 2 L † LU + LU σU † L † U σ -1 , (2)
d dt σ = -1 2 (U † L † LU σ + σU † L † LU ) + U † LU σU † L † U + 1 m Tr L † (I n -U U † )L U σU † I m . (3)
Notice that H only appears in (2) and not in (3), a fact that is particularly appropriate when H dominates L, in which case (3) may be understood as a slow evolution as compared to the dynamics (2). Then the projection of the original Lindblad dynamics (1) onto the tangent space to the set D m of density matrices of rank m takes the explicit form

d dt ρ LR = -i[H, ρ LR ] -1 2 (L † Lρ LR + ρ LR L † L) + Lρ LR L † -(I n -P ρLR )Lρ LR L † (I n -P ρLR ) + Tr(LρLRL † (In-Pρ LR )) m P ρLR , (4)
where the orthogonal projection on the image of ρ LR , P ρLR = U U † , only depends on ρ LR . Notice that the rightmost term of (4) allows Tr (ρ LR ) to be preserved in time.

In the sequel, we denote by L the right-hand side of (1), L that of (4), and by L ⊥ = L -L . We have described in details in [2] how system (2)-(3) may be efficiently simulated and then provides an accurate approximation of (1) in the case when the rank can be actually reduced. In that work, the rank m was prescribed beforehand. Our purpose here is to explain how the approach can be amended so as to allow for a dynamical adaptation of the rank m of the reduced system.

The adaptation we suggest is based on the evaluation, and update, of the projection error

L ⊥ (ρ LR ) = (I n -P ρLR )Lρ LR L † (I n -P ρLR ) - Tr Lρ LR L † (I n -P ρLR ) m P ρLR
committed when replacing (1) by ( 4). This error may be reduced upon adding one dimension to the m-dimensional subspace Im(U ) associated to the projector P ρLR . The best possible such dimension to add is that for which the projection error is minimal. The rank-m projector

P ρLR = U U † ∈ R m×m is modified into the rank-m+1 projector P ρLR +Q = U U † +V V † ∈ R m×m where Q = V V † , with V ∈ R n ,
is the one-dimensional projector associated to the one dimension added. Denoting by

G = (I n -P ρLR )Lρ LR L † (I n -P ρLR ), a straightforward calculation yields Tr 2 L ⊥ P +Q (ρ LR ) = (I n -Q)G(I n -Q) 2 + 1 m+1 Tr 2 ((I n -Q)G) .
One may then prove (see the details in [START_REF] Roussel | Numerical simulation of high dimensional open quantum systems[END_REF]) that the directions

V = Argmin V ⊥U, V =1 L ⊥ Pρ LR +Q (ρ LR ) (5)
minimizing the projection error are the eigenvectors of the symmetric matrix G associated to its largest eigenvalue. The matrix G being of large size n, determining its largest eigenvector is challenging computationally. To this end, we notice that

V † GV = V † (I n -P ρLR )LU √ σ 2 ,
and that the range of (I n -P ρLR )LU is of dimension less or equal to m and is orthogonal to the range of ρ LR . Thus, denoting by (φ j ) 0≤j<r an orthonormal basis of the range of (I n -P ρLR )LU with r ≤ m, it is sufficient to consider V as linear combination of the φ j to get an eigenvector of G with largest eigenvalue: V = r j=1 v j φ j where the rdimensional vector v of component v j corresponds to the eigenvector with largest eigenvalue of

K = Φ † (I n -P ρLR )Lρ LR L † (I n -P ρLR )Φ,
with Φ the n × r matrix formed by the r vectors φ j . Since K is of size r ≤ m, this provides an effective manner to determine the optimum in [START_REF] Lubich | From quantum to classical molecular dynamics: reduced models and numerical analysis[END_REF].

We have performed a comprehensive series of test of the approach. The practical implementation of the dynamical adaptivity of the rank is performed as follows.

We denote by θ(ρ

LR ) ≡ L ⊥ (ρ LR ) L (ρ LR )
, fix a maximal angular error θ max and update the rank

• increasing m by 1, when θ(ρ LR ) > θ max and then complementing U via the solution of ( 5);

• reducing m by 1, when the smallest eigenvalue λ min of σ is such that θ(ρ LR ) + λ min < 1 2 θ max .

As an illustrative example, we consider a qubit resonantly coupled to a quantized harmonic damped oscillator:

d dt ρ = Ω 0 2 [a † σ --aσ + , ρ] -κ(nρ/2 + ρn/2 -aρa † ), (6)
with Ω 0 > 0 the vacuum Rabi pulsation, a the photon annihilation operator, σ -the qubit lowering operator (σ -= |g e|), σ + = σ - † , n = a † a the photon-number operator and 1/κ > 0 the oscillator damping time. In Figure 3.20, page 156 of [START_REF] Haroche | Exploring the Quantum: Atoms, Cavities and Photons[END_REF], numerical simulations of quantum collapse and revivals are presented, for κ = 0, when the qubit is initially in the excited state |e and the oscillator in a coherent state with n = 15 photons. We consider here the same system with κ = Ω 0 /500, a value small compared to Ω 0 . This corresponds to a photon life-time 1/κ around 10 times as large as the revival time T r = 4π √ n Ω0 . In this simple case, we can perform the full-rank simulation with high precision. This simulation illustrated on Figure (1) provides us with a reference calculation. As shown in [2], a constant rank of 4 is sufficient to compute accurately the solution t → ρ t for t between 0 and 2T r . For intermediate values of t, larger than 2T r but not excessively larger, the quantum state ρ gets more and more mixed. For t very large, ρ becomes again pure, since its limit for t → +∞ is the lowest energy state (qubit in the ground state |g , oscillator with zero photon). This behavior, illustrated by the numerical simulations of Figure 2, is well captured by our adaptive approach.

DENOISED MONTE-CARLO APPROACH

In the case when the Lindblad equation ( 1) is genuinely high-dimensional, the model reduction previously described is likely to be either ineffective or inaccurate, while the direct integration of the equation is out of reach. In this situation, the classical approach is to use a Monte-Carlo sampling. One derives a stochastic dynamics on the wave function |ψ t such that the density ρ(t) = E (|ψ t ψ t |) constructed from |ψ t (E stands for expectation value, i.e. ensemble average) solves the Lindblad equation (1). In practice, M independent trajectories, ψ (k) t with k = 1, . . . , M , are then simulated using that stochastic dynamics and the estimator of the mean

ρMC (t) = 1 M M k=1 ψ (k) t ψ (k) t (7)
is used as an approximation of ρ(t). The practical difficulty of Monte-Carlo approaches is, as briefly mentioned, above, the variance, that is, the noise intrinsically present in the approach.

In principle, there are infinitely many dynamics on |ψ t that are consistent with the Lindblad dynamics. Interestingly, a straightforward calculation from (7) shows that

E Tr (ρ MC -ρ) 2 = 1 -Tr ρ 2 M . (8)
The variance of the estimator of the mean ρMC is therefore independent of the specific unravelling choice, namely the stochastic dynamics set on |ψ t , provided that (color online) High precision full-rank numerical simulation of t → ρ(t) for the quantum collapse and revivals of a qubit resonantly coupled to a slightly damped quantum harmonic oscillator governed by (6). Top plot: the solid curve corresponds to the excited population e|ρ|e versus time; the normalized time corresponds to t/Tr where Tr is the revival time; the oscillator damping time is around 10Tr; the initial state ρ0 corresponds to the qubit in excited state |e and oscillator in a coherent state with n = 15 photons; a truncation of the number of photons to a maximum of 2n = 30 yields an underlying Hilbert space of dimension 62. Bottom plot: zoom of top plot for a normalized time between 0 and 10. state |ψ t remains normalized. This is easily seen in the proof of ( 8). This property is, of course, a remarkable peculiarity of the present context. And the discretization in time of the process can of course slightly affect that property. The classical unravelling choice (see [START_REF] Gisin | [END_REF][START_REF] Breuer | The Theory of Open Quantum Systems[END_REF][START_REF] Barchielli | Quantum Trajectories and Measurements in Continuous Time: the Diffusive Case[END_REF]) is the Wiener process defined by

d |ψ t = D 1 (|ψ t ) dt + D 2 (|ψ t )dW t , (9)
with the drift term

D 1 (|ψ ) = -iH |ψ + D 0 1 (|ψ ) where D 0 1 (|ψ ) = 1 2 L + L † |ψ L -L † L - 1 4 L + L † 2 |ψ |ψ , (10)
and the diffusion We have used the notation A |ψ = ψ|A|ψ . We emphasize that other choices of dynamics on |ψ , all consistent with the Lindblad equation through ρ(t) = E (|ψ t ψ t |), could be made. In particular, Poisson processes could be considered instead of Wiener processes (a choice that can be seen as more natural given the applications addressed in the present work where photons are emitted or absorbed). In any event, given the property (8) and assuming that the Poisson process remains normalized, the variance remains identical. We indeed double-checked that, in actuality, using Poisson processes does not bring further practical variance reduction, see [START_REF] Roussel | Numerical simulation of high dimensional open quantum systems[END_REF] for more details.

D 2 (|ψ ) = L - 1 2 L + L † |ψ |ψ . (11)
Variance is an issue for the numerical simulation and affects the accuracy of the results. It is thus desirable to come up with further, dedicated variance reduction approaches that may reduce the computational cost at accuracy fixed, or improve the accuracy for a given com-putational cost. An approach that has proved effective in many engineering sciences for reducing variance of Monte-Carlo simulations is that of control variate, see e.g. [14, page 54]. In short, the approach consists in concurrently simulating the original system under consideration and a system correlated to that original system so as to minimize the variance in the simulation of the former system. Intuitively, the approach works by "cancellation" of the noise because the same random draws are used for both systems. More specifically, we consider here as control variate the low-rank approximation ρ LR of previous section. Even though that low-rank system is not a correct approximation of the original system (which we have deliberately assumed here high-dimensional), it is sufficiently correlated to that system to provide an efficient variance reduction. Practically, we construct an ("Control-Variate") estimated density depending on the adjustable scalar parameter λ,

ρCV = ρMC + λ(ρ LR -ρMCLR ), (12)
as a combination of the original Monte-Carlo estimated density ρMC (constructed from the simulation of ( 7)-( 9)), the estimated density ρMCLR constructed from a low-rank dynamics ("Monte-Carlo-Low-Rank", see below ( 13)), and the density ρ LR obtained upon solving the corresponding low-rank Lindblad equation. Since by construction

E[ρ MCLR (t)] = ρ LR (t),
the approximation method ( 12) is unbiased (taking the expectation of both sides of [START_REF] Breuer | The Theory of Open Quantum Systems[END_REF]

yields ρ = E[ρ CV ] = E[ρ M C ]).
The scalar parameter λ is adjusted so as to minimize the variance E Tr (ρ CV -ρ)

2

, a polynomial in λ of degree 2, in such a way that

E Tr (ρ CV -ρ) 2 E Tr (ρ MC -ρ) 2
so that the simulation of ρCV is eventually more effective than that of ρMC . The more correlated the reduced model and the original model, the closer λ to one, the smaller this variance and thus the more efficient the denoising. In passing, we notice that, although this will not be the case in the actual numerical experiments we perform, the low-rank dynamics ρ LR (t) could itself be chosen with an adaptive rank, as in the previous section. Likewise, we could pick as control variate another dynamics than the low-rank dynamics, if a more convenient one is available.

The remaining question is to derive a stochastic dynamics on a wave function |ψ LR such that E (|ψ LR ψ LR |) = ρ LR , given that the dynamics of ρ LR is known since easy to compute via (U, σ) solutions of ( 2) and (3). For this purpose, it is a natural idea to seek |ψ LR under the form |ψ LR = U |ν where the reduced wave function |ν ∈ R m is a stochastic process to be determined. It can be shown (see [START_REF] Roussel | Numerical simulation of high dimensional open quantum systems[END_REF]) that the correct dynamics to consider reads

d |ψ LR = P ρLR D 1 (|ψ LR ) dt + P ρLR D 2 (|ψ LR )dW t + (I n -P ρLR ) - 1 2 L † L + LU σU † L † U σ -1 U † |ψ LR dt + Tr (I n -P ρLR )Lρ LR L † 2m U σ -1 U † |ψ LR dt, (13)
with D 1 and D 2 defined in [START_REF] Roussel | Numerical simulation of high dimensional open quantum systems[END_REF] and [START_REF] Gisin | [END_REF].

Since ρMCLR = 1 M M k=1 ψ (k) LR ψ (k)
LR , a simple computation, exploiting the fact that both ψ (k) and ψ

(k) LR are normalized , shows that M E Tr (ρ CV -ρ) 2 = 1 -Tr ρ 2 LR λ 2 + 2 E | ψ|ψ LR | 2 -Tr (ρρ LR ) λ + 1 -Tr ρ 2 .
Thus the optimal denoising choice for λ reads

λ = Tr (ρρ LR ) -E | ψ|ψ LR | 2 1 -Tr (ρ LR 2 ) . (14)
In practice, and also in the simulations of Figure (3), the adjustable parameter λ is given by ( 14) where the ρ is replaced by ρMC (we do not have access to ρ itself) and where E | ψ|ψ LR | 2 is replaced by the estimated mean

1 M M k=1 ψ (k) |ψ (k) LR 2
.

Figure (3) illustrates the interest of this denoising method for the Lindblad system simulated in Figure (2). We take M = 400 trajectories. The dynamics of |ψ LR is based on a low-rank approximation ρ LR of constant rank m = 2. On the bottom plot of Figure 2, we observe that, at normalized time 2, an accurate approximation of ρ must be of rank 10 or larger. Nevertheless, Figure 3 indicates that, at normalized time 2, a reduction of the standard deviation around 50% is obtained with ρCV instead of ρMC . Such a variance reduction corresponds, with a classical quantum Monte-Carlo method, to an increase of the number of trajectories by a factor 4. Such a gain confirms the definite interest of combining deterministic low-rank approximations with quantum Monte-Carlo trajectories for the numerical simulation of high-dimensional Lindblad equation.

ACKNOWLEDGMENTS

The first two authors are partially supported by the ANR, Projet Blanc EMAQS ANR-2011-BS01-017-01. The solid black curve corresponds to Tr ((ρ -ρCV) 2 ), the dashed blue one to Tr ((ρ -ρMC) 2 ) and the dotted red one to λ. The number of trajectories M is set to 400. The lowrank approximation ρLR used for the control variate governed by ( 13) is maintained at the constant rank m = 2.

  FIG. 1. (color online)High precision full-rank numerical simulation of t → ρ(t) for the quantum collapse and revivals of a qubit resonantly coupled to a slightly damped quantum harmonic oscillator governed by(6). Top plot: the solid curve corresponds to the excited population e|ρ|e versus time; the normalized time corresponds to t/Tr where Tr is the revival time; the oscillator damping time is around 10Tr; the initial state ρ0 corresponds to the qubit in excited state |e and oscillator in a coherent state with n = 15 photons; a truncation of the number of photons to a maximum of 2n = 30 yields an underlying Hilbert space of dimension 62. Bottom plot: zoom of top plot for a normalized time between 0 and 10.

FIG. 2 .

 2 FIG. 2. (color online)Comparison between the full-rank trajectory t → ρ(t) illustrated on figure1with the adaptive rank trajectory t → ρLR(t) governed by(2)(3). Top plot: the solid curve corresponds to Tr ((ρ -ρLR) 2 ) versus time with scale on the left. The dashed curve corresponds to the adaptive rank m with scale on the right; the initial low rank state ρLR coincides with ρ0; the initial rank m is thus set to one; it evolves according to the maximal angular error θmax = 1 1000 . Bottom plot: zoom of top plot for a normalized time between 0 and 10.

FIG. 3 .

 3 FIG. 3. (color online)Variance reduction corresponding to the quantum collapse and revivals in the bottom plot of Figure2. The solid black curve corresponds to Tr ((ρ -ρCV) 2 ), the dashed blue one to Tr ((ρ -ρMC) 2 ) and the dotted red one to λ. The number of trajectories M is set to 400. The lowrank approximation ρLR used for the control variate governed by (13) is maintained at the constant rank m = 2.