Aurélien Grabsch

Alexis Poncet

Pierre Rizkallah

Pierre Illien

Olivier Bénichou

Exact closure and solution for spatial correlations in single-file diffusion

INTRODUCTION

Single-file transport, where particles diffuse in narrow channels with the constraint that they cannot bypass each other, is a fundamental model [1,2] for tracer subdiffusion in confined systems. The very fact that the initial order is maintained at all times leads to the subdiffusive behavior ⟨X 2 t ⟩ ∝ √ t of the position of a tracer particle (TP) [3], in constrast with the regular diffusion scaling ⟨X 2 t ⟩ ∝ t. This theoretical prediction has been experimentally observed by microrheology in zeolites, transport of confined colloidal particles, or dipolar spheres in circular channels [4][5][6].

The Symmetric Exclusion Process (SEP) is an essential model of single-file diffusion. Particles, present at a density ρ, perform symmetric continuous-time random walks on a one-dimensional infinite lattice with unit jump rate, and with the hard-core constraint that there is at most one particle per site (Fig. 1(a)). The SEP has become a paradigmatic model of statistical physics and it has generated a huge number of works in the mathematical and physical literature (see, e.g., Refs. [1,2,7,8]). A major recent advance has been achieved with the calculation of the large deviation function of the position X t of a tracer in the long time limit [9,10]. It gives access to all the long-time cumulants of X t , which are in particular found to behave anomalously as √ t [9,11]. Similarly, the cumulants of the time integrated current through the origin Q t have been shown to also scale as √ t, and the large deviation function has been determined [8].

This collection of anomalous behaviors in the SEP originates from the strong spatial correlations in the singlefile geometry, which makes them determining quantities. Even if this has been recognized qualitatively for long and that the case of dense and dilute limits have been recently studied [12], up to now there is no quantitative determination of the bath-tracer correlations at arbitrary density, despite extensive effort. Indeed, although the SEP has been studied for more than 40 years, analytic formulas for these functions are still missing. The calculation of these correlations in the SEP actually constitutes an open many-body problem, which we solve here. More generally, we put forward bath-tracer correlations as fundamental quantities to analyze single-file diffusion, since we show that they satisfy a strikingly simple exact closed equation. The central role of this key equation as a novel tool for interacting particle systems is further demonstrated by showing that it applies to out-of equilibrium situations, other observables and other representative single-file systems.

We consider a SEP of average density ρ, with a tracer, of position X t at time t, initially at the origin. The bath particles are described by the set of occupation numbers η r (t) of each site r ∈ Z of the line at time t, with η r (t) = 1 if the site is occupied and η r (t) = 0 otherwise (see Fig. 1(a)). The statistics of the position of the tracer is described by the cumulant-generating function, whose expansion defines the cumulants κ n of the position:

ψ(λ, t) ≡ ln e λXt ≡ ∞ n=1 λ n n! κ n (t). (1)
Its evolution equation is given by (see Eq. (S20) in Supplementary Information (SI)):

dψ dt = 1
2 (e λ -1)(1 -w 1 ) + (e -λ -1)(1 -w -1 ) , (2) where the generalized density profiles (GDP) generating function is defined by w r (λ, t) ≡ ⟨η Xt+r e λXt ⟩/⟨e λXt ⟩.

Note that, besides controlling the time evolution of the cumulant-generating function, w r (together with ψ) com-

FIG. 1. Models and bath-tracer correlations. (a)

The simple exclusion process (SEP). The tracer is located at position Xt, and Qt is the integrated current of particles through the origin up to time t (i.e. the total flux of particles between sites 0 and 1). The occupation numbers of the sites are denoted ηr. (b) In this article, we put forward bath-tracer correlations ⟨ηX t +r e λX t ⟩/⟨e λX t ⟩ (or bath-current correlations in the case of the integrated current Qt) as fundamental quantities to analyze single-file diffusion, as they satisfy the simple closed integral equation (6). Besides being key technical tools, α) they characterize the bath-tracer coupling (see Eq. ( 4)), β) they quantify the response of the bath of particles to the perturbation induced by the tracer (see item (iv) after Eq. ( 13)) and γ) in turn they control the large deviations of the subdiffusive motion of the tracer (see Eqs. (1), ( 2) and ( 13)). Other representative models of single-file systems: (c) the Random Average Process: at exponential times, particles on a line can jump in either direction to a random fraction of the distance to the next particle; (d) the Kipnis-Marchioro-Presutti model: each site carries a continuous energy variable. At exponential times, the total energy of two neighbouring sites is randomly distributed between them.

