Joint distribution of currents in the symmetric exclusion process

,

The Symmetric Exclusion Process (SEP) is a paradigmatic model of single-file diffusion [1,2], which has been the object of several recent and important developments [3][4][5][6][7][8].In this model, particles perform symmetric random walks in continuous time on an infinite one-dimensional lattice, with the constraint that there can be at most one particle per site.In this context, several quantities have attracted attention: (i) the integrated current through the origin Q t (defined as the number of particles which have crossed the origin from left to right, minus those from right to left, up to time t) [8][9][10][11]; (ii) the position X t of a tracer [5][6][7][12][13][14][15][16][17][18], initially placed at the origin; (iii) the generalised current J t which counts the number of particles which cross a moving boundary 1 at position x t (counted positively from left to right, and negatively from right to left) [3,4,12].This latter observable actually provides an alternative way to study the displacement of a tracer, since its position X t corresponds to the value of x t for which J t = 0, because the order of the particles is conserved [3,4,12].The statistical properties of the current Q t have been fully characterised by the computation of its cumulant generating function using Bethe ansatz [9].Concerning the position X t of a tracer, its fluctuations have first been quantified by the computation of the variance [12].The full distribution has later been computed, first in the high-density limit [19], then in the lowdensity limit [16,17,20,21], and finally at arbitrary density [3,4] by relying on tools from integrable probabilities (which give a microscopic solution).These latter works [3,4] actually provide the full cumulant generating function of the generalised current J t , from which the statistical properties of the position of the tracer is deduced.
Recently, combining microscopic and macroscopic approaches, it has been shown that, in the long time limit, all these results can be easily recovered from the solution of a simple integral equation [6,7].This equation is satisfied by generalized density profiles [5,6], which characterise the correlations between the observable under consideration (Q t , J t or X t ) and the density of particles in the SEP.On top of providing a more direct way to obtain the statistical properties of these observables, this equation constitutes a strikingly simple closure of the infinite hierarchy of equation satisfied by these generalized density profiles [6,7].These results are part of a context of intense activity around exact solutions for one-dimensional interacting particle systems [8,[22][23][24][25].
Although the individual properties of Q t , J t or X t have been characterised, the determination of their correlations remains an open question.These observables are indeed expected to be strongly correlated since, for instance, if the tracer (initially at the origin) moves to a position X t to the right, then the current Q t through the origin can only be positive.The determination of these correlations is the main goal of this article.
Here, relying only on a macroscopic description of the SEP, we show that the integral equation of [6,7] can be generalised to describe the joint correlations between the currents Q t , J t and the density of particles of the SEP in the long time limit.As a consequence, we deduce the joint statistical properties of Q t and J t , and, as a byproduct, those of Q t and X t .Importantly, this equation is completed by boundary conditions, which were derived from microscopic considerations in [5][6][7].Here, we provide a macroscopic derivation of these boundary relations, and extend them beyond the SEP to any single-file system.Furthermore, thanks to this generalisation, we give a clear physical meaning to these relations.
The article is organised as follows.We first present in Section 2 a summary of our main results, followed by a discussion of these results and their consequences in Section 2.3.We then present in Section 3 the Macroscopic Fluctuation Theory (MFT) [26][27][28], which gives a large scale description of diffusive systems, and is our starting point.In Section 4 we first illustrate on the simpler case of low density how this approach can be used to study joint properties of different observables.We give in Section 5 the derivation of our main results, which rely both on the MFT and on the inverse scattering technique [29] which has recently been applied to MFT and related problems [8,[22][23][24][25].We give in Section 6 a perturbative solution of our main equations, from which the first joint cumulants of Q t , J t and X t are obtained.We finish by several concluding remarks in Section 7.

Summary of the main results
We consider a SEP in which particles hop with rate 1.We describe the state of the system by the set of occupation numbers {η i (t)} i∈ , with η i (t) = 1 if site i is occupied at time t, and 0 otherwise.Initially, we consider that each site is filled independently with probability ρ + for i > 0 and ρ − for i ≤ 0.
The integrated current Q t counts the total number of particles that cross the origin from left to right, minus the number from right to left, up to time t.It can be written explicitly in terms of the occupation numbers by comparing the number of particles to the right of the origin at times t and 0, Similarly, the generalised current J t counts the number of particles that cross a fictitious moving boundary located at x t at time t.It can be expressed in terms of the occupation numbers, by comparing the number of particles to the right of x t at time t, and the number of particles to the right of x 0 = 0 at initial time.Explicitly, it can be written as which can equivalently be expressed in a slightly different manner by subtracting the mean density at infinity ρ + and separating the sums.This gives the definition We have chosen a specific expression for x t , such that x t ∼ t for large t since the system is diffusive.We will consider only the case ξ > 0, but the case ξ < 0 can be obtained similarly.
Our main results concern the joint cumulant generating function of the two currents Q t (1) and J t (3), in the long time limit and their correlation with the density of surrounding particles, encoded in the generalised density profiles [5] w r (λ, The profile Φ also depends on λ, ν and ξ, but we omit the dependency on these variables to simplify the notations.These profiles actually correspond to the average occupations of the sites in the so-called tilted or canonical path ensemble [30][31][32].This ensemble describes the dynamics of the system under a bias due to the exponential terms in (5), which favors the realisations in which Q t and J t differ from their expected values.The profiles w r , and thus Φ, measure the response of the density of particles to this bias.In particular, if for instance λ > 0, the exponential in (5) favors realisations in which Q t is larger than its average.In these realisations, we expect an accumulation of particles to the right of the origin, and a depletion of the left, due to the particles that have crossed the origin from left to right.This is indeed what we observe in Fig. 1, left.The same interpretation holds for J t .Furthermore, expanding (5) in powers of λ and ν, w r generates all the correlations η r (t)Q n t J m t c between the currents Q t , J t and the occupations η r (t).In particular we discuss below the lowest order correlations.
We also characterise the correlations between the currents and the initial density of particles, encoded in the initial profile which has a similar interpretation to the profile at final time (5).For instance, for λ > 0, the configurations favored by the exponential are the ones which yield a larger current Q t .At t = 0, this is realised by a fluctuation of the initial condition with more particles on the left of the origin, and less on the right.Then, the dynamics will make these particles flow through the origin, yielding a large current Q t .This is what is shown in Fig. 1, right.
We have obtained equations satisfied by the profiles Φ and Φ, from which the cumulant generating function ψξ is deduced.We first present these main equations, before giving some of their consequences on the correlations between the different observables and finally discussing the status of these results and relations with the recent works [6,7].Figure 1: Profiles Φ (left) and Φ (right) obtained from the numerical resolution of the integral equations (8,9) with the conditions (11-17) (solid red lines), compared to the numerical resolution of the MFT equations (blue points).Here, we have set ξ = 2 2, λ = 2 and ν = 5.These can be interpreted as a mean density profile, under the condition that Q t = 0.83 and J t = 0.27.The numerical schemes used to obtain these curves are described in Appendix C and Appendix D.

Equations for the correlations and the cumulants
Instead of the profiles at final time Φ and initial time Φ, we found that it is their derivatives which verify closed equations.More precisely, we define the functions with multiplicative constants a − , a 0 , a + , b − and b + to be determined.We find that the functions Ω and Ω satisfy the following integral equations where Θ is the Heaviside step function, and the kernels which involve additional parameters α and β.The multiplicative constants in the definitions (7) are determined by the following boundary conditions which involve the chemical potential µ.For the SEP, it takes the simple form We also have boundary conditions at infinity, The last constants α and β are determined by conservation equations where coincides with the single parameter identified in the SEP [9,10,33], but with two parameters λ and ν.Together, Eqs.(8)(9)(10)(11)(12)(13)(14)(15)(16) fully determine the profiles Φ and Φ.Their knowledge allows to deduce the joint cumulant generating function by using (see below)

Consequences on the observables
We only present the results for the case of a flat initial density ρ + = ρ − = ρ because they take a simpler form, but the results from the general case ρ + ̸ = ρ − can be obtained similarly.From the above equations, we recover the known results on Q t [9] and J t [4], for instance lim where erf(z) = 2 π z 0 e −x 2 dx.In addition, we obtain the joint statistical properties of these two observables, such as their covariance where erfc(z) = 1 − erf(z) is the complementary error function.This shows that Q t and J t are strongly correlated for ξ → 0, and their correlation decay when ξ → ∞, as Remarkably, these behaviours do not depend on the density ρ of the particles, but only on ξ.This is not expected to hold for a step of density ρ + ̸ = ρ − .The knowledge of the statistical properties of J t allows to deduce those of the position of a tracer X t [3,4].For instance, we recover the known variance [12], and obtain the covariance from which we deduce that these quantities are fully correlated, since This is due to the fact that, Q t = ρX t at leading order in time, which results from However, the two observables Q t and X t are not simply proportional, even in the long time limit, since for instance This means that the variance of Q t − ρX t is nonzero, since the higher order cumulants do not vanish.However, the expression of this variance is out of reach of our approach, since we can only describe the leading t with the MFT.Note that (28) is not symmetric in ρ and 1 − ρ, since following a tracer (which is a particle) breaks the particle-hole symmetry of the SEP.In particular, Eq. ( 28) vanishes faster when ρ → 1 compared to ρ → 0. In the dense limit, both ln e λQ t and ln e χ X t are of order 1 − ρ, while (28) vanishes as (1 − ρ) 3 .This shows that Q t = X t in this limit, as can be seen from the fact that they have identical cumulants in that case [9,19].On the other hand, in the dilute limit, X t = O(1/ρ) and Q t = O(1), with ln e λQ t = O(ρ).This time (28) does not vanish faster than the cumulants of the current, and thus Q t ̸ = ρX t , as illustrated from the fact that their cumulants differ [9,17] in that case.
In addition to the joint statistical properties of Q t and J t , we also obtain the profiles Φ and Φ which quantify the correlation between the density of particles and the observables, at final and initial times.At lowest orders in λ and ν, we recover the profiles obtained previously [6,7] with x = r/ t.We additionally obtain the joint profiles, such as To understand the meaning of this expression, we can rewrite it as The blue line is the result of a direct numerical simulation of the SEP with 1000 sites, 250 particles, measured at time t = 1000 and averaged over 10 7 realisations.
is the empirical covariance.This means that this profile measures the covariance between the density of particles on one hand, and the correlations between Q t and J t on the other hand.From (30), we see that it is positive for ρ < 1 2 and negative for ρ > 1  2 .This means that adding particles when ρ < 1  2 increases the correlation between Q t and J t , while it decreases it when ρ > 1  2 .This behavior is expected, since the maximal currents are reached for ρ = 1 2 .In addition, ( 30) is extremal near x = 0 − and x = ξ + .In other words, a change of the number of particles in these sectors affects strongly the correlation between Q t and J t .This is due to the fact that the particles in these regions are more likely than distant particles to cross the two "walls" x = 0 and x = ξ and thus affect both Q t and J t .Conversely, the profile (30) vanishes between 0 and ξ, indicating that adding more particles in that region does not affect the correlation between Q t and J t .This is also expected, since these particles can only cross one "wall"2 either at x = 0 or x = ξ, and can thus affect only one of these observables.This is confirmed by numerical simulations of the SEP, as shown in Fig. 2.
Beyond the perturbative expansion in λ and ν, we can also plot the profiles for finite values of the these parameters by solving numerically the integral equations (8,9).These profiles have the interpretation of mean density profiles under the condition that Q t and J t take given values [6].They are represented in Fig. 1.For instance, the plot in Fig. 1 (left) corresponds to having currents Q t and J t larger than their mean values.This is why there is an accumulation of particles to the right of the two "walls" x = 0 and x = ξ, and a depletion to the left.Conversely, for the corresponding initial profile Φ (Fig. 1, right), there is an increase of particles to the left of the origin, and a depletion to the right, so that with the diffusive time evolution, these particles will cross the origin, and ξ, to contribute to a larger value of Q t and J t .

Discussion
Integral equations.-The Wiener-Hopf integral equation obtained previously in the case of a single observable (Q t , J t or X t ) [6,7] can be recovered from ( 8) and ( 9) by setting either λ = 0 or ν = 0.This corresponds to α = 0 or β = 0 respectively in Eqs (8,9), so that the two profiles at initial and final time satisfy the same equation, as noted previously [7,8].
Note that if we set ξ = 0, we recover the equations of [6,7].This is expected since then, J t = Q t and thus Φ and Φ involve a single observable.
In the case α ̸ = 0 and β ̸ = 0, the equation for the profile at final time ( 8) is more complicated than the one obtained in the case of a single observable [6,7], as it now involves a double convolution.On the other hand, the equation for the initial profile ( 9) keeps a simpler Wiener-Hopf structure, but with a more complicated kernel which involves the solution at final time (10).
These integral equations extend the ones discovered in [6,7] in the case of a single observable (current Q t or J t , tracer position X t ).This further emphasises the key role of such strikingly simple integral equations involving partial convolutions in interacting particle systems.
Boundary conditions.-We stress that the boundary equations (11,12) hold for any single-file system, and not only for the SEP.Furthermore, these equations take a simple physical form in terms of the chemical potential µ, which can be written as in terms of the diffusion coefficient D and the mobility σ, which describe the system at large scale [26][27][28].For the SEP, D(ρ) = 1 and σ(ρ) = 2ρ(1 − ρ), which gives the expression (14) for the chemical potential of the SEP.Explicitly, (11) states that the chemical potential in the system is discontinuous at x = 0 and x = ξ, with a discontinuity given by the parameters λ and ν of the joint generating function (4).From a physical point of view, this can be understood as follows.The parameters λ and ν play the role of conjugate variables (in the sense of thermodynamics) to the integrated currents Q t and J t , which count particles.The conjugate quantity to the particle number being the chemical potential, it is expected that λ and ν are related to the chemical potential.Finally, λ and ν have an effect on the density of particles.For instance, when λ > 0, the exponentials in the definitions of the profiles (5) give more weight to the realisations in which Q t > 0, and thus we expect an increase of the number of particles to the right of the origin.This results from a higher chemical potential to the right of the origin compared to the left, as described by (11).
Remarkably, the equations (11,12) obtained for the case of the currents Q t and J t can be extended to other observables.For instance, it has been recently shown that the current Q t in a single-file model can be mapped onto the position of a tracer in a dual single-file model [34].Under this mapping, the relations (11,12) become, for the new system (see Appendix A) where P is the pressure3 and Φ is now the long time limit of the correlation between the position X t of the tracer and the density in the reference frame of the tracer, Finally, in the case of the SEP, these relations (11,12) and (33) reduce to the ones obtained from microscopic considerations in [5][6][7].The precise relation with the corresponding microscopic equations is given in Appendix B.

Hydrodynamic description using macroscopic fluctuation theory
The Macroscopic Fluctuation Theory (MFT) gives an effective description at large scales of a diffusive system [26][27][28].Introducing a scaling factor T , which corresponds to the large observation time, the density of particles is defined as In the limit T → ∞, this density converges to a continuous stochastic function ρ(x, t).The probability of observing an initial profile ρ(x, 0) evolves to another profile ρ(x, 1) at time t = 1 (corresponding to the large time T in the SEP) takes a large deviation form [10] where H is a conjugate field, and S is the MFT action For the SEP, D(ρ) = 1 and σ(ρ) = 2ρ(1 − ρ).This result was first proved for the SEP [35], and was later extended to arbitrary 1D diffusive systems, which can be described by other transport coefficients D(ρ) and σ(ρ) [26][27][28].See for instance [34] for a list of models, and their corresponding coefficients.Initially, the system is described by the random density ρ(x, 0), which fluctuates around a given density ρ 0 (x), of distribution [10] For the moment, we consider an arbitrary initial condition ρ 0 (x), which we will later specify to the case ρ 0 (x) = ρ + Θ(x)+ρ − Θ(−x).Microscopically, it corresponds to picking independently, for each site i of the SEP, an occupation number η i (0) = 1 with probability ρ 0 (i/ T ).In the continuous limit, this becomes (38).
The two currents Q t (1) and J t (3) can be expressed in terms of the density ρ(x, t) (35) as Within this formalism, the joint moment generating function of Q T and J T reads In the long time limit T → ∞, these integrals can be evaluated by a saddle point method, which yields where we have denoted (q, p) the saddle point of (ρ, H).It can be determined by minimising the terms in the exponential, which yields the MFT equations completed by the boundary conditions The expression of the cumulant generating function can be simplified by taking a derivative with respect to λ or ν, and using that the saddle point solution (q, p) is the minimum of the action: These are standard relations in the context of large deviations [36,37].
The profiles w r (5) can be obtained from the same procedure, for instance, Performing again the saddle point estimate, we obtain, Similarly, the correlation with the initial occupations wr (6) reads wr (λ, ν, T ) ≃ The MFT profile q at initial and final time actually coincides with the correlations w r and wr in the long time limit, as shown in [5,6].Furthermore, the joint cumulant generating function ψξ is fully determined by the knowledge of the profile q, thanks to (45).Our goal is thus to determine these profiles.

The example of the low density limit
The MFT equations for the SEP (42,43) being rather complicated, we first focus on the simpler case of the low density limit.In this limit, the SEP becomes equivalent to a model of reflecting Brownian particles on the real line [3].The MFT equations reduce to These equations can be reduced to diffusion equations by the Cole-Hopf transform P = e p and Q = qe −p [10,21], so that The initial and final conditions (44) become We straightforwardly obtain the solution q = QP, and thus the profiles both at initial and final times, where we have assumed that ξ > 0 for Φ.The expression is similar in the case ξ < 0. The cumulant generating function can be obtained from the relations (45), which gives Integrating with the initial value ψ(0, 0) = 0, we get Note that this expression is compatible with the very recent study [38].
Finally, the density profiles of the SEP assume simple explicit forms in the low density limit.

Derivation of the main equations
We now address the case of arbitrary density of the SEP.To obtain the equations satisfied by the initial and final profiles Φ and Φ, we will rely on the inverse scattering approach which has recently been applied to solve systems of equations related to (42,43), in the context of the KPZ equation or MFT [8,[22][23][24][25]39].As we will see below, this formalism is powerful to obtain the bulk equations for Φ and Φ, but introduces unknown constants which can be tricky to determine.Here, we will obtain these constants by making use of boundary conditions which are deduced from the MFT equations (42-44).

Boundary conditions
We first derive the boundary conditions (11)(12)(13), which are direct consequence of the MFT equations (42)(43)(44).These equations will take a simple form, in terms of physical quantities.The equation satisfied by p (43) is an antidiffusion, with no singularity in the r.h.s. for t < 1, except a discontinuity for q at t = 0. Therefore, the solution p(x, 0) is a smooth function of x, and in particular at x = 0.The boundary conditions are thus straightforwardly deduced from the initial condition (44), which takes the from, where we have introduced the chemical potential µ(ρ), defined by (32).Evaluating (58) at x = 0 + and x = 0 − , and taking the difference and using the continuity of p(x, 0) at x = 0, we get the first boundary condition for q(x, 0) ≡ Φ(x) Similarly, writing the continuity of the first derivative of p(x, 0) at x = 0, we deduce from (58) For ρ 0 (x) = ρ + Θ(x) + ρ − Θ(−x), these relations become (13).
To obtain the conditions at final time, we rely on a time-reversal mapping, which extends the time-reversal symmetry discussed in [10] for the case of the current Q t .In that work, the MFT action was found to be invariant under time reversal symmetry ρ(x, t) → ρ(x, 1 − t) and j(x, t) → − j(x, 1 − t), with a density ρ and a current j satisfying the conservation relation ∂ t ρ + ∂ x j = 0.At the saddle point of the MFT action, the density becomes q and the current becomes j = −D(q)∂ x q + σ(q)∂ x p [10].The time reversal symmetry then becomes q(x, t) → q(x, Here, we do not have the time-reversal symmetry, because at final time the currents are measured at positions 0 and ξ, which are different from the initial position 0. Nevertheless, we define two new fields q and p by Integrating the second relation gives with c a constant.Inserting these relations into the MFT equations (42,43), we find that q and p obey the same equations: This was already noticed in [40].The initial and final conditions (44) become p(x, 0) = µ(q(x, 0)) These conditions are different from the original ones (44), and they are the source of the breaking of time-reversal symmetry. 4We can however still use a similar argument as we used above at t = 0.The conjugate field p obeys an antidiffusion equation, which is not singular for t < 1. Therefore p(x, 0) is smooth.From (64) left, this straightforwardly yields the conditions for q(x, 0) = q(x, 1) ≡ Φ(x), and from the derivative, These are the relations (11,12) announced above.Note that these equations for Φ hold for any initial density profile ρ 0 .
Important remark: We have derived the boundary conditions for Φ in the case of an annealed initial condition.One could also consider a quenched initial condition, which corresponds to q(x, 0) = ρ 0 (x).In this case, the mapping (61) can still be performed.One obtains the same MFT equations (42,43), but with the initial and final conditions p(x, 0) = µ(q(x, 0)) The second relation involves the unknown function p(x, 0), which is smooth, but the first relation is identical to the annealed case.The same argument as above applies, and the boundary conditions (65,66) still hold in the quenched case.

Mapping to the AKNS equations
We adapt the inverse scattering approach that was applied to the case of the integrated current Q t in the SEP in [8] to the case of the joint distribution of Q t and J t .The first step is to introduce the new functions Under this transformation, the MFT equations for the SEP (42,43) become the AKNS equations [41] These equations are integrable and can be solved using the inverse scattering transform [29].
Before entering the resolution in more details, let us study the initial and final conditions for u and v. From the conditions on p and q (44) and the transformation (68), we obtain From now on, we consider the case of a step initial density ρ 0 (x) = ρ + Θ(x) + ρ − Θ(−x).In Eq. ( 70), the term ∂ x ρ 0 thus gives another δ(x) term, but with an unknown prefactor, because ρ 0 is discontinuous at 0. Even in the case of a constant density ρ + = ρ − , the prefactor of the remaining δ function is unknown, because q is discontinuous at x = 0. Therefore, we can only write with an unknown constant c 0 .Similarly, for v(x, t = 1) at final time, we get with two other unknown constants c 1 and c 2 , coming from the fact that the term in the exponential is not well defined since p and q are discontinuous at x = 0 and x = ξ.Actually, we can get rid of one of these unknown constants by using the invariance of the AKNS equations (69) under the transformation u(x, t) → u(x, t)/K and v(x, t) → K v(x, t).Choosing K = c 0 , we have the initial and final conditions with α = c 1 c 0 and β = c 2 c 0 .We will see below how we can determine the constants α and β.The simplicity of these conditions will allow for an explicit solution of the AKNS equations at initial and final times.Furthermore, this solution will yield the desired equations for the profiles, since ∂ x p(x, 1) is a sum of δ functions as seen from (44), ) , for x < 0 , a 0 ∂ x q(x, 1) , for 0 < x < ξ , a + ∂ x q(x, 1) , for x > ξ , (74) with different proportionality constants a − , a 0 and a + for each domain x < 0, 0 < x < ξ and x > ξ because p(x, 1) is discontinuous at both x = 0 and x = ξ, hence the value of the exponential differs in each interval (by an unknown factor, since q(x, 1) is also discontinuous at these points).Similarly, We can simplify this equation in the following way.We take the derivative of the boundary condition (44) at t = 0, which gives, with new constants c 3 and c 4 (which we will not need to determine).Indeed, combining with (75), we get with two different proportionality constants b − and b + for x < 0 and x > 0.
To summarize, the solutions u(x, t) and v(x, t) of the AKNS equations are directly related to the derivative of the profiles at initial (48) and final times (47), (78) with constants a − , a 0 , a + , b − and b + that will be determined by the boundary conditions derived in Section 5.1.The other relations are given by (73), with the constants α and β that remain to be determined.

Solution using the scattering technique
Our goal is now to obtain integral equations verified by Ω and Ω.To solve the AKNS equations we rely on the standard approach [29], recently used in [8,22,24], and introduce the auxiliary linear problem for the two-component vector Ψ, The compatibility condition between the first two equations, 69).The idea is therefore to solve the simpler linear problem (79), and deduce the solution for u(x, t) and v(x, t).Since ∂ x q → 0 and ∂ x p → 0 for x → ±∞, u(x, t) and v(x, t) decay to 0 at ±∞.The matrix U then becomes diagonal at ±∞, and therefore Ψ is a superposition of plane waves in this limit.We introduce two independent solutions φ and φ, defined by their behaviour at −∞, where we have placed the factors e ±2k 2 t so that φ and φ satisfy the time evolution equation (79) at −∞.For x → +∞, we can write the solution as the superposition of the same two plane waves, This defines a scattering problem, in which plane waves at −∞ are scattered by the potentials u(x, t) and v(x, t) into a superposition of plane waves at +∞.The coefficients a, ā, b, b are called the scattering amplitudes.All the information on the functions u(x, t) and v(x, t) are encoded in the scattering amplitudes, so that u(x, t) and v(x, t) can be reconstructed from a, ā, b, b.This is called the inverse scattering procedure, and is quite complicated to do in practice.Here, we will follow a different route, used in [8, 22-25, 39, 42]: we will determine the scattering amplitudes at initial and final times in terms of the functions u(x, 1) = Ω(x) and v(x, 0) = Ω(x), and relate them using the time evolution (79) to obtain integral equations satisfied by these functions.Indeed, one strength of the scattering approach is that it transforms the complicated time evolution of the AKNS equations (69) into a very simple time dependence of the scattering amplitudes.Their time evolution can be computed using the matrix V (+∞) = Diag(2k 2 , −2k 2 ) in (79) at +∞, which directly gives and thus There only remains to determine the scattering amplitudes at t = 0 and t = 1.For this, we solve the spatial equation involving the matrix U (79).Let us first write this equation at t = 0.For the second component of φ = (φ 1 φ 2 ) T , we get and the same equation holds for φ2 .Integrating this equation with the boundary conditions at −∞ (80), we obtain Using these expressions in the equations for the first components φ 1 and φ1 (79) yields Integrating with the boundary conditions (80), we obtain From these expressions, we deduce the expressions at x = 0, Combining the results (85,87,88) with the asymptotic behaviour (81), we deduce the scattering amplitudes The amplitudes a(k, 0) and ā(k, 0) can also be determined, but we will not need them in the following, so we do not write their expressions explicitly.We can proceed similarly at final time t = 1, this time starting with the equations for the first components, as it is the one that involves the δ functions, with again the same equation for φ2 .Integrating with the asymptotic at −∞ (80) yields Inserting these expressions into the equations for φ 2 and φ2 , we get Integrating with the asymptotic behaviour (80), we obtain We therefore deduce From the solutions (92a,94), we can read the asymptotic behaviours (81) and deduce the scattering amplitudes Again, a and ā can be obtained similarly but we will not need them here.The last step is to relate the scattering amplitudes b and b at initial time (90a) with the ones at final time t = 1 (99,100) by using the time evolution (83a).This gives the following equations for Ω and Ω: We can obtain equations in real space by taking the inverse Fourier transform.More precisely, multiplying (101) by e 2ik x /π and integrating over k yields the equation for Ω (8).Similarly, multiplying (102) by e −2ik x /π and integrating over k, we obtain the equation for Ω (9).These integral equations (8,9) clearly admit a unique solution when α = β = 0, given by Ω(x) = K(x) and Ω(x) = 0.One can then use this starting point to look for a perturbative solution for small α and β, as it is done in Section 6 below.This leads us to expect that these equations admit a unique solution, at least for small α and β.
This concludes our derivation of the integral equations (8,9).There only remains to determine the constants α and β.

Determination of the remaining constants α and β
We now turn to the determination of the last constants α and β which appear in the integral equations (8,9).A first equation can be obtained from the conservation of the number of particles in the SEP between initial and final time, i.e., This equation can be derived from the MFT equation ( 42), by integration on x from −∞ to +∞, which yields ∞ −∞ ∂ t q(x, t)dx = 0, since ∂ x p and ∂ x q decay to 0 at ±∞.The second equation can be determined by following the approach used in [8], which relies on the scattering formalism.The scattering amplitudes defined in (81) can be equivalently defined by regrouping the two vectors φ and φ in a single matrix, so that where the matrix M (x; k, t) satisfies with the matrices U and V given in (79).Remarkably, the spatial equation for M can be explicitly solved when k = 0 by using the specific form of the functions u and v in terms of p and q (68).The solution is given in the Supplemental Material of Ref. [8], and reads where K is the constant we introduced above Eq.( 73) to get rid of the constant in front of the δ function in the initial condition for u.Using the asymptotic behaviours q(x, t) where we have denoted C = e − ∞ −∞ q∂ x p .Using these expressions in (104), we obtain a simple expression for the product of the diagonal elements at k = 0, with ω the single parameter identified for the SEP in [9,33], still with two parameters λ and ν.Using the expressions of b and b at t = 0 (90a), and the bulk equation for Ω (102), this last equation yields Note that the same equation can be obtained from the expressions at t = 1 (99,100).Comparing with the definition of the kernel K (10), we notice that this equation can also be written in the compact form With this last derivation, we now have all the equations needed to determine the profiles Φ and Φ and thus deduce the cumulants from (45).

Perturbative solution for the first joint cumulants
We do not have an explicit solution of the equation for Ω (8), so we will rely on a perturbative solution in λ and ν.

For the currents
The equations for Ω (8) and Ω (9) only involve the parameters α and β.We thus first write the solutions of these equations perturbatively in α and β, and in a second step express them in terms of λ and ν by using relations ( 11)-( 16).We denote the expansions of Ω and Ω as Inserting these expansions into the integral equations (8,9), we obtain To deduce Φ and Φ, we integrate Ω and Ω (7) with respect to x, with the boundary conditions at infinity (15), (117) For the above expressions, these integrals can be computed using the tables in [43].We will also need the integral of Φ and Φ in Eqs.(16,17), which correspond to double integrals of Ω.This is not convenient to compute in practice, and it is more practical to use integration by parts which can now be computed using the tables in [43].
Next, we expand all the parameters in powers of λ and ν, with Z ∈ {α, β, a − , a 0 , a + , b − , b + , d 0 }.Inserting these expansions into the boundary conditions at x = 0 and x = ξ (11,13) and into the conservation equations (17,16), we obtain the coefficients of these expansions up to order 4 included for α and β and up to order 3 included for b + , b − , d 0 , 1 a − , 1 a 0 and 1 a + in the case ρ + = ρ − = ρ.This difference of orders come from the fact that α and β begin at order 2 in λ and ν, Figure 3: Profile Φ at final time obtained from the numerical solution of the integral equations (8,9) with the conditions (11-17) (dashed red lines), compared to the perturbative solution (121) up to order 4 in λ and ν (solid blue line).Remarkably, the perturbative solution is in good agreement even for values λ = ν = 1 which are not small (left), but ultimately differs from the correct solution when λ or ν is increased (right).
therefore so does Ω, while Ω already has non zero terms at order 0. On the other hand, Φ and Φ have terms at first order in λ and ν.This is why the expansions of b + and b − begin at order 1, while those of a + , a 0 and a − have terms in 1/λ or 1/ν.Consequently, we obtain the lowest orders of the profiles Φ and Φ.For instance, where T(z, a) is Owen's T-function, defined by [43] T We compare this perturbative expression to the result obtained from the numerical resolution of the integral equations (8,9) in Fig. 3.Although the perturbative solution is expected to be in good agreement with the actual solution for values of λ and ν which are small, we observe a reasonnable agreement up to λ = ν = 1.To accurately describe larger values of λ or ν, one should include higher orders the perturbation series.
From the expressions of Φ and Φ, we deduce ψξ from (19), which yields One can check that, for ν = 0, we recover the first orders5 of the cumulant generating function of [9], for λ = 0 it gives the first orders of the one obtained in [3], and for ξ = 0, since J t = Q t , we recover the first orders of the one for Q t [9], evaluated at λ + ν.Taking derivatives of this expression with respect to λ and ν, we obtain the different joint cumulants The first cumulants are written in Section 2.

For the current/tracer correlations
The distribution of the position of a tracer can be obtained from the distribution of the current J t [3,4].The idea is that, since the particles remain in the same order, the number of particles to the right of the tracer is conserved, the tracer is located at the position X t such that J t (X t ) = 0.This relation is not quite exact, since there could be several values for which J t (x) = 0. Actually, X t corresponds to the smallest of these values.However, in the long time limit, this indeterminacy becomes a subdominant correction, and this relation becomes exact at leading order in t.This implies that (X t = x) = (J t (x) = 0).We can directly extend this relation to the joint distribution of Q t and X t , We have computed the joint cumulant generating function of the two currents, which implies e λQ t +νJ t (x t ) = q j e λq+ν j (Q t = q & J t (x) = j) ≃ t→∞ e t ψξ (λ,ν) .( 126) We can take the inverse Laplace transform in ν, which can be evaluated by a saddle point approximation for large t, which gives a mixed distribution/generating function where ϕ ξ is given by the Legendre transform Using the relation between J t and X t (125), we deduce We can obtain the joint cumulant generating function of the current Q t and the position X t by another Legendre transform, This procedure extends the one of [3,4] to the joint distribution of Q t and X t .It can be carried out explicitly starting from the expression of ψξ at lowest orders (123), and yields This directly gives the first joint cumulants of Q t and X t , which are given in Section 2. In particular, setting χ = −ρλ, this gives the generating function of Remarkably, there is no term in λ 2 , indicating that the variance of Q t − ρX t is smaller than t for large t, indicating strong correlations between the current Q t and the positions X t of the tracer.However, this does not indicate that Q t and X t are proportional, since the higher order cumulants grow as t.This indicates that Var(Q t − ρX t ) is nonzero, but grows slower than t.

Conclusion
We have studied the joint distribution of the current Q t through the origin and the current J t through a moving boundary in the SEP, as well as their correlations with the density of particles.These correlations are described by generalised density profiles.We have obtained integral equations satisfied by these generalised density profiles.These integral equations extend the ones discovered in [6,7] in the case of a single observable (current Q t or J t , tracer position X t ).This further emphasises the key role of such strikingly simple integral equations involving partial convolutions in interacting particle systems.In the case of a single observable, the integral equations naturally obtained are bilinear, but surprisingly they are equivalent to linear equations at the expense of introducing analytic continuations [6,7].An important open question is whether the equation obtained here, which is trilinear, can be reduced to such linear equations (for which an explicit solution can be obtained).
We have also obtained simple boundary conditions for the generalised density profiles.These boundary conditions take a simple physical form, in terms of the chemical potential, and can be applied to any model of single-file diffusion.This extends the relations that have been obtained for the SEP from microscopic considerations [5][6][7].
As a consequence of these equations, we have characterised the joint statistics of the current through the origin Q t and the position of a tracer X t , initially at the origin.These variables are strongly correlated, and even become equal in the high density limit.
This work opens the way to the study of more than two observables, such as multiple currents or tracers, in the SEP and other models of single-file systems.

B Relation between the physical boundary conditions and the microscopic equations for the SEP
In the case of the SEP, other boundary conditions have been obtained for the different observables considered here (individually), from microscopic considerations [5][6][7].Note that they have been obtained with a different choice for the time scale, with D(ρ) = 1 2 and σ(ρ) = ρ(1 − ρ).Here, we show that they are equivalent to the physical boundary conditions (11,12,33) obtained in this article.
For the case of the current Q t , they take the form [6] Φ(0 where Φ has been defined in [6] with a slightly different scaling which is indeed Eq. ( 12).
In the case of a tracer at position X t , different relations have been obtained This is identical to the case of the current (B.4), hence it yields again (33), right.

C Numerical resolution of the integral equation
To solve the integral equations (8,9), we discretize the integrals by using the trapezoidal rule .By doing so, the integral equations (8,9) become a finite system of nonlinear equations in Ω(−L + i δx) and Ω(−L + i δx) with 0 ≤ i ≤ N = ⌊2L/δx⌋.This system can be solved using standard gradient descent algorithms.
The profiles Φ and Φ are then obtained from these solutions by discrete integration of (7).The parameters a − , a 0 , a + and b − , b + are determined from the boundary conditions (11,12,13).This procedure gives the profiles Φ and Φ, given α and β as input.Performing a final gradient descent, α and β are determined from the parameters λ, ρ − and ρ + from (16,17).

D Numerical resolution of the MFT equation
The coupled MFT equations have a forward/backward structure: p obeys an antidiffusion equation (43) with a final condition (44) (left); while q obeys a diffusion equation (42) with an initial condition (44) (right).To solve this system, we use the standard scheme described for instance in [11].which we briefly summarize here.
2. Then, use the resulting solution q(x, t) to solve the equation for p (43).
3. Iterate the process, by replacing at each step either p or q by the newly obtained function.
After a few iterations (∼ 3), it is usually helpful to replace the functions by a linear combination of the last two, e.g., q(x, t) ←− αq new (x, t) + (1 − α)q old (x, t) , (D.1) to avoid oscillating between two solutions.For instance, we use α = 0.75.
4. Repeat until the difference between the last two solutions is small enough, which can be measured for instance by [q new (x, t) − q old (x, t)] 2 dxdt .(D.2) To use standard methods for partial differential equations, we regularise the Heaviside Θ function by with a ∼ 150.This regularization causes a small discrepancy between the numerical solution and the exact one near the discontinuity of the function q(x, 1).

5. 2 18 6 7 24 A Mapping the boundary conditions for other observables 24 B 25 C Numerical resolution of the integral equation 26 D
Bulk equations 14 5.2.1 Mapping to the AKNS equations 14 5.2.2 Solution using the scattering technique 15 5.3 Determination of the remaining constants α and β Perturbative solution for the first joint cumulants 19 6.1 For the currents 20 6.2 For the current/tracer correlations 22 Conclusion Relation between the physical boundary conditions and the microscopic equations for the SEP

Figure 2 :
Figure2: Correlation function 〈η r Q t J t 〉 as a function of r/ t, for a mean density of the SEP ρ = 0.25, and for ξ = 1.5.The solid red line is the analytical result(30).The blue line is the result of a direct numerical simulation of the SEP with 1000 sites, 250 particles, measured at time t = 1000 and averaged over 10 7 realisations.