pletely characterizes the joint cumulant-generating function of (η Xt+r , X t ) and thus the bath-tracer correlations [12]. The GDP-generating function is therefore a key quantity, and the next step consists in writing its evolution equation from the master equation describing the system. However, similarly to Eq. ( 2), it involves higher-order correlation functions. In fact, we are facing an infinite hierarchy of evolution equations, which is the rule for tracer diffusion (and for other observables such as the integrated current through the origin) in interacting particle systems [9,[13][14][15], and whose closure has remained elusive up to now. We provide below a closed equation which allows the determination of the GDP-generating function in the hydrodynamic limit (large time, large distance).

In this limit, the position of the tracer satisfies a large deviation principle [9,10,16], which implies that the cumulant-generating function scales as ψ ∼ ψ√ 2t. In fact, this anomalous behavior originates from the more general scaling form

w r (λ, t) ∼ t→∞ Φ λ, v = r √ 2t ≡ ∞ n=0 λ n n! Φ n (v) (4)
of the GDP-generating function, where the coefficient Φ n gives the large-scale limit of the joint cumulant ⟨η Xt+r X n t ⟩ c of the tracer's position X t and the occupation number η Xt+r measured in its frame of reference (Fig. 1(b)). In the following, we will drop the argument λ of Φ for convenience.

RESULTS

We report here (see Materials and Methods and Section II.A of SI for details) that the two functions (rescaled derivatives of the profiles)

Ω ± (v) ≡ ∓2 ψ Φ ′ (v) Φ ′ (0 ± ) defined for v ≷ 0 (5)
are entirely determined by the closed Wiener-Hopf integral equations [17] with a Gaussian kernel

Ω ± (v) = ∓ω e -(v+ξ) 2 +ξ 2 -ω R ∓ Ω ± (z)e -(v-z+ξ) 2 +ξ 2 dz, (6)
where ξ ≡ d ψ dλ and we have analytically continued Ω + to v < 0 and Ω -to v > 0. The parameter ω is determined by the boundary conditions (see Eq. ( 5))

Ω + (0) = -Ω -(0) = -2 ψ, (7)
so that the functions Ω ± (v) are parametrized by ψ. At this stage, the expression of ψ(λ) has not been determined yet, but it can be obtained in the following way. First, Φ is deduced by integration of Ω ± , with

Φ(±∞) = ρ, ( 8
)
by definition, and the boundary conditions

Φ ′ (0 ± ) ± 2 ψ e ±λ -1 Φ(0 ± ) = 0. ( 9
)
The resulting Φ are at this stage parametrized by ψ and λ. Then, by using the large time limit of Eq. ( 2),

1 -Φ(0 -) 1 -Φ(0 + ) = e λ , (10)
ψ can be written as a function of λ, and we finally obtain the desired GDP-generating function Φ(λ, v).

DISCUSSION

Several comments are in order. (i) We show in SI (Section I.E) that the boundary condition ( 9) is exact; furthermore, we argue below that the bulk equation ( 6), and thus the obtained GDP-generating function Φ, are also exact. (ii) Importantly, the Wiener-Hopf equations ( 6) can be solved explicitly in terms of the one-sided Fourier transforms [17]:

R ± Ω ± (v)e ikv dv = ± (1 -exp[-Z ± ]) , (11)
where

Z ± ≡ 1 2 n≥1 (-ω √ π e -1 4 (k+2iξ) 2 ) n n erfc ± √ n ξ - ik 2 .
(12) (iii) As a byproduct, our approach yields the cumulant generating function ψ (or equivalently the large deviation function of the tracer's position),

ψ = - 1 2 √ π Li 3 2 (- √ πω), Li ν (x) ≡ n≥1 x n n ν , ( 13
)
which is shown in SI (Section II.E) to be identical to the exact expression obtained using the arsenal of integrable probabilities in [9,10]. (iv) Additionally, we obtain a full characterization of the spatial bath-tracer correlations and in particular analytical expressions of the Φ n by using the procedure described above (see Section II.C of SI for explicit expressions, which extend at arbitrary density the expressions given in [12] in the dilute (ρ → 0) and dense (ρ → 1) limits, and Fig. 2(a) and (b) for comparison with numerical simulations). (v) Our approach also provides the conditional profiles ⟨η Xt+r |X t = x⟩ defined as the average of the occupation of the site X t + r given that the tracer is at position x. Indeed, in the hydrodynamic limit, ⟨η

Xt+r |X t = x⟩ ∼ t→∞ Φ(λ * , v) where v = r/ √ 2t, λ * is defined by ξ(λ * ) = d ψ dλ (λ * ) = x/
√ 2t and Φ is the GDP-generating function determined above. While the (unconditional) profiles ⟨η Xt+r ⟩ are flat, the conditional profiles allow to probe the response of the bath of particles to the perturbation created by the displacement ξ of the tracer: in particular, for ξ > 0, it leads to an accumulation of bath particles in front of the tracer and a depletion behind (see Fig. 2(c)), quantified by the simple conservation relation

+∞ 0 (Φ(v) -ρ)dv - 0 -∞ (Φ(v) -ρ)dv = ρξ, (14)
which is a consequence of (11,12). Another striking feature is the non-monotony of the conditional profile Φ(λ * , 0 + ) in front of the tracer as a function of the rescaled position of the tracer ξ(λ * ), see Fig. 2(d). Surprisingly, Φ(λ * , 0 + ) does not saturate to 1 as ξ → +∞, but instead returns to its unperturbed value ρ, which results from a global displacement of bath particles induced by the tracer. Importantly, Equation (6) describes several other situations of physical relevance. (i) First, it applies to the out-of-equilibrium situation of an initial step of density ρ + for x > 0 and ρ -for x < 0, with the tracer initially at the origin. This paradigmatic setup has attracted a lot of attention [8-10, 18, 19] since it remains transient at all times and never reaches a stationary state. The GDP-generating function Φ is obtained from the solution (11) by following the procedure described above, upon only changing the boundary condition (8) into Φ(±∞) = ρ ± . Again, we recover the results of [9,10] on the cumulant generating function ψ. Additionally, we obtain the complete spatial structure of the bath-tracer correlations (see Fig. 3(a) and Section II.D.1 of SI for explicit expressions). (ii) Second, and strikingly, it also gives access to the statistics of other observables, as exemplified by the case of the integrated current through the origin Q t (see Section II.D.3 of SI for the application to the generalized current, which is an extra observable), defined as the total flux of particles between sites 0 and 1 during a time t. This quantity has been the focus of many studies, both in the context of statistical physics [8,13,18,20,21] and mesoscopic transport [22][23][24], in particular in the nonequilibrium situation ρ -̸ = ρ + [8,18]. Note that while the statistics of tracer diffusion and integrated current are easily related in the case of quenched initial conditions [25], the relation is more entangled for the annealed case considered here due to the fluctuations of the initial condition. The quantities introduced previously (1-3) on the example of tracer diffusion are naturally adapted by substituting Q t for X t . The corresponding profiles Φ Q are then obtained as a particular case of Equation ( 6) by setting ξ = 0, completed by modified boundary conditions (9,10) derived from the microscopic model (see Section II.D.2 of SI). In particular, the resulting Eq. ( 13) gives back the exact cumulant generating function of Q t obtained in [8] by Bethe ansatz, since in this case we find that

ω Q √ π = ρ -(1 -ρ + )(e λ -1) + ρ + (1 -ρ -)(e -λ -1), ( 15
)
which coincides with the single parameter involved in [8,18]. Additionally, the Φ Q determined here provides the associated spatial structure (see Fig. 3(b) and Section II.D.2 of SI for explicit expressions). These profiles have been introduced and studied numerically in [26] for an infinite system (see also [14] for a finite system between two reservoirs), but no analytical expressions were available until now. (iii) Finally, beyond the SEP, it applies to other representative single-file systems of interacting particles with average density ρ (see also [12], which is however limited to the calculation of the first order Φ 1 ). Such systems can be described at large scale by two quantities: the diffusivity D(ρ) and the mobility σ(ρ) [27]. The case of the SEP considered above corresponds to D(ρ) = 1/2 and σ(ρ) = ρ(1 -ρ). Equation ( 6) (with adaptations of equations (9,10) given in Section III.C of SI) more generally applies to single-file systems with D(ρ) = 1/2, σ ′′ (ρ) constant and σ(0) = 0, by replacing the -ω that multiplies the integral in (6) by ωσ ′′ (ρ)/2. Important cases covered by our approach include (see Fig. 1 for definitions and Sections III.D and III.E of SI for explicit expressions): (a) the model of hard Brownian particles (σ(ρ) = ρ) for which the GDP-generating function of [12] is recovered; (b) the Kipnis-Marchioro-Presutti (KMP) model [28] (σ(ρ) = ρ 2 , see Fig. 3(c)), which describes situations as varied as force fluctuations in packs of granular beads [29], the formation of clouds and gels, self-assembly of molecules in organic and inorganic materials and distribution of wealth in a society (see [30] and references therein). (c) the Random Average Process (RAP), which appears in a variety of problems such as force propagation in granular media, models of mass transport or models of voting systems [29,[31][32][33]. Although σ ′′ (ρ) is not constant in this case, the GDP-generating function can be deduced from our results thanks to a mapping between tracer diffusion in the RAP and the integrated current in the KMP model [34] (see Fig. 3(d)).

We finally argue that the central equation ( 6) is exact for the following reasons. (i) We show in SI (Section IV) that the Macroscopic Fluctuations Theory (MFT) [35] can be used to determine perturbatively the first coefficients Φ n analytically. These coefficients computed up to order n = 5 (which is the highest order for which we managed to determine the integrals involved) coincide with those obtained by our approach. Furthermore, the agreement holds also non-perturbatively in λ, as displayed in Fig. 2(c) and Fig. 3(c4),(d4) where the numerical solution of the MFT equations is compared to the analytical solution (11,12). Moreover, and as mentioned above, the exact expression of the cumulant generating functions of (ii) the tracer position of [9] and (iii) the integrated current of [8] are contained in our approach, including the case of an initial step of density.

All together, we have determined analytically the spatial correlations in the SEP, which allowed us to fully quantify the response of the bath to the perturbation induced by a tracer. Besides being paramount physical observables, these correlations have been shown to be fundamental technical quantities, since they satisfy a strikingly simple closed equation and control large deviations in single-file diffusion. This very same equation applies to a variety of situations involving single-file transport, which makes it a novel and promising tool to tackle interacting particle systems.

MATERIALS AND METHODS

Analytical calculations for the SEP

Details on analytical calculations are provided in SI.

We sketch here the main steps that led to the closed equation ( 6) for the SEP. The starting point is a master equation describing the time evolution of the complete system (bath and tracer in the SEP), from which we obtain the time evolution of the cumulant generating function ψ and the GDP-generating function (3). The main difficulty is that the latter involves higher-order correlation functions.

The next step consists in using the scaling (4) of the GDP-generating function and ψ(λ, t) ∼ ψ(λ) √ 2t to derive the hydrodynamic limit of the problem (details given in Section I.E of SI). The obtained bulk equation, valid at arbitrary density, is still not closed. We explain in SI that a closed equation obeyed by Φ has to satisfy the following constraints: (i) it must reduce to the known equations obtained in the limits of high and low density in [12]; (ii) it should also reproduce, as a byproduct, the cumulants of the tracer's position derived recently in [9,10]; (iii) additional constraints concern the way the different parameters appear in the equation (see Section II.A of SI for details); (iv) finally, the equation we write should have a "proper scaling" with time determined in Section I.E of SI.

Following these ideas and constraints, we obtain a first closed equation which holds at lowest orders in λ, see (S40) in SI, which properly reproduces the known cumulant κ n for n ≤ 6. This equation is conveniently rewritten by introducing the new functions Ω ± (v) defined by Eq. ( 5). Extension of this equation to arbitrary order in λ and then its resummation yields the closed equations (S48,S49) of SI. Several technical steps detailed in Section II.A of SI allow us to finally transform them into the closed Wiener-Hopf integral equation (6), which is our central result.

Numerical simulations for the SEP

The numerical simulations of the SEP are performed on a periodic ring of size N , with M = ρN particles at average density ρ. The particles are initially placed uniformly at random. The jumps of the particles are implemented as follow: one picks a particle uniformly at random, along with one direction (left and right with equal probabilities). If the chosen particle has no neighbor in that direction, the jump is performed, otherwise it is rejected. In both cases, the time of the simulation is incremented by a random number picked from an exponential distribution of rate N . We keep track of one particle (the tracer) and compute the moments of its displacement and the generalized density profiles. The averaging is performed over 10 8 simulation.

Extensions

Analytical (Section III of SI) and numerical (Section V of SI) extensions (other systems than the SEP, other

I. GENERAL EQUATIONS AND HYDRODYNAMIC LIMIT

In this Section, we mostly recall the equations and results obtained in [1] in order to provide a self-contained document. We additionally introduce in Section I C the conditional profiles discussed in the main text after Eq. (13).

A. Master equation of the SEP

We consider the symmetric exclusion process (SEP) with a tracer at position X. We denote the configuration of the system by η = {η r } r∈Z with η r ∈ {0, 1} the occupation number of site r by the bath particles (η r = 1 if site r is occupied, 0 otherwise). At time t, the system is characterized by a probability law P (X, η, t).

We initially start from the equilibrium distribution of the occupations, and the tracer at the origin (with the convention that the site occupied by the tracer is empty of bath particles):

P (X, η, 0) = δ X,0 δ η0,0 r∈Z * δ ηr,γr , (S1)
where γ r are independent Bernouilli variables with parameter ρ (density of the system) and the site 0 is treated independently because it is occupied by the tracer, and not a bath particle.

The probability P (X, η, t) obeys the following master equation,

∂ t P (X, η, t) = 1 2 r̸ =X,X-1 P (X, η r,+ , t) -P (X, η, t) + 1 2 µ=±1 (1 -η X )P (X -µ, η, t) -(1 -η X+µ )P (X, η, t) , ( S2
)
where η r,+ is the configuration η in which the occupations of sites r and r + 1 are exchanged. The first term corresponds to the jumps of the bath particles while the second one takes into account the displacement of the tracer.

B. Observables and large-times scalings

We consider the cumulant-generating function of the displacement of the tracer,

ψ(λ, t) ≡ ln e λXt . (S3) At large time t, it scales as √ t [2, 3], ψ(λ, t) ∼ t→∞ ψ(λ) √ 2t , ψ(λ) = n≥1 κ n λ n n! , ( S4
)
where κ n is the n th cumulant of the tracer's position (rescaled by √ 2t). We also consider the generalized profiles,

w r (λ, t) ≡ ⟨η Xt+r e λXt ⟩ ⟨e λXt ⟩ . (S5)
Their expansion in powers of λ gives the cross-cumulants ⟨η Xt+r X n t ⟩ c between the occupations and the position of the tracer. At large time, they satisfy a diffusive scaling,

w r (λ, t) ∼ t→∞ Φ v = r √ 2t , λ , Φ(v, λ) = n≥1 Φ n (v) λ n n! . ( S6
)
This scaling is based on observations originating from numerical simulations. In addition it is compatible with the known scaling of the cumulant generating function (S4). Finally, we consider the "modified centered correlations",

f µ,r (λ, t) ≡ (1 -η Xt+µ )η Xt+r e λXt ⟨e λXt ⟩ - (1 -w µ )w r-µ if µr > 0 , (1 -w µ )w r if µr < 0 , (S7)
At large time, the leading term is in t -1/2 and the sub-leading term in t -1 with the same diffusive scaling as for the profiles,

f µ,r (λ, t) = 1 √ 2t F µ v = r √ 2t , λ + 1 2t G µ v = r √ 2t , λ + O(t -3/2 ) . (S8)
This scaling ansatz relies again on numerical observations.

C. Equivalent description: large deviations and conditional profiles

Alternatively, we can also consider the distribution of the position X t of the tracer at time t, which satisfies a large deviation principle [2][3][4],

P(X t = x) ≃ t→∞ e - √ 2tϕ(ξ) , ξ = x √ 2t , ( S9
)
where ψ is the large deviation function. The moment generating function of the position of the tracer is the Laplace transform of the distribution. We can show that these functions are simply related by writing

e λXt = x e λx P(X t = x) . (S10)
Using the scaling form (S9), and taking the continuous limit, we have the following integral representation for large t:

e λXt ≃ dξ e √ 2t(λξ-ϕ(ξ)) . (S11)
The integral can be estimated with a saddle point approximation, which gives

e λXt ≃ e √ 2t(λξ * -ϕ(ξ * )) , ϕ ′ (ξ * ) = λ . (S12)
The two functions ϕ and ψ are thus related by a Legendre transform

ψ(λ) = λξ * (λ) -ϕ(ξ * (λ)) , ϕ ′ (ξ * (λ)) = λ , (S13)
or equivalently,

ϕ(ξ) = λ * (ξ)ξ -ψ(λ * (ξ)) , ψ′ (λ * (ξ)) = ξ . (S14)
In this language, it is natural to introduce the conditional profiles

⟨η Xt+r |X t = x⟩ = P(η Xt+r = 1|X t = x) (S15)
which give the probability that a site located at a distance r from the tracer is occupied, given that the tracer is located at x. We can show the equivalence between the generalised profiles (S5) and the conditional profiles (S15) by writing

w r (λ, t) = η Xt+r e λXt ⟨e λXt ⟩ = x e λx P(η Xt+r = 1|X t = x)P(X t = x) x e λx P(X t = x) . ( S16
)
Defining the large t asymptotic form of the conditional profiles as

P(η Xt+r = 1|X t = x) ≃ t→∞ Φ ξ = x √ 2t , v = r √ 2t , ( S17
)
we obtain

w r (λ, t) ≃ Φ(λ, v) = dξ Φ(ξ, v) e √ 2t(λξ-ϕ(ξ)) dξ e √ 2t(λξ-ϕ(ξ)) ≃ Φ(ξ * (λ), v) , (S18) or equivalently, Φ(ξ, v) = Φ(λ * (ξ), v) . (S19)
These two functions being equivalent, we will drop the variables λ and ξ and simply denote them both by Φ(v).

D. Equations at arbitrary time

Using Eqs. (S2), one obtains the following equations for the time-evolution of the cumulant-generating function and of the generalized profiles.

∂ t ψ = 1 2 (e λ -1)(1 -w 1 ) + (e -λ -1)(1 -w -1 ) , ( S20
)
∂ t w r = 1 2 ∆w r -B ν ∇ -ν w r + 1 2 µ=±1 e µλ f µ,r+µ -f µ,r (r ̸ = ±1) (S21) ∂ t w ±1 = 1 2 ∇ ± w ±1 + B ± w ±1 + 1 2 e ±λ f ±1,±2 -f ∓1,±1 (S22)
where ν is the sign of r, the gradients are ∇ µ u r = u r+µ -u r , ∆u r = u r+1 -2u r + u r-1 and

B ± = ∂ t ψ e ±λ -1 . ( S23
)
In addition, the generalized profiles at large distance are equal to the density: lim r→±∞ w r = ρ.

E. Hydrodynamic equations at large time

Using the scalings of Section I B into the equation (S20), we first obtain at order 0 in t,

(e λ -1)(1 -Φ(0 + )) + (e -λ -1)(1 -Φ(0 -)) = 0 , ( S24
)
which we can rewrite as

1 -Φ(0 -) 1 -Φ(0 + ) = - e λ -1 e -λ -1 = e λ . ( S25
)
Similarly, using these same scalings into (S21), we first obtain at order 1/ √ t,

F -1 (v) = e λ F 1 (v) , ( S26
)
as well as the following hydrodynamic equations for the generalized profiles at order 1/t,

Φ ′′ (v) + 2(v + b ν )Φ ′ (v) + C(v) = 0 , b ν = ν ψ e νλ -1 (S27) C(v) = (e λ -1)F ′ 1 (v) + µ=±1 (e µλ -1)G µ (v) , (S28)
with ν the sign of v. Finally, using the scalings of Section I B into Eq. (S22), we obtain the boundary condition

Φ ′ (0 ± ) ± 2 ψ e ±λ -1 Φ(0 ± ) = 0 , ( S29
)
which is completed by a second boundary condition at infinity,

Φ(v) -----→ v→±∞ ρ . (S30)
We stress that these equations are exact in the hydrodynamic limit considered here.

F. Known results on the GDP generating function

In this Section, we focus on the bulk equation (S27). We recall here the results of [1] both in the high and low density regimes, in which this equation simplifies. This will be the starting point to tackle the arbitrary density case.

High density

In the high density limit ρ → 1, we define

Φ(v) = lim ρ→1 Φ(v) 1 -ρ . (S31)
In this limit, the higher order correlations C(v) can be neglected, and the bulk equation (S27) becomes [1] Φ′′ (v) + 2v Φ′ (v) = 0. (S32)

Low density

In the opposite limit of low density, ρ → 0, one should keep z = ρr and τ = ρ 2 t constant. With these scalings, we define

Φ(v, λ) = lim ρ→0 Φ(v, λ = ρ λ) ρ , β( λ) = lim ρ→0 ±b ± (λ = ρ λ) ρ = lim ρ→0 ψ(ρ λ) ρ λ . ( S33
)
In this limit, the bulk equation (S27) is not closed. In [1], a closure relation was found, which gives

Φ′′ (v) + 2(v + ξ) Φ′ (v) = 0, ( S34
)
with ξ the (rescaled) derivative of the cumulant-generating function with respect to its parameter,

ξ ≡ β + λ dβ d λ = d d λ ( λβ) . ( S35
)

Lowest order

At order 1 in λ, the equation (S27) is also closed [1]

Φ ′′ 1 (v) + 2vΦ ′ 1 (v) = 0, ( S36
)
and its solution gives

Φ 1 (v ≷ 0) = ± 1 -ρ 2 erfc |v|. ( S37
)

II. SEP AT ARBITRARY DENSITY

A. A closed integral equation

The bulk equation (S27), valid at arbitrary density, is not closed: in addition to the functions ψ and Φ(v) of interest, it involves the unknown function C(v). We thus look for a closed equation for Φ, aiming to extend the bulk equations obtained at high density (S32) and low density (S34) to arbitrary density. More precisely, we guess that this closed equation is of the form

Φ ′′ (v) + 2(v + ξ)Φ ′ (v) = ? , ( S38
)
with a right hand side to be determined, and which vanishes in both limits ρ → 0 and ρ → 1. For this equation to be closed, ? should be expressed in terms of Φ(v) and the parameters λ and ξ only. Furthermore, we should also have ? = 0 at order 1 in λ, because of (S36). We thus expect that, at order n, the r.h.s. ? will act as a source term for the determination of Φ n , by involving only the profiles Φ m with m < n. Furthermore, the resulting equation, combined with the boundary conditions (S25,S29,S30), should also reproduce the cumulants of the tracer's position obtained recently in [2,3]. An interesting feature of these cumulants is that they involve nontrivial factors √ 2 (for the fourth cumulant κ 4 ) and √ 3 (for κ 6 ), which cannot be produced by products or derivatives of Φ 1 (S37). These factors can however be obtained by considering half-convolutions of Φ 1 with itself, such as

∞ 0 dz Φ ′ 1 (v + z)Φ ′ 1 (-z) = (1 -ρ) 2 √ 2π e -v 2 2 erfc v √ 2 . (S39)
We also have some constraints on how the different parameters (λ, ψ and ξ) should appear in the desired equation.

For instance, λ is explicitly involved in the hydrodynamic equations of Section I E only through expressions of the form e ±λ -1, so we expect that only these kind of expressions appear. In the low density equation (S34), ψ does not appear explicitly, only its derivative ξ = d ψ dλ is involved, so we expect the same to happen at arbitrary density. Finally, the equation we write should have a "proper scaling" with time. Indeed, the bulk equation (S27) (which we aim to replace) is obtained by expanding at order 1/t the microscopic equation (S21), so it should be the same for this new equation. For instance, the functions Φ, Φ ′ and Φ ′′ are respectively of orders t 0 , t -1/2 and and t -1 . The scaling argument v = r/ √ 2t is of order t -1/2 , and the same scaling holds for ξ. One can thus check that, with these scalings, the l.h.s. of Eq. (S38) indeed has the correct scaling 1/t. The same should hold for the r.h.s. ? .

Following these ideas and constraints, we obtained that the equation (valid for v > 0)

Φ ′′ (v) + 2(v + ξ)Φ ′ (v) = e -λ -1 Φ(0 -) ∞ 0 dz Φ ′′ (-z)Φ ′ (v + z) + 2ξ e -λ -1 Φ(0 -) ∞ 0 dz Φ ′ (-z)Φ ′ (v + z) + (e λ -1)(e -λ -1) Φ(0 + )Φ(0 -) ∞ 0 dz ∞ 0 dz ′ Φ ′′ (-z -z ′ )Φ ′ (z + v)Φ ′ (z ′ ) - (e -λ -1) 2 Φ(0 -) 2 ∞ 0 dz ∞ 0 dz ′ Φ ′′ (z + z ′ + v)Φ ′ (-z)Φ ′ (-z ′ ) + 2ξ (e λ -1)(e -λ -1) Φ(0 + )Φ(0 -) ∞ 0 dz ∞ 0 dz ′ Φ ′ (-z -z ′ )Φ ′ (z + v)Φ ′ (z ′ ) -2ξ (e -λ -1) 2 Φ(0 -) 2 ∞ 0 dz ∞ 0 dz ′ Φ ′ (z + z ′ + v)Φ ′ (-z)Φ ′ (-z ′ ) + • • • , (S40)
properly reproduces the known cumulant κ n for n ≤ 6. The equation for v < 0 is deduced from the symmetry Φ(-v, λ) = Φ(v, -λ). This leads us to introduce the two functions

Ω + (v) = e λ -1 Φ(0 + ) Φ ′ (v) = -2 ψ Φ ′ (v) Φ ′ (0 + ) for v > 0 , (S41a) Ω -(v) = e -λ -1 Φ(0 -) Φ ′ (v) = 2 ψ Φ ′ (v) Φ ′ (0 -) for v < 0 , (S41b)
so that (S40) takes the more compact form

Ω ′ + (v) + 2(v + ξ)Ω + (v) = ∞ 0 dz Ω ′ -(-z)Ω + (v + z) + 2ξ ∞ 0 dz Ω -(-z)Ω + (v + z) + ∞ 0 dz ∞ 0 dz ′ Ω ′ -(-z -z ′ )Ω + (z + v)Ω + (z ′ ) + 2ξ ∞ 0 dz ∞ 0 dz ′ Ω -(-z -z ′ )Ω + (z + v)Ω + (z ′ ) - ∞ 0 dz ∞ 0 dz ′ Ω ′ + (z + z ′ + v)Ω -(-z)Ω -(-z ′ ) -2ξ ∞ 0 dz ∞ 0 dz ′ Ω + (z + z ′ + v)Ω -(-z)Ω -(-z ′ ) + • • • (S42)
We can rewrite this expression in terms of the matrix operator L defined as

L = L ++ L +- L -+ L -- , (S43) with (L ++ f )(v) = ∞ 0 dz Ω -(-z)f (v + z) , (L +-f )(v) = ∞ 0 dz Ω + (v + z)f (-z) , ( S44a
) (L -+ f )(v) = - ∞ 0 dz Ω -(v -z)f (z) , (L --f )(v) = - ∞ 0 dz Ω + (z)f (v -z) . ( S44b
)
Applying this operator to the column vector (Ω + 0) T , we get

L Ω + 0 =     ∞ 0 dz Ω -(-z)Ω + (v + z) - ∞ 0 dz Ω -(v -z)Ω + (z)     , (S45)
whose first component appears in our equation (S42). Applying L 2 to the same vector, we get for the first component,

∞ 0 dz ∞ 0 dz ′ Ω + (z + z ′ + v)Ω -(-z)Ω -(-z ′ ) - ∞ 0 dz ∞ 0 dz ′ Ω -(-z -z ′ )Ω + (z + v)Ω + (z ′ ) , ( S46
)
which corresponds to the next terms in Eq. (S42). Finally, after some integration by parts, and using that

Ω + (0) = -Ω -(0) , ( S47
)
we can write (S42) as

2vΩ + (v) + (∂ v + 2ξ) (1 + L ) -1 Ω + 0 (v) 1 + Ω + (v) (1 + L ) -1 Ω + 0 (0) 1 = 0 . (S48)
Similarly, the equation for v < 0 takes the form

2vΩ -(v) + (∂ v + 2ξ) (1 + L ) -1 0 Ω - (v) 2 + Ω -(v) (1 + L ) -1 0 Ω - (0) 2 = 0 . (S49)
In order to simplify these equations, we introduce another set of functions M ±,± defined by

(1 + L ) -1 Ω + 0 0 Ω - = M ++ M +- M -+ M -- . ( S50
)
Expanding the above expression in powers of L , we notice that

Ω + (v) = M ++ (v) -M +-(v) , Ω -(v) = M --(v) -M -+ (v) . (S51)
Using these results, we can rewrite the equations (S48,S49) as

(∂ v + 2(v + ξ) + M ++ (0))M ++ = (2v + M ++ (0))M +-, (S52) (∂ v + 2(v + ξ) + M --(0))M --= (2v + M --(0))M -+ . (S53)
We can obtain a closed system of equations for M by rewriting (S50) as

L M ++ M +- M -+ M -- + M +-M +- M -+ M -+ = 0 , (S54)
where the operator L can be written in terms of M only by using (S51). The two columns are identical, therefore this only gives two equations, which are

M +-(v) = - ∞ 0 dz [M ++ (v + z)M --(-z) -M +-(v + z)M -+ (-z)] , for v > 0 , (S55) M -+ (v) = ∞ 0 dz [M --(v -z)M ++ (z) -M +-(z)M -+ (v -z)] , for v < 0 . (S56)
In order to proceed further, it is instructive to look for a perturbative solution of these equations. By definition (S41), Ω ± (and thus M ±,± ) is small when ψ is small. We find that the solutions of (S52,S53,S55,S56) at first orders in λ can be conveniently expressed as

M ++ (v) = -ω e ξ 2 Ω (1)
+ (v) + 2 -ω e ξ 2 2 Ω (2)
+ (v) + 3 -ω e ξ 2 3 Ω (3)
+ (v) + O(ω 4 ) , (S57a) M +-(v) = -ω e ξ 2 2 Ω (2) + (v) + 2 -ω e ξ 2 3 Ω (3)
+ (v) + O(ω 4 ) , ( S57b
)
where we introduced the parameter ω defined from ψ as

ω = 2 ψ + √ 2π ψ2 + 2π 9 (9 -4 √ 3) ψ3 + O( ψ4 ) ⇒ ψ = ω 2 - √ π 4 √ 2 ω 2 + π 6 √ 3 ω 3 + O(ω 4 ) (S58)
and Ω

(1)

± (v) = ±e -(v+ξ) 2 , (S59a) Ω (2)
± (v) = ± 1 2 π 2 e -1 2 (v+2ξ) 2 erfc ± v √ 2 , (S59b) Ω (3)
± (v) = ± π 4 √ 3 e -1 3 (v+3ξ) 2 erfc ± 2 3 v + erfc ± v √ 6 -4T v √ 3 , √ 3 , ( S59c
)
with T the Owen's T-function defined as [5] T(h, a) = 1 2π

a 0 e -h 2 2 (1+x 2 ) 1 + x 2 dx . (S60)
This leads us to write the general form as

M ++ (v) = n≥1 n -ω e ξ 2 n Ω (n) + (v) , M +-(v) = n≥2 (n -1) -ω e ξ 2 n Ω (n) + (v) , (S61) so that, from (S51), Ω + (v) = n≥1 -ω e ξ 2 n Ω (n) + (v) , (S62)
and similarly for Ω -. Plugging these expressions into (S55), we obtain

Ω (n) + (v) = - n-1 p=1 ∞ 0 dz Ω (p) + (v + z)Ω (n-p) - (-z) . ( S63
)
Multiplying by (-ωe ξ 2 ) n and summing over n, we get

Ω + (v) --ω e ξ 2 Ω
(1)

+ (v) = - ∞ 0 dz Ω + (v + z)Ω -(-z) . (S64)
From the expression of Ω

+ (S59a), this becomes

Ω + (v) + ω e -(v+ξ) 2 +ξ 2 + ∞ 0 dz Ω + (v + z)Ω -(-z) = 0 , (S65)
and similarly we obtain the equation for Ω -from (S56):

Ω -(v) -ω e -(v+ξ) 2 +ξ 2 - ∞ 0 dz Ω -(v -z)Ω + (z) = 0 . (S66)
Equations similar to (S65,S66) can be found in Ref. [6]. In this paper, the authors show that these two coupled equations are equivalent to the two independent linear equations

Ω + (v) = -ω e -(v+ξ) 2 +ξ 2 -ω 0 -∞ dz Ω + (z) e -(v-z+ξ) 2 +ξ 2 , (S67) Ω -(v) = ω e -(v+ξ) 2 +ξ 2 -ω ∞ 0 dz Ω -(z) e -(v-z+ξ) 2 +ξ 2 , ( S68
)
which correspond to the Eq. ( 6) given in the main text. These equations give Ω ± for all v ∈ R (including the analytic continuations of Ω + to v < 0 and Ω -to v > 0 which both appear explicitly in the integrals).

B. Solution

Equations similar to (S67,S68) are solved in [7], but restricted to v > 0, and with Ω ± (∓v) = 0 for v < 0. We can nevertheless use these results to express Ω ± in this domain. The results are given in terms of the one-sided Fourier transforms:

Ω(+) ± (k) = ∞ 0 dv Ω ± (v) e ikv , Ω(-) ± (k) = - 0 -∞ dv Ω ± (v) e ikv . ( S69
)
From Ref. [7], we obtain the Fourier transforms of the analytic continuations of Ω ± :

Ω(-) + (k) = 1 -exp   1 2 n≥1 (-ω √ π e -1 4 (k+2iξ) 2 ) n n erfc - √ n ξ - ik 2   , (S70a) Ω(+) -(k) = 1 -exp   1 2 n≥1 (-ω √ π e -1 4 (k+2iξ) 2 ) n n erfc √ n ξ - ik 2   . ( S70b
)
We can obtain the Fourier transforms on the original functions by taking the Fourier transform of (S67), which gives

Ω(+) + (k) = K(k) + Ω(-) + (k)(1 -K(k)) , with K(k) = -ω √ π e -1 4 (k+2iξ) 2 . (S71)
We finally obtain

∞ 0 Ω + (v)e ikv dv = 1 -exp -Z + ω, ξ - ik 2 , (S72) 0 -∞ Ω -(v)e ikv dv = exp -Z -ω, ξ - ik 2 -1 , ( S73
)
where we have introduced

Z ± (ω, ξ) = 1 2 n≥1 (-ω √ π e ξ 2 ) n n erfc ± √ nξ . (S74)
We can obtain the values of Ω ± (0) = ∓2 ψ from these expressions. We set k = is in (S67) and let s → ∞, this gives

Ω + (0) = -2 ψ = 1 √ π Li 3 2 (-ω √ π) , ( S75
)
which expresses ψ in terms of ω. Note that this expression is consistent with the first orders obtained previously (S58).

C. Expansion in λ 1. Computation of the cumulants

Plugging the expansion of ψ (S4) into (S75), we can deduce the expansion of ω in powers of λ, in terms of the cumulants κ n . Combining the solution (S72) (for k = 0) with the definition of Ω + (S41), and the boundary condition (S29), we obtain

1 -exp   - 1 2 n≥1 (-ω √ π e ξ 2 ) n n erfc √ nξ   = (e λ -1) ρ -Φ(0 + ) Φ(0 + ) . ( S76
)
Since ξ = d ψ dλ , this equation gives Φ n (0 + ) in terms of the cumulants κ j for j ≤ n. We can proceed similarly for Φ n (0 -) using (S73). Finally, using the last relation (S25), we can determine the cumulants κ n and thus Φ n (0 ± ). Due to the symmetry v → -v with λ → -λ, all the odd order cumulants vanish. For the even order ones, we get for instance,

κ 2 = 1 -ρ ρ √ π , ( S77
)

κ 4 = (1 -ρ) π 3/2 ρ 3 12(1 -ρ) 2 -π(3 -3(4 - √ 2)ρ + (8 -3 √ 2)ρ 2 ) , ( S78
)

κ 6 = (1 -ρ) π 5/2 ρ 5 30π 2 9 √ 2 -20 ρ 2 + 60 -18 √ 2 ρ -15 (1 -ρ) 2 -π 2 8 -17 + 15 √ 2 -5 √ 3 ρ 4 + 480 -300 √ 2 + 80 √ 3 ρ 3 +5 -114 + 45 √ 2 -8 √ 3 ρ 2 -45 √ 2 -6 ρ -45 + 1020(1 -ρ) 4 , (S79)
which coincide with the cumulants obtained from the CGF computed in [2]. This is expected because the equation (S40) has been constructed in order to reproduce these cumulants. We have further checked with Mathematica, up to n = 10, that the next cumulants κ n obtained by our procedure also coincide with those obtained from [2]. This provides a nontrivial validation of our integral equations (S67,S68).

Computation of the generalized profiles

Having determined Φ n (0 ± ) and κ j for j ≤ n, and thus Φ ′ n (0 ± ) from the boundary condition (S29), we can express Φ ′ (v) in terms of Ω ± (v). Expanding the explicit solutions (S72,S73) in powers of ω and computing the inverse Fourier transform, we recover the expansion (S62) with Ω (n) ± (v) given by (S59) for n ≤ 3, but we can also access higher orders via an inverse Fourier transform. Having expressed ω in terms of the cumulants determined previously via (S75), we thus have the expansion of Φ ′ (v) in powers of λ. After integration, we obtain in particular

Φ 0 (v) = ρ , (S80a) Φ 1 (v) = 1 -ρ 2 erfc(v) , (S80b) Φ 2 (v) = (1 -ρ)(1 -2ρ) 2ρ erfc(v) - 2 π (1 -ρ) 2 ρ e -v 2 , (S80c) Φ 3 (v) = (1 -ρ) 2(3 + π)ρ 2 -(12 + π)ρ + 6 + πρ(1 -ρ) 2πρ 2 erfc(v) + 3(1 -ρ) 2 2(1 -ρ)v - √ π(1 -2ρ) π 3/2 ρ 2 e -v 2 - 3(1 -ρ) 2 4ρ erfc v √ 2 2 , (S81) Φ 4 (v) = (1 -ρ)(1 -2ρ) 24(1 -ρ) 2 -π(3(1 -2ρ) 2 + 4(1 -ρ)(1 -2ρ) + 18ρ(1 -ρ) -4) 2πρ 3 erfc(v) + 4(1 -ρ) 2 3 √ π(1 -ρ)(1 -2ρ)v -4(1 -ρ) 2 (4 + v 2 ) + π(5(1 -ρ)(1 -2ρ) + 3(1 + √ 2)ρ(1 -ρ) -2) π 2 ρ 3 e -v 2 + 12 √ 2 (1 -ρ) 3 πρ 2 e -v 2 2 erfc v √ 2 -3 (1 -ρ) 2 (1 -2ρ) 2ρ 2 erfc v √ 2 2 . (S82) Φ 5 (v) = - (1 -ρ) π 2 ρ 4 + 440(1 -ρ) 4 + 30π 2ρ 2 √ 2 -3 ρ -2 √ 2 + 6 -3 (1 -ρ) 2 240π 2 ρ 4 erfc(v) + (1 -ρ) 2 24π 5/2 ρ 4 4 √ π(1 -ρ) 2 19ρ + (4ρ -2)v 2 -11 + 12π ρ 2 √ 2 -5 ρ -2 √ 2 + 8 -2 (1 -ρ)v -3π 3/2 (1 -2ρ) 2ρ √ 2 -3 ρ - √ 2 + 5 -3 + 4(1 -ρ) 3 v 2v 2 + 27 e -v 2 + (1 -ρ) 3 ( √ π(1 -2ρ) -2(1 -ρ)v) 2 √ 2π 3/2 ρ 3 e -v 2 2 erfc v √ 2 + (1 -ρ) 2 π 2ρ 2 -4ρ + 1 -12(1 -ρ) 2 32πρ 3 erfc v √ 2 2 - (1 -ρ) 3 4 √ 3π ρ 2 +∞ v dz e -z 2 3 erfc 2 3 z + erfc z √ 6 -4T z √ 3 , √ 3 , (S83)
with the Owen-T function defined in (S60).

Conservation relation

Using the results above, the conservation relation

∞ 0 (Φ(v) -ρ)dv - 0 -∞ (Φ(v) -ρ)dv = ρξ (S84)
holds up to O(λ 6 ), and non perturbatively in λ (numerically).

D. Extensions

Step density profile

Our formalism can be extended to the case of an initial step density ρ -for v < 0 and ρ + for v > 0, by only changing the boundary condition at infinity (S30) into Φ(±∞) = ρ ± . Unlike the constant density case, the odd order cumulants do not vanish here. We can still apply the procedure described in Section II C, which gives that κ 1 is solution of

ρ - 1 + √ π κ 1 e κ 2 1 erfc(-κ 1 ) = ρ + 1 - √ π κ 1 e κ 2 1 erfc(κ 1 ) , ( S85
)
and the higher order cumulants are expressed in terms of κ 1 . For instance,

κ 2 = κ 2 1 2πe 2κ 2 1 κ 1 erfc √ 2κ 1 + 4πe 2κ 2 1 κ 1 ρ + ρ 2 + -3ρ -ρ + + 2ρ - (ρ --ρ + ) 3 - √ 2π , (S86)
These expressions coincide with those given in [2]. We additionally obtain the profiles Φ n (v), for instance

Φ 0 (v) = ρ + 2 erfc(-v -κ 1 ) + ρ - 2 erfc(v + κ 1 ) , (S87) Φ 1 (v) = 2 √ πe κ 2 1 κ 1 ρ -(1 -ρ + ) -(ρ -ρ + ) 2 2 (ρ --ρ + ) erfc (κ 1 + v) - e -(κ1+v) 2 4 √ πe 2κ 2 1 κ 3 1 (ρ --ρ + ) 3 erfc √ 2κ 1 + 2ρ + ρ 2 + + ρ -(2 -3ρ + ) -2 √ 2κ 2 1 (ρ --ρ + ) 3 2 (ρ --ρ + ) 2 + 1 2 √ πe κ 2 1 κ 1 (ρ + -ρ -) 4T √ 2κ 1 , κ 1 + v κ 1 -4T 2κ 1 + v, v 2κ 1 + v + erfc 2κ 1 + v √ 2 -erfc (κ 1 ) . (S88)
for v > 0 and κ 1 > 0.

Another observable: the current through the origin

We now consider another observable, which is the flux of particles through the origin during a time t [8,9], which we can write as

Q t = r≥1 (η r (t) -η r (0)) . (S89)
We consider the cumulant generating function

ψ Q (λ) = √ 2t ψQ = ln e λQt , ( S90
)
and the generalised profiles for the current

w Q;r (t) = η r e λQt ⟨e λQt ⟩ . (S91)
And for large t, we have

w Q;r (t) ≃ t→∞ Φ Q (v) = n≥0 Φ Q;n (v) λ n n! , v = r -1 2 √ 2t . ( S92
)
We still define the functions Ω ± (v) as

Ω ± (v) = ∓2 ψQ Φ ′ Q (v) Φ ′ Q (0 ± ) . ( S93
)
It can be checked that these functions again verify the integral equations (S67,S68), but with ξ = 0:

Ω + (v) = -ω Q e -v 2 -ω Q 0 -∞ dz Ω + (z) e -(v-z) 2 , (S94) Ω -(v) = ω Q e -v 2 -ω Q ∞ 0 dz Ω -(z) e -(v-z) 2 , (S95) which still imply that ψQ = - 1 2 √ π Li 3 2 (-ω Q √ π) . (S96)
In order to obtain the boundary conditions satisfied by Φ, we write the time evolution of the CGF from the master equation,

∂ t ln e λQt = 1 2 (e λ -1) η 0 (1 -η 1 )e λQt ⟨e λQt ⟩ + (e -λ -1) η 1 (1 -η 0 )e λQt ⟨e λQt ⟩ . ( S97
)
We proceed similarly for the generalized profiles

∂ t w Q;0 = (e -λ -(e -λ -1)w Q;0 ) ∂ t ln e λQt e -λ -1 + w Q;-1 -w Q,0 2 , ( S98
)
∂ t w Q;1 = (e λ -(e λ -1)w Q;1 ) ∂ t ln e λQt e λ -1 + w Q;2 -w Q;1 2 . ( S99
)
Taking the hydrodynamic limit, we get at leading order

(e λ -1)Φ Q (0 -)(1 -Φ Q (0 + )) + (e -λ -1)Φ Q (0 + )(1 -Φ Q (0 -)) = 0 , (S100) Φ ′ Q (0 -) = 2 ψQ 1 1 -e λ -Φ Q (0 -) , Φ ′ Q (0 + ) = -2 ψQ 1 1 -e -λ -Φ Q (0 + ) . (S101)
These boundary conditions, combined with the solution of the equations (S94,S95) yield

ω Q √ π = ρ -(1 -ρ + )(e λ -1) + ρ + (1 -ρ -)(e -λ -1) . (S102)
With (S96), this allows us to recover the result of Derrida and Gerschenfeld [8,9] on the cumulant generating function ψ. Additionally, we obtain the profiles Φ Q;n (v). For instance

Φ Q;1 (v) = ρ(1 -ρ) 2 erfc(v) , (S103a)
Φ Q;2 (v) = ρ(1 -ρ)(1 -2ρ) 2 erfc(v) , (S103b)
Φ Q;3 (v) = ρ(1 -ρ)(1 -3ρ + 3ρ 2 ) 2 erfc(v) -3 ρ 2 (1 -ρ) 2 4 erfc v √ 2 2 . ( S103c
)

Another observable: a generalized current

For a given position x, we consider the generalized current

J t (x) = r≥1 (η r+x (t) -η r (0)) (S104)
studied in [2,3]. It measures the difference between the number of particles on the positive axis at t = 0 and the number of particles at the right of x at time t. In this Section, we show that the cumulant generating function of the generalized flux which has been recently computed in [2,3] can be retrieved from our approach. One subtlety is that at t = 0 the observable J 0 (x) ̸ = 0, which results in additional difficulties due to the contribution of this (random) initial value. To circumvent this difficulty, we consider the observable J t (x t ), with x t = ⌊ξ √ 2t⌋ the integer part of ξ √ 2t, which now verifies J 0 (x 0 ) = 0. The time evolution of the cumulant generating function of this observable combines two contributions: the jumps of x t at times t n = (n/ξ) 2 /2 and the evolution at x t fixed between two jumps. Combining these two contributions, it can be shown that, for large t, Similarly, we obtain that the profiles satisfy

∂
∂ t w J;0 = (e -λ -(e -λ -1)w J;0 ) ∂ t ln e λJt(xt) -ξ √ 2t ln[1 + (e -λ -1)w J;1 (t)] e -λ -1 + w J;-1 -w J;0 2 + ξ √ 2t
e -λ w J;1 1 + (e -λ -1)w J;1 -w J;0 , (S107) and an analogous expression holds for w J;1 . For large times, we have the scalings

ln e λJt(xt) ≃ √ 2t ψJ (λ, ξ) , w r (t) ≃ t→∞ Φ J v = r -1 2 √ 2t , ξ, λ , (S108)
which we will denote by Φ J (v) for simplicity. Using these scalings, we get from (S105,S107):

(e λ -1)Φ J (0 -)(1 -Φ J (0 + )) + (e -λ -1)Φ J (0 + )(1 -Φ J (0 -)) = 0 , (S109) Φ ′ J (0 ± ) = ∓2Ψ 1 1 -e ∓λ -Φ J (0 ± ) , Ψ = ψJ -ξ ln[1 + (e -λ -1)Φ J (0 + )] . (S110) If we define Ω ± (v) = ∓2Ψ Φ ′ J (v) Φ ′ J (0 ± ) , (S111)
the solution of the integral equations (S67,S68) combined with the boundary conditions gives

ω J √ πe ξ 2 = ρ -(1 -ρ + )(e λ -1) + ρ + (1 -ρ -)(e -λ -1) , ( S112
)
and we obtain for the cumulant generating function of J t (x t ):

ψJ (λ, ξ) = - 1 2 √ π Li 3 2 (-ω J √ π) + ξ ln[1 + (e -λ -1)Φ J (0 + )] , ( S113
)
which we can combine with the expression of Φ J (0 + ) obtained by integration of Ω + to obtain

ψJ (λ, ξ) = ξ ln[1 + (e -λ -1)ρ + ] - n≥1 (-ω J √ πe ξ 2 ) n 2n e -nξ 2 √ nπ -ξ erfc( √ nξ) , ( S114
)
which coincides exactly with the expression obtained in [2,3]. This supports the exactness of our main equations (S67,S68), along with the definition (S111) of Ω ± in this case.

E. Comparison with Imamura et al for the position of the tracer

Obtaining the cumulant generating function of the tracer's position from the generalized current

We have shown in Section II D 3 that our main equations (S67,S68) allow to recover the cumulant generating function ψJ (S114) of the generalized current J t (S104) recently obtained by Imamura et al in [2,3]. In these papers, the authors deduce from this result the cumulant generating function ψ of the position of the tracer as follows.

First, taking a Legendre transform, one gets the distribution of the current J t (x), which takes the form

P(J t (x) = J) ∼ t→∞ e - √ 2t φ J (ξ,j) , with j = J √ 2t , ξ = x √ 2t , ( S115
)
where

φ J (ξ, j) = χ * (ξ, j)j -ψJ (χ * (ξ, j), ξ) , ∂ ψJ ∂χ χ=χ * = j , ψJ (χ, ξ) = lim t→∞ 1 √ 2t ln e χJt(x) . (S116)
We have replaced the parameter λ in (S114) by χ to avoid confusion with the argument λ of ψ(λ). Second, using that the number of particles to the right of the tracer is conserved, the position X t of the tracer verifies J t (X t ) = 0. The distribution of the position of the tracer can therefore be obtained as

P(X t = x) = P(J t (x) = 0) ∼ t→∞ e - √ 2t φ J (ξ,0) . (S117)
Finally, taking a Legendre transform of φ J (ξ, 0) yields the cumulant generating function ψ of X t . Having recovered the expression of ψJ (S114) derived in [2,3], we have therefore the same cumulant generating function ψ of X t as [2,3] from the procedure above.

An alternative parametrization for the cumulant generating function

Additionally, we have obtained an alternative parametrization of the cumulant generating function of X t , which takes the simple form (S75) in terms of the parameter ω. This parameter is deduced from ρ and λ as follows.

Using the relation (S29) between Φ ′ (0 ± ) and Φ(0 ± ), combined with the solutions (S72,S73) we obtain the relation (S76) which allows to determine Φ(0 + ), and a similar one for Φ(0 -). Combining these results with (S25), we obtain the first equation 2(cosh λ -1) = ρ(e λ -1) 2 e λe -Z+(ω,ξ) + ρ(e -λ -1) 2 e -λe -Z-(ω,ξ) , (S118) which relates ξ and ω to ρ and λ.

The second equation needed to fully determine ξ and ω is the relation ξ = d ψ dλ . Using this expression in (S118), together with (S75), we obtain a nonlinear differential equation, whose solution yields ψ(λ).

A more convenient parametrization can be obtained by using the conservation relation (S84) instead of ξ = d ψ dλ . Combined with the solutions (S72,S73), this equation can be rewritten as

ξ = 1 2 ∂ ξ ln e λ -e -Z+(ω,ξ) e -λ -e -Z-(ω,ξ) , ( S119
)
which can be written as an algebraic equation using the expression of Z ± (S74).

To summarize, given λ and ρ, the parameters ξ and ω are obtained by solving (S118) and (S119). The cumulant generating function ψ(λ) is then straightforwardly deduced from (S75). Furthermore, since the large deviations function ϕ(ξ) of the tracer's position is deduced from ψ by the Legendre transform (S14), we straightforwardly obtain from this parametrization that ϕ(ξ) = λξ -ψ(λ). This is an alternative route to the one followed in [2,3]. We have checked numerically that the two parametrizations give the same result, validating our approach.

F. Comparison with Derrida et al for the current through the origin

Combining our intermediate results (S96,S102) on the current through the origin Q t , we obtain the cumulant generating function

ψQ (λ) = lim t→∞ 1 √ 2t ln e λQt = - 1 2 √ π Li 3 2 -ρ -(1 -ρ + )(e λ -1) -ρ + (1 -ρ -)(e -λ -1) , (S120)
which is the result of [8,9]. This gives another validation of our approach.

III. EXTENSION TO OTHER SINGLE-FILE SYSTEMS

A. Description of single-file systems in terms of two transport coefficients

In the language of fluctuating hydrodynamics [10], a single-file system can be described at large distance and large time by a fluctuating density field ρ(x, t) that is shown to obey the following equation,

∂ t ρ(x, t) = ∂ x D(ρ(x, t))∂ x ρ(x, t) + σ(ρ(x, t))η(x, t) , ( S121
)
where η(x, t) is a normalized Gaussian white noise uncorrelated in space and time. The quantities D(ρ) and σ(ρ) were first defined from the microscopic details of a lattice gas [10]. It is nevertheless more intuitive to define them for a system of size L between two reservoirs at densities ρ a and ρ b [11]. The number of particles transferred from left to right at time t is denoted by Q t and is shown to satisfy

lim t→∞ ⟨Q t ⟩ t = D(ρ) L (ρ a -ρ b ) if (ρ a -ρ b ) is small, lim t→∞ Q 2 t t = σ(ρ) L if ρ a = ρ b = ρ. (S122)
This connection further extends to the generalized density profiles associated to these two observables. In order to show it explicitly, it is more convenient to use the conditional profiles discussed in Section I C. We consider for the RAP the mean density conditioned on the tracer's position

⟨ρ(y + x 0 (t), t)|x 0 (t) = x⟩ ≃ t→∞ Φ(RAP) v = y √ t , ξ = x √ t . (S129)
Similarly, for the KMP model, we introduce the mean density conditioned on the value of the current

⟨ν(z, t)|Q t [ν] = Q⟩ ≃ t→∞ Φ(KMP) Q u = z √ t , q = Q √ t . ( S130
)
Since ν(z, t) is the spacing between the particles labelled by z and ρ(x, t) the density of particles at position x, the two conditional profiles are related by

Φ(RAP) (v(u), ξ) = 1 Φ(KMP) Q (u, -ξ) , ( S131
)
where v(u) is the position of the particle labelled by u, in the reference frame of the tracer particle. It can be obtained by writing that

x k (t) -x 0 (t) = k 0 ν(z, t)dz (S132)
from the definition of ν(z, t). Taking the average, conditioned on x 0

(t) = -Q t [ν] = ξ √ t, we obtain v(u) = u 0 Φ(KMP) Q (u ′ , -ξ)du ′ . (S133)
Equations (S131,S133) give a parametric expression for the conditional profiles of the RAP. An analogous parametrization is given in [21] for the average density in the presence of a biased tracer, but without conditioning (in our case, this would correspond to a flat density profile). Similarly to the demonstration of Section I C, we can show that these profiles are equivalent to the joint cumulants generating functions (S124,S126),

w x (λ, t) = ρ(x + x 0 (t), t)e λx0(t) e λx0(t) ≃ t→∞ Φ (RAP) (v, λ) = Φ(RAP) (v, ξ * ) , ξ * = d ψ(RAP) dλ , (S134) w Q;z (λ, t) = ν(z, t)e λQt[ν] e λQt[ν] ≃ t→∞ Φ (KMP) Q (v, λ) = Φ(KMP) Q (v, ξ * ) , ξ * = d ψ(KMP) Q dλ . (S135)
Therefore, we finally have the parametrization

Φ (RAP) (v(u), λ) = 1 Φ (KMP) Q (u, -λ) , v(u) = u 0 Φ (KMP) Q (u ′ , -λ)du ′ . (S136)
Expanding this expression in powers of λ, we can obtain the profiles Φ (RAP) n (v) of the RAP from those associated with the flux in the KMP model. This is done in Section III E below.

C. Modified equations for the GDP-generating function

All the models discussed in this paper (including the RAP via the mapping to the KMP model) are described by the situation D(ρ) = D 0 with σ ′′ (ρ) constant and σ(0) = 0, so we restrict ourselves to this case.

For the position of the tracer

At large times, the cumulant generating function of the tracer's position scales as

ψ(λ, t) ∼ t→∞ ψ(λ) 4D 0 t , (S137)
and the GDP generating function (S124) as

w r (λ, t) ∼ t→∞ Φ v = r √ 4D 0 t , λ . (S138)
We still define the functions Ω ± (v) as

Ω ± (v) = ∓2 ψ Φ ′ (v) Φ ′ (0 ± ) . ( S139
)
Our main equation becomes

Ω ± (v) = ∓ω e -(v+ξ) 2 +ξ 2 + σ ′′ (0) 4D 0 ω R ∓ Ω ± (z) e -(v-z+ξ) 2 +ξ 2 dz , ( S140
)
with ξ = d ψ dλ , and the boundary conditions

Φ ′ (0 ± ) ∓ ψ σ ′′ (0) 2D 0 Φ(0 ± ) e ∓ σ ′′ (0)λ 4D 0 -1 = 0 , 2σ ′ (0) + σ ′′ (0)Φ(0 + ) 2σ ′ (0) + σ ′′ (0)Φ(0 -) = e σ ′′ (0)λ 4D 0 . (S141)
The solution of the integral equation (S140) can be easily deduced from (S72,S73). In particular, we get from the expression of Ω + (0) that

ψ = 2D 0 σ ′′ (0) √ π Li 3 2 σ ′′ (0) √ π 4D 0 ω . (S142)

For the current through the origin

Our results on the tracer's position can be extended to the current through the origin (S125), generalizing the discussion of Section II D 2 to more general single file systems. The cumulant generating function scales as ψ Q (λ, t) = ln e λQt ∼ t→∞ ψQ (λ) 4D 0 t , (S143) and the GDP generating function (S126) as

w Q;r (λ, t) ∼ t→∞ Φ Q v = r √ 4D 0 t , λ . (S144)
Defining again Ω ± as in (S93), these functions satisfy (S140) with ξ = 0 and ω replaced by ω Q , which is related to ψQ by ψQ = 2D 0 σ ′′ (0)

√ π Li 3 2 σ ′′ (0) √ π 4D 0 ω Q . ( S145
)
The boundary conditions now become (S148c)

Φ ′ Q (0 ± ) = ∓2 ψQ σ ′ (0) 2D 0 1 1 -e ∓ σ ′ (0) 2D 0 λ + σ ′′ (0) 4D 0 Φ Q (0 ± ) , Φ Q (0 + )(2σ ′ (0) + σ ′′ (0)Φ Q (0 -)) Φ Q (0 -)(2σ ′ (0) + σ ′′ (0)Φ Q (0 + )) = e σ ′ (

For the current through the origin

In the case of the current, we obtain the general expression of the cumulant generating function ψ(KMP)

Q (λ) = D 0 √ πσ 0 Li 3 2 σ 0 ρλ 2D 0 2 , ( S149
)
which coincides with the one given in [9], and also the profiles (S151)

Φ (KMP) Q;1 (v) = ρ 2 σ 0 4D 0 erfc(v) , (S150a)
Similarly, the profiles for the RAP can be deduced from the one associated with the current in the KMP model (S150) by setting D 0 = µ 1 /2 and σ 0 = µ 1 µ 2 /(µ 1 -µ 2 ), and using the parametrization (S131). This gives We now compare our results with those obtained in the formalism of the Macroscopic Fluctuation Theory [22][23][24][25][26]. We are interested in the study of the single-file system at a large time T . We introduce a new density ρ, constructed from ρ(x, t) introduced in Section III A by rescaling the time and position:

Φ (RAP) 1 (v) = µ 2 2(µ 1 -µ 2 ) erfc(v) , (S152a) Φ (RAP) 2 (v) = µ 2 2 2πρ(µ 1 -µ 2 ) 2 π erfc(v) 2 -2π 1 + v e -v 2 √ π erfc(v) -2e -v 2 + 2e -2v 2 , (S152b) Φ (RAP) 3 (v) = 3 4π 2 µ 3 1 ρ 2 (µ 1 -µ 2 ) 3 π 2 erfc(v)
ρ(x, t) = ρ u = x √ T , τ = t T . ( S153
)
The probability to evolve from an density ρ0 at t = 0 to a density ρ(u, 1) at t = T is given by [9]: For large T , the integral in (S158) is dominated by the minimum of S + F -λY , taken as a function of (ρ, H). We denote this minimum (q, p). These functions satisfy the evolution equations [16] ∂ τ q = ∂ u [D(q)∂ u q] -∂ u [σ(q)∂ u p] , (S160)

∂ τ p = -D(q)∂ 2 u p - 1 2 σ ′ (q)(∂ u p) 2 , ( S161

FIG. 2 .

 2 FIG.2. SEP. Generalized density profiles (GDP) of order: (a) n = 1 and n = 3 at density ρ = 0.5. Solid lines: numerical simulations at t = 3000. Dashed lines: theoretical predictions (equations (S79,S80) in SI); (b) n = 2 at densities ρ = 0.25, 0.5 and 0.75. Solid lines: numerical simulations at t = 3000. Dashed lines: theoretical predictions (equation (S79) in SI). Note that, interestingly Φ2 exhibits a minimum at a distance v > 0 for ρ < 1/2, which disappears for ρ > 1/2. (c) GDP-generating function (or equivalently, conditional profiles) obtained from the theoretical prediction (6) (dashed) compared to the numerical resolution of the MFT equations (solid line) described in Section IV.B in SI. (d) Non-monotony of the conditional profiles Φ(0 + ) in front of the tracer, as a function of the rescaled tracer's position ξ, for different densities ρ.

FIG. 3 .

 3 FIG. 3. Extensions: (a) The out of equilibrium situation of an initial step density. (b) Another observable, the current through the origin. (c),(d) Other representative single-file systems: the KMP model (c) and the RAP (d). The dashed lines correspond to the predictions obtained from the central equation (6). (a) GDPs at order n = 1 and n = 2 for the SEP with a step initial density ρ-= 0.7 and ρ+ = 0.2 at t = 1500 with 2000 sites. (b) GDPs for the current ⟨ηrQ n t ⟩c in the SEP for a density ρ = 0.5 at t = 900, for orders n = 1 and n = 3 (the profile for n = 2 is zero). (c) KMP model (D(ρ) = 1/2 and σ(ρ) = ρ 2 ) for ρ = 1. (c1)-(c3) First three orders of the GDPs at t = 900 and 500 sites. (c4) GDP-generating function (d) RAP for ρ = 1 (for a uniform jump distribution, D(ρ) = 1/(4ρ) and σ(ρ) = 1/ρ). (d1)-(d3) First three orders of the GDPs at t = 4000 with 5000 particles (solid lines). (d4) GDP-generating function compared to the numerical resolution of the MFT equations (solid line).

  and cumulants for the KMP modelApplying the procedure described above in Section II C to the case of the KMP model, with D(ρ) = D 0 and σ(ρ) = σ 0 ρ 2 , we can obtain the cumulants and profiles for this model.

1 .

 1 For the position of the tracerFor the tracer's position, we obtain for instance at first orders κ

Φ 1 4D

 1 and cumulants for the RAPWe straightforwardly obtain the cumulant generating function ψ(RAP) of the position of a tracer in the RAP from the one of the current in the KMP model (S149) via the relation (S128), which yieldsψ(RAP) (λ) = lim t→∞ (ρ)t ψ (RAP) (λ, t) = ρ(µ 1 -µ 2 )

P(ρ 0 1 0 1 )

 011 (u) -→ ρ(u, 1)) = D[ρ(u, τ )]D[H(u, τ )] e - √ T S[ ρ,H] , (S154)where the action S readsS[ρ, H] = du dτ H∂ τ ρ + D(ρ)∂ u ρ∂ u H -σ(ρ) 2 (∂ u H) 2 . (S155)The distribution of the initial condition ρ0 isP[ρ 0 ] ≃ e - √ T F [ ρ0] ,, the moment generating function of the tracer's position is given by[16] eλX T ≃ D ρ0 D[ρ(u, τ )]D[H(u, τ )] e - √ T (S[ ρ,H]+F [ ρ0]-λY [ ρ]) ,(S158)where Y [ρ] = X T / √ T is the rescaled position of the tracer, which is deduced from (S123): -ρ(u, 0)) du . (S159)

  )with the terminal condition for pp(u, τ = 1) = BΘ(u -Y ) , B = λ q(Y, 1) , (S162)and the initial condition for q, expressed in terms of p(u, 0):p(u, 0) = BΘ(u) + [1], our generalised density profiles (S124) can be deduced from the MFT solution q(u, τ ), sincew r (λ, T ) ≃ D ρ0 D[ρ(u, τ )]D[H(u, τ )] ρ(Y [ρ] + r/ √ T , 1) e - √ T (S[ ρ,H]+F [ ρ0]-λYT [ ρ]) D ρ0 D[ρ(u, τ )]D[H(u, τ )] e - √ T (S[ ρ,H]+F [ ρ0]-λYT [ ρ]) , (S164)

observables, nonequilibrium situations), following these lines, are described in SI.

We list below the expressions of D(ρ) and σ(ρ) for a few models considered here.

Model

D(ρ) σ(ρ)

Symmetric exclusion process [12]

Hard Brownian particles [12] D 0 2D 0 ρ Kipnis-Marchioro-Presutti [13] D 0 σ 0 ρ 2

Random average process [14,15]

D 0 is the diffusion coefficient of an individual particle, σ 0 = 2aD 0 with a the lattice constant of the KMP model [13], and µ k are the moments of the probability law of the jumps in the RAP [15].

We are interested in the position X t of the tracer, which we can define as [16] Xt

by expressing that the number of particles to the right of the tracer is conserved. We define the associated cumulant generating function and the generalized profiles ψ(λ) = ln e λXt , w r (λ, t) = ρ(X t + r, t)e λXt ⟨e λXt ⟩ . (S124)

We will also consider the current through the origin, which is expressed as

and the associated profiles ψ Q (λ) = ln e λQt , w Q;r (λ, t) = ρ(r, t)e λQt ⟨e λQt ⟩ . (S126)

B. Mapping the RAP to the KMP model

The random average process [14,17,18] consists of particles on an infinite line, placed at positions x k (t) with initial density ρ. The particles are allowed to move to a random fraction of the distance to the next one, either to the left or to the right with rate 1 2 . Only the first two moments µ 1 and µ 2 of the distribution of this random fraction are relevant in the hydrodynamic limit, in which the system is described by the coefficients D(ρ) and σ(ρ) given in the table above.

This model can be mapped onto the Kipnis Marchioro Presutti model [19,20], which describes a one dimensional lattice where each site contains an energy [21]. At random times, the total energy of two neighbouring sites is randomly redistributed on these sites. In the hydrodynamic limit, we can replace the discrete index k by a continuous variable z and consider the density of the spacings (or energies) ν(z, t), which averages to 1/ρ. This system is then described by [15]. The displacement of the tracer particle (initially at the origin) is then given by

where Q t [ν] is the current through the origin in the KMP model. The cumulant generating function of the tracer's position in the RAP is thus directly related to the one of the current in the KMP model:

which yields from a saddle point estimate:

where we have used here the scaling with time introduced above in the case D(ρ) = D 0 for consistency, but this relation between Φ and q also holds for a general D(ρ) (without rescaling the positions with D 0 ). Aside from a few specific cases (such as the hard Brownian particles), the MFT equations cannot be solved analytically for arbitrary λ.

We thus rely both on a perturbative solution and a numerical resolution of these equations to compare with our results.

A. Perturbative expansion in λ for the SEP

We come back to the case of the SEP, corresponding to

The MFT equations can be solved perturbatively by expanding them in powers of the parameter B, defined in (S162), which appears explicitly in the equations, as

(S166)

This procedure was carried out for the first orders in [16]. The difficulty is then to relate B and λ, because the solution q(u, τ = 1) is discontinuous at u = Y , which makes it impossible to use the definition (S162). The relation between these parameters can still be found by treating B and λ independently, and minimizing the resulting cumulant generating function with respect to B. This procedure was used in [16]. Here, we use a shortcut: since Y = √ 2ξ with ξ = d ψ dλ and ψ has been determined in [2], we can use this relation to obtain B as a function of λ after inversion of (S166). The main difficulty is now to compute q n (u, τ ) at a given order n.

The solution for the first two orders has been computed in [16], and reads

where

is the heat kernel. Starting from order 3, the resolution becomes more complex. In [16], the solutions at order 3 were written in the form

where the two additional functions p3 and q3 satisfy the following inhomogeneous heat equations

with the boundary conditions

However, these equations cannot be solved analytically, and were studied numerically in [16] in order to compute the fourth cumulant κ 4 . One can indeed obtain exact integral representations of the solution q3 (u, τ ), in terms of space-time convolutions of the r.h.s. with the heat kernel. This allow for precise numerical estimate of this function.

Here, we are only interested in its value at τ = 1 because of the relation (S165). Furthermore, because of the expected form of the equation (S38), the combination

should take a simpler form (the difference in the factors with (S38) comes from the fact that u = √ 2v). We can thus write an integral representation for this expression, instead of q3 , which can be computed numerically with an arbitrary precision for a large number of points (∼ 500). Fitting these points with the functions we expect from the equation (S40), we find that for u > 0 The coefficients obtained from this fit are extremely stable: they do not change by more that 10 -12 when the interval or the number of points are changed, or when adding other functions to fit with (these functions then get very small coefficients < 10 -10 ). Furthermore, we find that

Therefore,

Although we first obtained this result numerically as described here, it can actually be proved from the integral representation of the l.h.s. in terms of a space-time convolution: the spatial integral can be computed using [5], and the remaining time integration reduces to the above result after some manipulations. Unfortunately, this procedure can only be carried explicitly at this order, while the numerical evaluation can be performed at higher orders. Indeed, using this procedure, we also obtained:

, (S181)

Going back to the profiles Φ n (v) using the relation (S165), these expressions coincide perfectly with the results obtained from our integral equations (S67,S68).

B. Numerical solution at arbitrary λ

The MFT equations (S160,S161) have a forward/backward structure due to the terminal condition on p (S162) and initial condition on q (S163). One can still obtain a numerical solution using the scheme described in [27], which we briefly summarize here.

1. First solve the equation for q (S160) using the initial guess p(u, τ ) = BΘ(u -Y ), corresponding to the terminal condition (S162) extended to all times.

2. Then solve the equation for p (S161) using the newly obtained function q(u, τ ).

3. Iterate the process, replacing each time either p or q by the newly obtained function. After a few iterations (∼ 3), the stability of the algorithm can be improved by replacing the functions by a linear combination of the last two, e.g.,

For instance with α = 0.75.

The Heaviside Θ function must be regularized in order to use the standard methods for solving partial differential equations. We used the following approximation

with a ∼ 100. This regularization causes a small discrepancy between the numerical solution and the exact one near the discontinuity of the function q(u, 1). This algorithm uses explicitly as inputs the values of Y (the position of the tracer) and B, which are related in this problem by the conservation relation (S159). This relation is not satisfied for arbitrary values of both Y and B. Given Y as an input, we find the corresponding value of B by performing a dichotomy, until relation (S159) is verified. The value of the parameter λ is deduced from the definition of B (S162). This gives the solution q(u, τ ) of the problem for a given position of the tracer Y , and thus the profiles Φ(v) from (S165). The plots for different models given in the main text are in excellent agreement with the solution of our main equations (S67,S68).

V. NUMERICAL SIMULATIONS

A. Symmetric exclusion process

The simulations of the SEP are performed on a periodic ring of size N , with M = ρN particles at average density ρ. The particles are initially placed uniformly at random. The jumps of the particles are implemented as follow: one picks a particle uniformly at random, along with one direction (left and right with equal probabilities). If the chosen particle has no neighbor in that direction, the jump is performed, otherwise it is rejected. In both cases, the time of the simulation is incremented by a random number picked from an exponential distribution of rate N .

We keep track of one particle (the tracer) and compute the moments of its displacement and the generalized density profiles. The averaging is performed over 10 8 simulation.

B. Kipnis-Marchioro-Presutti model

We consider a periodic lattice of N = 500 sites, each one carrying a continuous energy variable ε i > 0. Initially the energy of each site is picked independently from a Boltzmann distribution at inverse temperature ρ = 1. At a random time picked from an exponential distribution of rate N , we pick uniformly a site n of the lattice. The total energy of sites n and n + 1 is randomly redistributed between these sites with a uniform distribution. This process is repeated until the final time t is reached.

The position of the tracer (initially chosen as 0) is defined as the boundary which delimits two regions where the energy is conserved (upon subtracting the flux through the periodic boundary). In the continuous limit, this is equivalent to the definition (S123).

The generalized density profiles are averaged over 1.2 • 10 9 simulations.

C. Random-average process

By construction of the RAP, if the density of the particles is denoted by ρ and if x and t are the spatial and temporal coordinates, the observables depend only on the two rescaled coordinates z = ρx and τ = ρt. For this reason, we only consider the RAP at density ρ = 1.

The simulations are performed on a periodic ring of length L = 10000, with N = 10000 particles at positions x i . We choose a uniform probability law for the jumps of the particles. The steady state of the RAP is non-trivial, and can be written in terms of the gaps g i = x i+1 -x i between two particles (with x N +1 = x 1 ) as [28] P N,L ({g n }) ∝

Denoting G n = √ g n , this corresponds to a uniform distribution of the vector (G 1 , . . . , G N ) on the N -dimensional sphere of radius √ L. This initial condition can be easily implemented by generating N i.i.d. Gaussian random variables X i with zero mean and unit variance, and computing

The observables of interest are then averaged over 8 • 10 6 simulations.


Emna Amri

Hermann Courteille

Alexandre Benoit

Philippe Bolon

Dominique Dubucq

Gilles Poulain

Anthony Credoz

Deep learning based automatic detection of offshore oil slicks using SAR data and contextual information

Keywords: Deep Learning, Image Segmentation, Neural Networks, Remote Sensing, Oil Slick Detection, SAR Images, Wind Speed

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

In the wide era of offshore data, the detection of anthropogenic or natural oil slicks has been a longstanding challenge. Their potential impact on ecosystems makes their detection a mandatory task. [START_REF] Girard-Ardhuin | Oil slick detection by sar imagery: potential and limitation[END_REF] State of the art highlights the usefulness of remote sensing data for their detection, particularly the Synthetic Aperture Radar(SAR) technology. [START_REF] Brekke | Oil spill detection by satellite remote sensing[END_REF] Its main advantage is providing global coverage, independent of sunlight, weather and cloud coverage. Nevertheless, slick detection remains a challenge, particularly because of the variability in slick shape and extent. In addition, slicks can be easily confused with other natural sea surface phenomena. Moreover, oil slicks detection depends on the wind conditions that modulate the contrast of the radar response between the oil and the surrounding sea surface. A further challenge is the ability to quickly and efficiently process the number Further author information: E.A, H.C A.B, P.B : E-mail: emna.amri@external.totalenergies.com, (hermann.courteille, alexandre.benoit, philippe.bolon)@univ-smb.fr, Telephone: +33 (0)7 51 39 45 33. D.D, G.P, A.C : E-mail: (dominique.dubucq, gilles.poulain, anthony.credoz)@totalenergies.com, Telephone: +33 (0)6 34 49 18 34. 1 of acquired images. Since oil slick monitoring requires high resolution and high revisit frequency, automatic detection becomes necessary to assist photo-interpreters.

To address this, advanced models based on Deep Neural Networks (DNNs) have been proposed, 34 with the prospect of higher accuracy and generalization capability compared to classical approaches. [START_REF]Sar imagery for detecting sea surface slicks: Performance assessment of polarization-dependent parameters[END_REF] Different strategies can be proposed. The choice depends on the operational aims, the data, the model relevance and the postprocessing strategies. This work considers an annotated dataset of real slick monitoring scenarios coupled with expert interpretation and wind estimation. We propose the following contributions: i) a comparison of two approaches for oil slick detection (slick instance detection and segmentation v.s. semantic segmentation), ii) a further performance analysis that considers contextual factors such as wind speed and slick size. To our knowledge, this is the first study evaluating the impact of these critical factors on an extensive data collection covering different world areas.

The organization of the paper is as follows: Section 2 presents the oil slicks phenomena, state of the art, and the importance of the wind information. Section 3 reviews the considered models. Section 4 describes the SAR data pre-processing and the performance metrics. A final section discusses the results and further work.

OIL SLICKS DETECTION RELATED WORKS

Oil slicks detection from SAR data

The phenomenon of oil slicks on the sea surface is observed relatively often. There are two types of oil slicks: anthropogenic (spills) and natural (seeps). Fig. 1(a). illustrates a sample of both types of oil slicks. Oil slicks dampen the waves on the sea surface, thus reducing the sea surface roughness and the corresponding radar backscatter. [START_REF] Miegebielle | Use of remote sensing techniques to survey detect and interpret hydrocarbon seeps and spills[END_REF] As a consequence, oil slicks appear on SAR images as dark patches contrasting with the brightness of the surrounding clean sea. [START_REF] Alpers | Oil spill detection by imaging radars: Challenges and pitfalls[END_REF] However, several natural phenomena (algal blooms, currents and areas of low wind), called lookalikes, can generate similar areas of low back-scattering areas, as shown in Fig. 1(b). These lookalikes are therefore known to cause false alarms. State of the art on oil slick detection is rich, we provide a brief classification of the approaches relying on the surveys proposed by Brekke et al. [START_REF] Brekke | Oil spill detection by satellite remote sensing[END_REF] and Alpers et al.: [START_REF] Alpers | Oil spill detection by imaging radars: Challenges and pitfalls[END_REF] 1. Entirely manual inspection, where photo-interpreters are trained to detect and characterize oil slicks of different origins on SAR images. A tedious task with high resource costs. In addition, it cannot handle the number of images acquired by the sensors.

2. Conventional approaches are based on three steps: detection and isolation of all dark patterns feature extraction 8910 and dark patterns classification. [START_REF] Fingas | Review of oil spill remote sensing[END_REF] Features may be related to wind speed, geometric parameters, contrast, proximity of ships and oil installations, and detection of nearby seeps in case of seeps detection. These approaches present poor generalization behaviours, remaining accurate only in specific configurations.

3. Semi-automatic approaches, where machine learning tools such as Neural Networks (NNs) are integrated into the steps of the conventional approach in order to overcome its limitations. [START_REF] Alpers | Oil spill detection by imaging radars: Challenges and pitfalls[END_REF] 4. End to End learning approaches such as DNNs. In this paper, particular emphasis is laid on the works carried out using SAR data.

Different studies on offshore oil slick detection, such as Al-Ruzouq et al. ?? report several NNs models applied to oil slick detection and segmentation. Table . 1 shows examples of recent NNs models using Sentinel-1 images, as well as the number and size of images used. Semantic segmentation models such as Unet and object detection methods such as Mask-RCNN have already been proposed. Moreover, it is noteworthy that few annotated SAR data are exploited, which thus limits supervised learning of large models and explains the dominance of Transfer Learning (TL) strategies that build on data representation learnt from other tasks and domains. The specificity of SAR data encourages the direct training of models with such data to avoid bias.

Contextual Data

Numerous studies report the interest of contextual information to better interpret SAR images for oil slick detection. The study proposed by Brekke et al. in 2 highlights the interest in weather conditions, distance from ships, and infrastructure positions (platform, pipeline, etc.). The correlation between the visibility of oil slicks with meteorological and oceanic conditions at particular wind speeds is highlighted. They affect the slick characteristics (size, shape, etc.) while moving oil across the sea surface. [START_REF] Espedal | Oil spill detection using satellite based sar, preoperational phase a[END_REF] Wind also impacts the backscattering contrast between sea and slick areas. At medium wind speeds, oil wave dampening is visible but gradually disappears as wind speed increases. On the opposite, locally low wind speed, such as that created by the lee of an island, can generate lookalikes with low back-scatter. [START_REF] Wang | Sar image despeckling using a convolutional neural network[END_REF] During the manual annotation of the photo-interpreters, an adjustment of their assessment is made taking into account the instantaneous wind speed. According to Fingas et al. [START_REF] Fingas | Review of oil spill remote sensing[END_REF] guidance, a minimum wind speed of 1.5 m/s is required to allow detectability, and a maximum wind speed of 6-10 m/s will further eliminate the oil signature. The most accepted limits are 1.5 m/s to 10 m/s. Further work, such as La et al. [START_REF] La | Detection of convective systems through surface wind gust estimation based on sentinel-1 images: A new approach[END_REF] and Brekke et al., 2 refines the wind speed range. A summary is highlighted in Table . 2. However, it may remain specific to localize small-scale studies that do not permit generalization. These experiments report the requirement for moderate wind speeds.

CONSIDERED OIL DETECTION MODELS

In this paper, we compare two representative neural network models among instance detection and semantic segmentation: Mask Region Based Convolutional Neural Networks (Mask-RCNN) [START_REF] He | Mask r-cnn[END_REF] with object segmentation Wind speed m/s Date 3 to 7-10 2 2005 3.5 to 7 [START_REF] Garcia-Pineda | Using sar images to delineate ocean oil slicks with a texture-classifying neural network algorithm (tcnna)[END_REF] 2009 2 to 7 [START_REF] La | Detection of convective systems through surface wind gust estimation based on sentinel-1 images: A new approach[END_REF] 2018 capability and Fully Convolutional DenseNet (FC-DenseNet) for instance segmentation. [START_REF] Jégou | The one hundred layers tiramisu: Fully convolutional densenets for semantic segmentation[END_REF] While the latter cannot differentiate between object instances, in the case of a two-class problem (sea/oil), a standard connected component algorithm applied on its segmentation masks adds the detection capability of the instances. Both models can thus detect oil slicks instances and related masks. However, their complexity differs strongly and impacts their training strategies, as detailed in the following. FC-DenseNet model can be trained directly on the target data. It relies on the identification of a non-linear function f s (x, θ s ) = y with parameters θ s optimized to map a sample SAR image x to a semantic mask y that differentiates the background, including lookalikes and oil slicks. It generates segmentation masks at the same image scale. This Encoder-Decoder architecture is based on dense blocks whose internal layers process the aggregation of all the preceding features. It maximizes the reuse of features and relies on significantly fewer parameters than models like the UNet vanilla model. It can thus be trained on less data with a more limited risk of over-fitting.

Regarding the Mask-RCNN model, this multitask approach generates target object region bounding boxes with object classification scores and segmentation masks. It is based on a feature extractor f b (x, θ b ) = r that extracts the features r from the x input and feeds several specialized heads f ti (f b (x, θ b ), θ ti ) = y ti so that a given head produces a task-specific output based on a unified data representation. However, its complex structure cannot be trained directly on the target data domain when few data are available. Hence the application of transfer learning is required. As stated in previous study, [START_REF] Emna | Offshore oil slicks detection from sar images through the mask-rcnn deep learning model[END_REF] we transfer parameters from a pre-training on COCO dataset [START_REF] Lin | Microsoft COCO: common objects in context[END_REF] and consider the fine-tuning of the features extractor and region proposal head and the training of new classification and segmentation heads on the SAR data.

USED DATA AND PERFORMANCE EVALUATION 4.1 Data Acquisition And Pre-processing

This work relies on ESA's Sentinel-1 SAR data, C-band frequency (5.405 GHz) and 10m per pixel resolution. We focus on two specific areas (Namibia and Lebanon) monitored by photo-interpreters over 2015-2019. A set of 1428 images with their slick annotations is selected, and the corresponding sea surface wind estimations are provided based on the method proposed by Mouche in. [START_REF] Mouche | Sentinel-1 ocean wind fields (owi) algorithm definition[END_REF] Data pre-processing relies on the workflow described by Amri et al. in [START_REF] Emna | Offshore oil slicks detection from sar images through the mask-rcnn deep learning model[END_REF] and consists of image enhancements that facilitate slicks visibility for photo-interpreters. Regarding the wind speed information, its estimation is based on empirical Geophysical Model Functions (GMFs) such as CMOD variants. It is used to extract wind speed from SAR images using the relationship between the sea surface roughness and the SAR back-scatter, considering many parameters, including incidence angle and radar wavelength. [START_REF] Mouche | Sentinel-1 ocean wind fields (owi) algorithm definition[END_REF] The wind speed is nevertheless underestimated over the slick areas because of the wave damping effect. Therefore, for performance evaluation, we consider the average wind speed on each annotated slick region from the wind speed estimated in the neighbourhood, relying on a radius of 50 meters. The characteristics of the data crops extracted from the SAR images are summarized in Table 3. The dataset is based on 512x512 pixels sized crops that are split into training and validation sets with respective ratios of 85%/15%, while 10 other large images are kept as a test set. Those images are captured from a third region never used for training/validation. This dataset contains 214 oil slicks and was carefully selected to give a broad representation of the slicks shape and type variation. Several lookalikes phenomena are also present, such as no wind areas. The cost of such data preparation limits the number of available samples but represents real monitoring and annotation scenarios. Since oil slicks are small targets on large images, they cover less than 4% of the total pixels in the dataset. Thus, the imbalance between the number of slick pixels and the number of sea pixels that include lookalikes is highlighted. Fig. 4 and Fig. 5 illustrate the oil slicks annotation distribution with respect to the wind speed and oil slick size ranges on the accumulated (blue and green) columns for the test. As for the training and validation dataset, small to medium sized slicks are reported. Most of them were annotated at average wind speeds (3m/s), few were detected at very low speeds and no slicks could be annotated at winds above 6m/s.

Performance Metrics

Since both approaches allow instance-level detection and semantic segmentation, the following metrics are considered: (i) semantic segmentation quality metrics (pixel level) such as Intersection over Union (IoU) of the oil slicks, (ii) standard object detection metrics (instance level) counting detection from a minimum intersection with the Ground Truth (GT), and therefore missed detections, False Alarms (FA) and their rates. These measures are evaluated based on contextual information: slick size and local wind speed. For each, a categorization is established. For the wind, a compromise is considered between the speed ranges proposed in the literature. For slick size ranges, the categorization is made according to the size of the slick in the data collection. One must note that our study relies on real monitoring scenarios. Under these conditions, the photo-interpreter cannot provide accurate annotations on object boundaries due to fuzzy contours of the slicks and the annotation rhythm. Thus, object detection metrics are more relevant in this case study.

RESULTS AND EVALUATION

These results are obtained on a test data set including several slicks of different nature, sizes and acquisition data. Besides, the test images are sourced from a new acquisition area, which is not seen during the learning phase. Changing the study area leads to variations in weather conditions and, therefore, sea surface and slicks characteristics. The performance results of both approaches are reported with respect to the wind speed level and the size of the slick.

Quantitative results at the oil slick instance level are presented in Fig. 4 and Fig. 5. Similar behaviour can be noted for both models, which are able to detect slicks over the whole range of wind speed (i.e. ]0; 6] m/s) and slick size (i.e. ]0; 10 6 ] hm²). The FC-DenseNet model has a better slick instance detection rate (92.5%) than Mask-RCNN (82.7%). This result may highlight the difference in the training process of the models; FC-DenseNet is trained directly on the SAR data while TL of multimedia data tunes Mask-RCNN. However, the FA rate of FC-DenseNet model is 1.5 higher than that of the Mask-RCNN. For both models, FA areas are greater than 20 hm² and mostly correspond to windless areas. Furthermore, FA rates decrease as wind speed increases for both models.

More globally, over the entire test dataset, the amount of well detected slick pixels for FC-DenseNet (13%) is significantly better than for Mask-RCNN (7%) with slick IoU values of 0.30 and 0.06 respectively. The slick IoU values remain low, which can be explained partly by the annotation uncertainty of the slick boundaries during the manual annotation step that is shown in Fig. 3.

Qualitative results are illustrated in Fig. 3. Four oil slicks have been selected from the test set. Both model predictions are reported as well as the GT annotations. One observes that Mask-RCNN omits one annotated oil slick while FC-DenseNet finds them all. The IoU results can be illustrated in Fig. 3, where Mask-RCNN shows less overlap with the GT once the detection is valid. Regarding FC-DenseNet, the higher FA rate can be explained visually: oil slicks masks can be fragmented into several instances that do not always intersect with the raw annotation. A post-processing step and refined annotation can solve such a problem. We notice that in the case of large slicks, the detection of FC-DenseNet is fragmented, as shown in the right prediction of Fig. 3. This may be caused by the lack of large slicks in the database, as they are divided into (512*512) crops, so the network is given a partial view. If a slick covers the entire crop, its detection will be hampered by the lack of contextual information. Morphological post-processing can be performed to optimize the detection of large slicks.

CONCLUSION AND FUTURE WORKS

This paper compared instance detection and semantic segmentation approaches for oil slick detection relying on popular deep neural network models. The contribution focuses on performance analysis regarding contextual information, particularly wind speed and slick size in a real oil slick monitoring scenario. Offering readymade predictions to the photo-interpreters presents a potential change in their task. This work transforms the manual and tedious task of oil slick detection on SAR images into a task of verification and validation of the obtained predictions, thus significantly improving the process. Future work will introduce contextual information as input of the models to improve detection performance and decrease false alarm rates.

2

 2

  (a) Sample of oil spill (left) and seep (right) outlined in green. (b) Lookalikes phenomenon samples, (right) low wind area, (middle) algal bloom, (left) algae.

Figure 1 .

 1 Figure 1. Example of oil slicks (a) and similar phenomena (b) from Sentinel-1 SAR images.

  (a) Input SAR image. (b) Output of FC-DenseNet. (c) Output of Mask-RCNN.

Figure 2 .

 2 Figure 2. An example of FC-DenseNet (b) and Mask-RCNN (c) prediction on a SAR image (a) containing slicks.

Figure 3 .

 3 Figure 3. An example of FC-DenseNet (blue) and Mask-RCNN (green) prediction on Sentinel-1 images where the photointerpretation slick annotation is marked by red polygons.

Figure 4 .

 4 Figure 4. Slick instance detection performance with respect to the wind speed(m/s). Green for well detection, blue for missed detection and gray for false alarms (the right log-scale axis).

Figure 5 .

 5 Figure 5. Slick instance detection performance with respect to the size of the GT (hm²), same legend as Fig. 4

Table 1 .

 1 List of different deep learning models for oil slicks segmentation using SAR data.

	Architecture	Crop number Crop size TL
	DeepLabv3+ 3	1002	321Ö321	-
	Autoencoders 12	-	256Ö256	-
	FCN/UNet 13	713	160Ö160	-
	Mask-RCNN 4	9302	512Ö512

Table 2 .

 2 The optimal wind speed range for oil slicks detection according to the literature.

Table 3 .

 3 Considered data collections.

		Slicks	Sea including lookalikes	Total
	Crops	2208	861	3069
	Pixels ∼ 2 * 10 7	∼ 7 * 10 8	∼ 7 * 10 8
	Ratio	3.4%	96.6%	-

ACKNOWLEDGMENTS

This work has been done thanks to the HPC resources (Pangea III) from Total and GENCI-IDRIS (Grant 2020-AD011011418). We also thank NVIDIA Corporation for the donation of one NVIDIA TitanXP GPU.