Correlations in nonequilibrium diffusive systems
Abstract
ACKNOWLEDGMENTS This work is part of the Project of I+D+i Ref. No. PID2020-113681GB-I00, financed by MICIN/AEI/10.13039/501100011033 and FEDER “A way to make Europe.”
Full text
PHYSICAL REVIEW E 106, 024107 (2022) Correlations in nonequilibrium diffusive systems P. L. Garrido * Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, E-18071 Granada, Spain (Received 23 May 2022; accepted 18 July 2022; published 8 August 2022) We study the behavior of stationary nonequilibrium two-body correlation functions for diffusive systems with equilibrium reference states (DSe). We describe a DSe at the mesoscopic level by Mlocally conserved continuum fields that evolve through coupled Langevin equations with white noises. The dynamic is designed such that the system may reach equilibrium states for a set of boundary conditions. In this form, we make the system driven to a nonequilibrium stationary state by changing the equilibrium boundary conditions. We decompose the correlations in a known local equilibrium part and another one that contains the nonequilibrium behavior and that we call correlation’s excess ¯ C(x,z). We formally derive the differential equations for ¯ C.Tosolvethem order by order, we define a perturbative expansion around the equilibrium state. We show that the ¯ C’s first-order expansion, ¯ C(1), is always zero for the unique field case, M=1. Moreover, ¯ C(1) is always long range or zero when M>1. We obtain the surprising result that their associated fluctuations, the space integrals of ¯ C(1),are always zero. Therefore, fluctuations are dominated by local equilibrium up to second order in the perturbative expansion around the equilibrium. We derive the behaviors of ¯ C(1) in real space for dimensions d=1and2 explicitly. Finally, we derive the two first perturbative orders of the correlation’s excess for a generic M=2 case and a hydrodynamic model. DOI: 10.1103/PhysRevE.106.024107 I. INTRODUCTION Particle systems are characterized by the dynamics they follow: classical or quantum for material particles, stochastic rules for models in ecology, biology, etc. Moreover, boundary conditions are an essential part of dynamics because they determine the values that some variables must take in some spatial regions. Today, we can establish from first principles, theories, and observations the dynamical rules of a given system with reasonable precision. However, to extract valuable information for our understanding, we need to solve coupled ordinary nonlinear differential equations, partial differential equations, or stochastic equations with many degrees of freedom. Our mathematical tools are minimal for this enormous task. However, we managed to get an idea of the properties of the system by simplifying the original dynamics by focusing on the aspects that we consider relevant to studying a particular observed phenomenon. The most common strategy is to adapt the modeling of the system to our mathematical knowledge. That allows us to use the tools that we master to extract some answers from those complex equations. This natural-looking scheme has some drawbacks. In our opinion, the most relevant is the robustness of the chosen model, whether or not small changes in the dynamic rule imply proportionally small changes in the observed behavior. This problem is far from trivial, but it is generally neglected because we know many relevant physical situations in which the models we build are robust by construction. For example, we know that most details of the structure of a molecule and the interaction between *[email protected] them are irrelevant to describing equilibrium macroscopic properties of a system as the equation of state. Moreover, for equilibrium systems we have thermodynamics and the ensemble theory that help us to design simple microscopic models that contain the phenomena we want to characterize with detail. In conclusion, we are in a controlled environment in many developed theories where we have models that are reasonably simple and typically robust. Let us mention as a counterexample that there are very relevant equilibrium systems, such as water, where we do not know a simple model that contains all the rich set of properties and phases [1] that has been observed. Nature is far from an equilibrium state. There are currents and flows of particles and energy, unbalanced chemical reactions, births, and deaths. Dynamic details and boundary conditions frequently determine the system’s overall qualitative behavior. Therefore, the modeling of these systems becomes a very subtle issue, and robustness is always under deep scrutiny. Fortunately, there are cases in which we have successfully managed all those issues. For instance, after centuries of observation, experiments, and theories, we derived a successful macroscopic theory as the Navier-Stokes equations for fluids. They have been the starting point to understanding many exciting phenomena associated with them such as turbulence and convection [2]. Moreover, many efforts have been made in ecology to determine the basic principles and build resilient models [3]. In other relevant cases, as in the evolutionary potential games, there are profound developments in models and methods for understanding the behavior of a set of players in a general sense [4]. In recent years, we have been interested in looking for a common theoretical framework that permits us to model different systems from diverse disciplines, each with its own 2470-0045/2022/106(2)/024107(28) 024107-1 ©2022 American Physical Society
P. L. GARRIDO PHYSICAL REVIEW E 106, 024107 (2022) particular dynamic rules. That lets us look for generic properties that can be of common interest. The first step in this direction was the Onsager-Machlup’s theory for irreversible processes [5] where a Markovian mechanism is proposed to explain how the thermodynamic variables relax and fluctuate toward and around their equilibrium value. They assume that the macroscopic variables evolve by a Langevin equation where its deterministic evolution is proportional to the causes that provoke it and call them thermodynamic forces. For instance, for fluids, the local heat current is proportional to the local temperature gradient (Fourier’s law), or the local particle current is proportional to the chemical potential gradient (Fick’s law). Moreover, the stochastic process is a white noise design, so the model fulfills the fluctuations at equilibrium. The Onsager and Machlup idea was recently extended by Bertini et al. [6] to nonequilibrium systems where the time reversibility is lost, which is typical in equilibrium. They developed the macroscopic fluctuating theory (MFT) in the context of diffusive systems because we know rigorous results about hydrodynamic limits and large deviation properties. Many MFT ideas were already developed for systems with a discreet number of degrees of freedom [7], and they are easily generalized to different models with or without local conserved quantities [8]. It is, in our opinion, the natural context to develop theoretical tools that help us understand the complex behavior of nonequilibrium systems. One exciting object in MFT is the quasipotential that defines the stationary measure in the weak noise limit. It is the nonequilibrium equivalent to the thermodynamic potential for systems at equilibrium. The quasipotential has been derived for some one-dimensional systems [9,10]. Also, there is some algebraic method that may help in getting them as the solution from a Hamilton-Jacobi equation [11]. The quasipotentials have a highly complex structure with a nonlocal behavior that strongly depends on the boundary conditions. That makes it very difficult to find regularities and generic behaviors to build, if possible, complete nonequilibrium thermodynamics beyond the one based on local-equilibrium assumptions [12]. Therefore, it is convenient to get some more insight into the system’s behavior by studying the correlations. We know that the correlations are just the inverse of the kernel coming from the second-order expansion of the quasipotential around the stationary state. Moreover, they contain precious physical information about the system’s physical structure. Correlations have been extensively studied in fluids by experiments and theories. Here, we can mention the fluctuating hydrodynamics that is an MFT [13]. Fluctuating hydrodynamics is built following an Onsager-Machlup’s type of assumption by adding a local equilibrium white noise to the deterministic NavierStokes equations. In this context, we may highlight a couple of classical works by Tremblay et al. [14] and by Mansour et al. [15] where they deeply study the correlations for fluid by linearizing the Navier-Stokes equations in different situations and approximations. Inspired by the classical works in fluids, we study the two-body correlations in a generic nonequilibrium model with three main properties: (1) the system is described by Mfields that are locally conserved by the dynamics, (2) the local currents are proportional to the local field’s gradients, and (3) the equilibrium state may be reached by the system for a given set of external parameters. We call these systems DSe: diffusive systems with reference equilibrium states. We will use the last property to have a reasonable definition for the noise term and, later, to make a perturbative expansion around the equilibrium to get precise results. Section II presents the model definition through the Langevin equations and its connection with the reference equilibrium state. We also point out the properties we will assume in the paper, for instance, an unique locally stable stationary state. Section III obtains the partial differential equations for the two-body equal time correlation functions from the Hamilton-Jacobi equation for the quasipotential. We also decompose the correlations in a local equilibrium contribution and a correlation excess that carries the nonequilibrium structure because it is equal to zero at equilibrium. Section III is devoted to extracting some general property by doing a perturbative expansion of the correlation excess around the equilibrium. For example, we find that all DSe systems with only one field have a zero first-order correction. We also find that, in general, all DSe with parallel plates as boundary conditions have their field fluctuations (integrals over the space of the two-body correlations) equal to zero at first order in the expansion around the equilibrium, although their correlation excess to such order being nonzero. In Sec. V, we focus on studying the basic correlation function F, which is the common part for the correlation excess at first order in the perturbation for any model. We study its behavior numerically in one dimension in real space after a nontrivial transformation. We also look at dimensions greater than one in the thermodynamic limit but near a system’s boundary. We see the rich power-law behaviors depending on how we do the long-distance limits. In Sec. VI, we show the first-order perturbation correlation excess in the case of two fields in dimensions one and two. Finally, Sec. VII is devoted to getting the correlation excess up to second order in the perturbation expansion for a two-dimensional particle model whose hydrodynamic equations have been derived recently [16]. Some comments, most of the detailed computations, and the math relations we have derived to get the results shown in the central part of the paper have been left to the six Appendixes. II. THE MODEL Let us define a mesoscopic system defined by Mconserved real fields φα(x,t), α=1,...,Min a d-dimensional region x∈⊂Rd. The fields evolve by the Langevin equation: ∂tφα(x,t)+∇Jα(x,t)=0.(1) Jαis the local vector current associated to the φαfield that it is composed by a deterministic part, JD α, and a fluctuating one, JR α: Jα(x,t)=JD α(x,t)+JR α(x,t).(2) We study in this paper diffusive systems with equilibrium reference states (DSe). That is, we impose two conditions on the form of the currents: (1) JDshould be linear combinations of the field’s gradients (diffusive system) and (2) it should describe an equilibrium system with the appropriate boundary conditions (equilibrium reference state). We will assume in this paper only spatially uniform equilibrium reference states, 024107-2
CORRELATIONS IN NONEQUILIBRIUM DIFFUSIVE SYSTEMS PHYSICAL REVIEW E 106, 024107 (2022) and we will not consider the action of external fields, like gravity, on the system. We may think of this model as the linear approximation around a given stationary state of a much more complex nonequilibrium conserved model as, for instance, the fluctuating hydrodynamics [13] that, as we know, contains equilibrium states as a part of its description. In this class of models, nonequilibrium stationary states are built by changing the boundary conditions without introducing any other external effect. Therefore, DSe’s currents have the form JD α(x)= β gα,β (φ(x))∇φβ(x)(3) and JR α,i(x,t)= d j=1 M β=1 σα,i;β,j(φ(x,t))ψβ,j(x,t)i=1,...,d. (4) All the sums over Greek symbols run from 1 to M(the number of fields), and the ones with Latin symbols from 1 to d(the spatial dimension). φ(x,t)≡{φα(x,t)}M α=1,ψα,i(x,t)isan uncorrelated white noise, ψα,i(x,t)ψβ,j(x,t)=−1δα,βδi,jδ(x−x)δ(t−t),(5) and 1 is a large parameter that characterizes the separation between the microscopic and macroscopic scales. The condition of having an equilibrium reference state implies a relation between gand σ. We know from MFT [6,8] that the deterministic current that describes a system at equilibrium should be of the form JD α,i(x)=−1 2 γ k χα,i;γ,j(φ(x))∂k δVeq[φ] δφγ(x),(6) where χα,i;β,j(φ;x)= γ k σα,i;γ,k(φ(x))σβ,j;γ,k(φ(x))(7) and Veq is the equilibrium mesoscopic potential that defines the equilibrium probability distribution: Peq[φ]≃exp [−Veq[φ]],→∞.(8) In order to get the linear form (3) from (6) we need to assume χα,i;β,k(φ)=2Lα,β (φ)δi,k, δVeq[φ] δφα(x)=−∂˜s(φ) ∂φαφ=φ(x)+cte,(9) where Lis a positive defined symmetric matrix by construction and ˜s(φ) is a function of Mvariables. Finally, we find that g=LS,(10) Sα,β (φ)=∂2˜s(φ)/∂φα∂φβ. In conclusion, the Langevin equation for DSe models is completely determined by giving the L symmetric matrix and the function ˜s(φ). From Eq. (9) we may deduce, by a simple integration, the particular form of the equilibrium potentials that give rise to this linear set of currents: Veq[φ]=− dx˜s(φ(x))−˜s(φeq) − α (φα(x)−φeq,α )∂˜s(φ) ∂φαφ=φeq ,(11) where we have made use of the known properties Veq[φeq]=0 and δV[φ]/δφα|φ=φeq =0. φeq are the equilibrium values of the fields. The deterministic current structure (3) reminds us of the macroscopic linear laws we observe in nature as Fick’s law for diffusion or the Fourier’s law for heat conduction. We know that there is a set of theories that describe how a system characterized by its mesoscopic variables fluctuates around its equilibrium state (see Einstein’s theory of fluctuations [17] and Appendix A) or how it relaxes toward the equilibrium (Onsager’s theory [5]). Both of them are contained in the so-called nonequilibrium thermodynamics [12]. Let us connect our description above with this classic point of view. Let s(φ) be the entropy per unit volume of a system in equilibrium with macroscopic observables φ=(φ1,φ 2,...,φ M). It is reasonable to think that a macroscopic system relaxing to the equilibrium state from a nearby initial state is locally at equilibrium with an entropy s(φ(x))at each macroscopic point x∈in the system. It is again assumed that s(φ) does not contain an explicit dependence on position xdue to the action of an external field such, for instance, gravitation. In this situation, the macroscopic currents associated to the conserved fields, φ, are found to have the form JD α,i(x)= β k ˜ Lα,i;β,k(φ(x))Xβ,k(φ(x)),(12) where ˜ L(the Onsager’s coefficients) is a symmetric matrix on the (α, i) index and X’s are the so-called thermodynamic forces that are defined from the local entropy s(φ): Xα,i(φ)=∂i˜yα(φ),˜yα(φ)=∂s(φ) ∂φα .(13) Observe that this classical description coincides with ours above if we identify ˜s(φ)=s(φ), the thermodynamic entropy, and ˜ Lα,i;β,k=Lα,βδi,k. Moreover, the expression (11) can be derived from equilibrium statistical mechanics (see Appendix A). Our model includes the classical description of how the macroscopic variables of systems perturbed from their equilibrium state evolve effectively toward it by assuming that local equilibrium is fulfilled. A final property is asked for JD: The equilibrium state is stable under small perturbations at the deterministic level. That is, the deterministic evolution starting from any initial set of fields near the equilibrium should relax toward it. The deterministic evolution equation is ∂tφD α(x,t)=−∇ β gαβ (φD)∇φD β.(14) 024107-3
P. L. GARRIDO PHYSICAL REVIEW E 106, 024107 (2022) Let us assume that φα(x,t)=φeq,α +θ(x,t) with θsmall and we expand the deterministic equation up to first order in θ: ∂tθα(x,t)=− β gαβ (φeq )∇θβ(x,t)+O(θ2).(15) We see that the evolution is characterized by the gmatrix evaluated at equilibrium. We rewrite this evolution equation for the Fourier transform of θ,ˆ θ: ∂tˆ θα(k,t)=k2 β gαβ (φeq )ˆ θβ(k,t) (16) and then expand ˆ θin the eigenvector basis of g: ˆ θα(k,t)= n an(k,t)vn,α,gvn=λnvn,(17) and the resulting evolution equation for a’s is given by ∂tan(k,t)=k2λnan(k,t) (18) whose solution is an(k,t)=an(k,0)exp[λnk2t].(19) The evolution of angoes to zero and then the equilibrium state is stable if and only if all the eigenvalues of ghave their real part negative: Re(λn)<0∀n.(20) Throughout the paper, we will consider only gmatrices that can be diagonalized and have positive the real part of all eigenvalues, and these cases are the most common ones. Nevertheless, further research is going to be needed to study other possibilities. These equations describe the dynamics of a DSe relaxing to the equilibrium state whenever the boundary conditions are compatible with such state, for instance, φ(x)=φeq ∀x∈ ∂. If we change such boundaries, the stationary distribution is no longer the equilibrium one defined by Veq[φ]. Moreover, the local equilibrium property is lost, and the stationary state’s quasipotential has a nonlocal structure that implies long-range correlations (see, for instance, Refs. [9–11]). This model permits us to answer some interesting and refined questions: What happens near the equilibrium? How is local equilibrium lost? What are the fluctuations of the observables? Once we have constructed the model, let us study its behavior in a generic nonequilibrium stationary state. First, the deterministic stationary solution would depend on x,φ∗(x), and it is solution of ∇JD α(x)=0⇒∇ β gαβ (φ∗(x))∇φ∗ β(x)=0.(21) The boundary conditions affects dramatically the system behavior. Typically they are assumed to be of Dirichlet type, φ(x)=φ0(x)∀x∈∂. As a helpful example, we discuss in Appendix Bthe conditions on the system’s boundaries and/or in the system’s dynamics when we require to have constant currents jα=βgαβ (φ∗) ∇φ∗ β. Moreover, there could also be global conservation laws of a field, dxφ∗ α(x)=cte, so that we also discuss their effect in the correlations in Appendix C. The fluctuating properties of such nonequilibrium stationary states are studied by using the Fokker-Planck equation associated to the above Langevin equation: ∂tP[φ;t]= α i dx∂i δ δφα(x)−JD α,i(x)P[φ;t] +1 β∂i δ δφβ(x)[Lα,β (φ)P[φ;t]].(22) The stationary distribution when →∞is of the form Pst [φ]≃exp [−V0[φ]],(23) where V0[φ] is called the quasipotential that it is solution of the Hamilton-Jacobi equation: 0= α i dx∂i δV0[φ] δφα(x) ×JD α,i(x)+ β Lα,β (φ)∂i δV0[φ] δφβ(x).(24) This simple derivation of the Hamilton-Jacobi equation hides a set of important quasipotential properties that we do not address here. We ask the reader to look at Refs. [6,7]fora complete description of them. The quasipotential contains all the relevant behavior about the system’s stationary state, but it is not easy to get explicit solutions from the Hamilton-Jacobi equation for generic cases [11]. However, let us show that from the Hamilton-Jacobi equation, we can derive a set of closed equations for the equal-time correlation functions of the stationary state. We know that these capture the essential features of the system’s spatial structure and are closely related to the quasipotential shape around the stationary state. III. EQUAL-TIME CORRELATION FUNCTIONS The correlations for our M-field model are defined as ˜ Cα1α2...αn(x1,x2...xn) ≡φα1(x1)−φα1(x1)st ...φαn(xn)−φαn(xn)st st , (25) where ·st =Dφ·Pst [φ]. In the weak noise limit (large values of ), we can use the quasipotentical V0to compute the correlations. It is a matter of algebra to show that Cα1α2...αn(x1,x2...xn)≡lim →∞n−1˜ Cα1α2...αn(x1,x2...xn) =− δF[φ∗[b],b] δbα1(x1)...δbαn(xn)b=0 ,(26) where F[φ,b]=V0[φ]− M α=1 dxbα(x)φα(x) (27) and φ∗[b]issolutionof δF[φ,b] δφα(x)φ=φ∗[b]=0⇔δV0[φ] δφα(x)φ=φ∗[b]=bα(x).(28) 024107-4
CORRELATIONS IN NONEQUILIBRIUM DIFFUSIVE SYSTEMS PHYSICAL REVIEW E 106, 024107 (2022) φ∗[0] =φ∗is the stationary solution of the Langevin equation without noise given by Eq. (21). We construct a set of closed equations for the correlations by using the HamiltonJacobi equation (24) with φ→φ∗[b] and then expanding the equation in powers of b’s (see for instance Ref. [8]fora detailed computation for the M=1 case). We get, at order b2, the general equations for the two-body correlations: β dy[Kαβ (x,y)Cβγ (y,z)+Kγβ(z,y)Cβα(y,x)] =2 i ∂xi∂zi[Lα,γ (φ∗(x))δ(x−z)],(29) where Kαβ (x,y)=δ∇JD α(x) δφβ(y)φ=φ∗ .(30) We can now substitute the JDcorresponding to the DSe (3) and we get β∇x[ bαβ (x)Cβγ (x,z)+gαβ (φ∗(x))∇xCβγ (x,z)] + β∇z[ bγβ(z)Cβα(z,x)+gγβ(φ∗(z))∇zCβα(z,x)] =2∂xi∂zi[Lαγ (φ∗(x))δ(x−z)],(31) where bαβ (x)= γ ∂gαγ (φ) ∂φβφ=φ∗∇xφ∗ γ(x).(32) In the equilibrium case φ∗(x)=φeq. Therefore, bαβ (x)=0 and the equations (31) become β gαβ (φeq )∇2 xC(0) βγ (x,z)+ β gγβ(φeq )∇2 zC(0) βα (z,x) =−2Lαγ (φeq )∇2 xδ(x−z) (33) whose solution for boundary conditions such that φ(x)= φeq ∀x∈∂ is C(0) αβ (x,y)=−(S−1)αβ (φeq )δ(x−y).(34) This result could be obtained directly from the equilibrium quasipotential Veq[φ]: C(0) αβ (x,y)=δ2Veq[φ] δφα(x)δφβ(y)−1 (φeq)δ(x−y).(35) In this paper, we are going to consider only the case of equilibrium fluctuating boundary conditions: Cαβ (x,z)= C(0) αβ (x,z)∀xor z∈∂. At this point, it is convenient to decompose the correlations in two terms, one that represents the local-equilibrium contributions (equilibrium correlations evaluated at each macroscopic point with the corresponding field values of the stationary state) and the rest that contains the strong nonequilibrium behavior: Cαβ (x,y)=−(S−1)αβ (φ∗(x))δ(x−y)+¯ Cαβ (x,y).(36) After substituting Eq. (36) into Eq. (33), we obtain the central equation for the two-body correlations: β∇x[ bαβ (x)¯ Cβγ (x,z)+gαβ[φ∗(x)]∇x¯ Cβγ (x,z)] + β∇z[ bγβ(z)¯ Cβα(z,x)+gγβ[φ∗(z)]∇z¯ Cβα(z,x)] =[∇x Aαγ (x)]δ(x−z)+[ Aαγ (x)− Aγα(x)]∇xδ(x−z), (37) where Aαγ (x)= β [ bαβ (x)(S−1)βγ +gαβ[φ∗(x)]∇x(S−1)βγ ] (38) with boundary conditions: ¯ Cαβ (x,z)=0∀xor z∈∂.This equation has the symmetry (α, x)↔(γ,z) and also that the ∇δterm does not exist in the one-field case. We see that these coupled equations for the correlation’s excess are highly nonlinear because it depends on the nonequilibrium stationary state φ∗(x), on the equilibrium reference state represented by the entropy Hessian S, and on the diffusive model g.We are interested in studying the role of the local equilibrium at the level of correlations. Therefore, we will expand these correlations near the equilibrium state to get some generic results on their properties. We should mention that our system is typically open because of the boundary conditions. However, we could think of models where some fields have global conservation constraints, for instance, in a system of particles enclosed in a container where only energy is exchanged at the boundaries. The field corresponding to the density is precisely conserved at any time, but, in contrast, the field associated with the energy is not strictly conserved. In Appendix C, we study the effect in the correlations of the existence of global conservation in some fields. We show there that the correlations CSC αβ (x,y) of a set of global conserved fields, ˜ M, can be expressed as combinations of the correlations corresponding to the nonconserved case, COB αβ (x,y): CSC αβ (x,y)=COB αβ (x,y)− ¯α¯ β∈˜ M dz1 dz2 ×COB α¯α(x,z1)(A−1)¯α¯ βCOB ¯ ββ (z2,y),(39) where Aαβ = dx dyCOB αβ (x,y)α,β ∈˜ M.(40) Therefore, global conservation do not introduce new complexities at this level and we just focus in cases where all the fields are globally nonconserved. IV. NONEQUILIBRIUM CORRELATIONS NEAR THE EQUILIBRIUM: TWO THEOREMS The DSe are driven from an equilibrium state to a nonequilibrium stationary state by changing the boundary conditions. Let us assume that the system’s stationary state is near the equilibrium. Therefore, a parameter 0 ⩽1 represents the distance of the values of its boundaries to their corresponding 024107-5
P. L. GARRIDO PHYSICAL REVIEW E 106, 024107 (2022) equilibrium ones. Then, let us assume that the deterministic stationary state, φ∗(x), can be analytically expanded: φα(x)∗=φeq,α +h(1) α(x)+2h(2) α(x)+O(3).(41) Then, from Eq. (21) we find that h(1,2) αare solutions of β gαβ (φeq)∇2 xh(1) β(x)=0,(42) β gαβ (φeq )∇2 xh(2) β(x)=− β∇xg(1) αβ (x)∇xh(1) β(x),(43) where g(1) αβ (x)= γ ∂gαγ ∂φβφeq h(1) γ(x) (44) with given boundary conditions. For instance, in a onedimensional system in a unit box [0,1], when φα(0) = φeq,α +and φα(1) =φeq,α then h(1,2) α(1) =0 and h(1) α(0) = 1, h(2) α(0) =0. When →0, the correlations tend to their equilibrium value Cαβ (x,y)→C(0) αβ (x,y) and therefore, ¯ Cαβ (x,y)→0in such limit. Thus, we can assume the existence of an analytic expansion for the correlation’s excess, ¯ C: ¯ Cαβ (x,y)=¯ C(1) αβ (x,y)+2¯ C(2) αβ (x,y)+O(3).(45) We substitute Eqs. (41) and (45)into(37) and we get a hierarchy of closed equations that for ¯ C(1) αβ (x,y) and ¯ C(2) αβ (x,y) are β gαβ∇2 x¯ C(1) βγ (x,z)+ β gγβ∇2 z¯ C(1) αβ (x,z) = a(1) αγ ∇xδ(x−z),(46) β gαβ∇2 x¯ C(2) βγ (x,z)+ β gγβ∇2 z¯ C(2) αβ (x,z) =∇ x A(2) αγ (x)δ(x−z)+ a(2) αγ (x)∇xδ(x−z) − β g(1) αβ (x)∇2 x¯ C(1) βγ (x,z)− β g(1) γβ(z)∇2 z¯ C(1) αβ (x,z) − β b(1) αβ (x)+∇xg(1) αβ (x)∇x¯ C(1) βγ (x,z) − β b(1) γβ(z)+∇zg(1) γβ(z)∇z¯ C(1) αβ (x,z),(47) where a(i) αγ (x)= A(i) αγ (x)− A(i) γα(x) (48) and A(1) αγ = ¯γ∇xh(1) ¯γ β∂gα¯γ ∂φβ−∂gαβ ∂φ¯γ(S−1)βγ ,(49) A(2) αγ (x)= σ∇xh(2) σ(x)+ ση∇xh(1) σ(x)h(1) η(x)∂ ∂φη × β (S−1)βγ η ∂Lαη ∂φβ Sησ ,(50) b(1) αγ (x)=∇ xg(1) αγ (x).(51) We have simplified the notation: gαβ ≡gαβ (φeq ) and ∂gαβ/∂φ¯γ≡∂gαβ (φ)/∂φ¯γ|φ=φeq . In general, after any operation, a functional that depend on φis considered to be evaluated at φeq. At this point we find the first general result: Therorem 1. All DSe systems with one field, M=1, have ¯ C(1) =0. That is, the excess of correlations is, at most, of order 2. That is due because A(1) =0 and the solution of Eq. (46)is a harmonic function whose maximum or minimum should be at the boundary that in our case is always zero: ¯ C(1)(x,z)=0 ∀xor z∈∂ and therefore ¯ C(1)(x,z)=0∀x,z. This property was already observed in two specific one-dimensional models, the symmetric simple exclusion process (SSEP) [9] and the Kipnis, Marchioro, and Presutti model (KMP) [10]. In these works, it is shown that ¯ C=2F(x,z) for any . In order to go forward, we need to give specific boundary conditions. Our natural choice is to place our system between two parallel plates placed at x1=0 and Lwhere the values of the fields are given: φα(0,x⊥)=φeq,α,φ α(L,x⊥)=φeq,α +φα,(52) where φαare given constants. Therefore, ¯ C(1) α,β (x,z)=0,xand/or z∈∂ ={(0,w⊥)}∪{(L,w⊥}. (53) We also assume periodic boundary conditions in the perpendicular d−1 directions: ¯ Cαβ (x1,x⊥±ηLaj)= ¯ Cαβ (x1,x⊥),∀j=2,...,d, where x=(x1,x⊥) and aiare the unit vectors on the principal directions. η>0 is a form factor. These boundary conditions have the advantage to give us a simple stationary state around the equilibrium (see Appendix B). In particular, a(1) αγ =a(1) αγ ˆıinEq.(46), with a(1) αγ constant. We apply to the ¯ C(1) αβ (x,z) functions the sinus Fourier transform to the x1,z1coordinates because they incorporate the boundary conditions and a normal Fourier transform to the perpendicular coordinates x⊥,z⊥in Eq. (46). (See Ref. [18] for general background on Fourier transforms and their applications. We also include in Appendix Dsome properties of the Fourier transform and Fourier series that are used in the paper.) Then ¯ C(1) αβ (x,z)=1 (L)d−1 n⊥∈Zd−1 ei2π Ln⊥·(x⊥−z⊥)∞ n=1 ∞ m=1 sin nπx1 L ×sin mπz1 Lˆ C(1) αβ (n,m;n⊥),(54) where L≡ηLand ˆ Cαβ functions are solution of the equations: n2+4 η2n2 ⊥ β gαβ ˆ C(1) βγ (n,m;n⊥) +m2+4 η2n2 ⊥ β gγβ ˆ C(1) βα (n,m;n⊥) =(n,m)a(1) αγ ,(55) 024107-6
CORRELATIONS IN NONEQUILIBRIUM DIFFUSIVE SYSTEMS PHYSICAL REVIEW E 106, 024107 (2022) FIG. 1. The basic structure function ˆ S(n,m;k,θ). It is shown only the (n,m) values where ˆ Sis nonzero. The sublattices with equal parity, i.e., both nand meven or odd, have ˆ S=0. where (n,m)=−4 π2[1−(−1)n+m]nm n2−m2(n= m) =0(n=m).(56) We observe that for a given set of values (n,m,n⊥), we have an ensemble of equations with the unknowns linearly related. We can express them in matrix notation: d(n)gˆ C(1) +d(m)ˆ C(1)gT=(n,m)a,(57) where d(n)=n2+4n2 ⊥/η2and we only show the arguments that change to simplify the notation. To solve these equations, we use the assumption that the matrix gcan be diagonalized or, in other words, there is an eigenvector basis that spans the M-dimensional space. We know that under these conditions gTw(s)=λ(s)w(s),Re(λ(s))<0∀s,(58) we multiply by w(s)Tthe left of equation (57) and by w(s) its right and we can isolate the ˆ C(1) matrix components: w(s)Tˆ C(1)w(s)=¯ Gss=(n,m)w(s)Taw(s) d(n)λ(s)+d(m)λ(s)(59) or, in Cartesian coordinates ˆ C(1) =(P−1)T¯ GP−1,(60) where Pis the matrix where its scolumn is the components of w(s): Pe(s)=w(s) with e(s)i=δs,ibeing the canonical orthonormal basis. Finally, we can write ˆ C(1) in components: ˆ C(1) αβ (n,m;n⊥)= σσ Gσσ;αβ ˆ F(n,m;n⊥;σ,σ),(61) where Gσσ;αβ =(P−1)σα(P−1)σβ s s a(1) ssPsσPsσ(62) and ˆ F(n,m;n⊥;σ,σ) =(n,m) (n2+4n2 ⊥/η2)λ(σ)+(m2+4n2 ⊥/η2)λ(σ),(63) where a(1) ss=a(1) ssˆı is given by Eq. (48). Observe that the property a(1) αγ =−a(1) γα implies Gσσ;αβ =−Gσσ;βα. Please, observe that ˆ C(1) αβ (n,m;n⊥) is a linear combination of the basic structure function ˆ F(n,m;n⊥;σ,σ). We show in Fig. 1the behavior of a related function that only depends on the relation between eigenvalues: ˆ S(n,m;k,θ σσ)≡λ(σ)ˆ F(n,m;n⊥;σ,σ),(64) where k2=4n2 ⊥/η2and θ2 σσ=λ(σ)/λ(σ). Let us stress the fact that we only need to study values with θσσ<1 because of the relation ˆ Sm,n;k,1 θσσ=−θ2 σσˆ S(n,m;k,θ σσ).(65) Finally, the correlations in real space given by Eq. (54) can be written ¯ C(1) αβ (x,z)= σσ Gσσ;αβF(x,z;σ,σ),(66) where we call F the basic correlation function and it is written as F(x,z;σ,σ)=˜ F(x,z;σ,σ)−˜ F(z,x;σ,σ) ≡λ(σ)−1S(x,z;θσσ) (67) 024107-7
P. L. GARRIDO PHYSICAL REVIEW E 106, 024107 (2022) and ˜ F(x,z;σ,σ)=−8 π2 1 (L)d−1 n⊥∈Zd−1 ei2π Ln⊥·(x⊥−z⊥)∞ n=1 sin (2n−1)πx1 L∞ m=1 sin 2mπz1 L ×(2n−1)2m (2n−1)2−(2m)2 1 [(2n−1)2+4n2 ⊥/η2]λ(σ)+[(2m)2+4n2 ⊥/η2]λ(σ).(68) The field’s spatial average fluctuations are interesting observables associated with the two-body correlations. For instance, at equilibrium, these magnitudes are related to other thermodynamic properties of the system as the Einstein relation between the system’s overall energy fluctuations and its specific heat. For a system composed of Mfields, we can define the fluctuations between the fields αand βat the stationary state as αβ =(eα−e∗ α)(eβ−e∗ β)ss,(69) where eαis the spatially averaged field φα: eα=1 || dxφα(x) (70) and e∗ αis its average value at the stationary state. Fluctuations can be written as the sum of correlations: αβ =1 ||2 dx dzCαβ (x,z).(71) This expression for the DSe is written as the sum of two contributions: αβ =leq αβ +neq αβ ,(72) where the local equilibrium contribution is leq αβ =1 ||2 dx (−S−1)αβ (φ∗(x))(73) and the remaining, the nonequilibrium part, is neq αβ =1 ||2 dx dz ¯ Cαβ (x,z).(74) When we -expand neq αβ through the correlation expansion, we get our second general result: Theorem 2. neq αβ =O(2) for all DSe with parallel plates as boundary conditions. In other words, the fluctuations for DSe systems with parallel plates as boundary conditions are, near to the equilibrium, at most of the order 2. The field’s global averaged values are very well described by the local equilibrium approximation whenever the stationary state is at the linear regime (order ). The nonequilibrium corrections appear at order 2despite their correlations that deviate from local equilibrium already at order , and they are long range. The proof of this theorem is straightforward. We compute explicitly neq αβ near to the equilibrium at first order in and in the case of parallel plates (see Sec. III above). We just substitute expression (61)into(54): neq,(1) αβ =1 π2(L)d−1 ∞ n=1 ∞ m=1 1 nm[1−(−1)n] ×[1−(−1)m]ˆ C(1) αβ (n,m;0).(75) We observe that the sums over nand min Eq. (75) run over odd values due to the factors in front of ˆ C(1). Moreover, ˆ C(1) αβ includes the factor (n,m) given by Eq. (56) that is different from zero whenever nand mhave different parities and therefore the overall result is zero. V. THE BEHAVIOR OF THE BASIC CORRELATION FUNCTION We observe that the correlations at the first order in the expansion are given by Eq. (68), which is a linear combination of Ffunctions (basic correlation function). Therefore, F contains the structural part of the nonequilibrium correlations in real space, and it is interesting to get some insight into it. Let us begin the study of Fwith the one-dimensional case (d=1). We see that we can get some idea of its behavior by doing numerically the sums in Eqs. (66)–(68) for given values of the ratio θσσ=λ(σ)/λ(σ). However, the sums converge very poorly due to the sinus functions. Therefore, we had to transform it to another form with a better numerical convergence behavior. After some algebra (see details in Appendix E), we transform Eq. (67)into S(x,z;θσσ)≡λ(σ)F(x,z;σ,σ)=1 π 1 1+θ2 σσ∞ m=1 arctan Aσσ(m;x,z)−∞ m=1 arctan Aσσ(m;z,x)+π 2sgn(x−z),(76) where Aσσ(m;x,z)=4 cos π 2¯xcosh π 2θσσ(2m−1)sinh π 2θσσ¯z cosh (πθσσ(2m−1))−cosh (πθσσ¯z)+2 cos2π 2¯x,(77) with θ2 σσ=λ(σ)/λ(σ), ¯x=2x/L−1, and ¯z=2z/L−1, and sign(x)=x/|x|when x= 0 and sign(0) =0. We show in Fig. 2the behavior of the S(x,z;σ,σ) versus (¯x=2x/L−1,¯z=2z/L−1) for θσσ=0.5. We obtain the 024107-8
CORRELATIONS IN NONEQUILIBRIUM DIFFUSIVE SYSTEMS PHYSICAL REVIEW E 106, 024107 (2022) FIG. 2. The S(x,z;θ)vs(¯x=2x/L−1,¯z=2z/L−1) for θ= 0.5. The red line shows the function S(x,x;θ). The black line is a reference (¯x,¯x,0). figure by computing Snumerically with Eq. (76). There are several points to remark. First, we see how S(x,z;θ) is zero for values (x,z) located at the boundaries. Moreover, let us observe a defined discontinuity along the line ¯x=¯zwhere two antisymmetric halves meet, forming a well-defined gap. There are some apparent rounding effects near the points (¯x,¯z)= (±1,±1) but these are just due to numerical computation difficulties. S(x,x;θ) is shown separately in Fig. 2byared line located at the middle of the gap. We can analytically compute the magnitude of the gap along the line (¯x,¯x) for a given θfrom Eq. (76) and we find S=lim →0[S(x0+,x0−;θ)−S(x0−,x0+;θ)] =1 1+θ2.(78) It seems remarkable that the size of the gap is independent of x0. Similarly, we also find the limiting value at ¯x=±1: S(L,L;θ)=−S(0,0; θ)=1 π 1 1+θ2[arctan θ−1−arctan θ]. (79) All these exact results are, of course, consistent with the numerical behavior obtained in Fig. 2. In the d=1 case, we got the qualitative behavior of Sby computing numerically part of the infinite sum (76) and/or extracting some analytical results from it. There we were lucky because we could obtain a fast converging expression that made possible its overall description. However, we have been unable to find a similar expression to Eq. (76) when d>1. We circumvent this difficulty by studying analytically the limit L→∞, where we can use the Riemann summation formula to substitute the sums by integrals. As we will see below, we pay the price of describing the correlations only in part of the domain, that is of arbitrary but not infinite size. Let us write from Eq. (67): S(x,z;θ)=˜ S(x,z;θ)−1 θ2˜ Sz,x;1 θ,(80) where ˜ S(x,z;θσσ)=λ(σ)˜ F(x,z;σ,σ) (81) and ˜ Fis given by Eq. (68). We can do the summation over msimilarly as we did for the d=1 case (see Appendix E). We get ˜ S(x,z;θ)=2 π 1 1+θ2 1 (L)d−1 n⊥∈Z exp i2π Ln⊥(x⊥−z⊥)∞ n=1 sin (2n−1)πx1 L ×2n−1 (2n−1)2+4n2 ⊥/η2cos (2n−1)πz1 L−sinh α(n,n⊥;θ)π 21−2z1 L sinh α(n,n⊥;θ)π 2(82) with α(n,n⊥;θ)=[θ2(2n−1)2+4(1 +θ2)n2 ⊥/η2]1/2.(83) We can do explicitly the first sum by using some of the Fourier sums that we derived in Appendix Dand we obtain ∞ n=1 sin (2n−1)πx1 L2n−1 (2n−1)2+4n2 ⊥/η2cos (2n−1)πz1 L =π 8cosh |n⊥|π η−1(1 −1 L(x1+z1)) cosh |n⊥|π η1−2x1+z1 L−1 L(x1+z1)−1 ×cosh |n⊥|π η2x1+z1 L−3+sign(x1−z1) cosh |n⊥|π η1−2|x1−z1| L.(84) 024107-9
P. L. GARRIDO PHYSICAL REVIEW E 106, 024107 (2022) FIG. 9. D12(x,z;g)forθ=0.5 from Eq. (116). Behavior for fixed value of x1=1 (top left) and fixed value of |x2−z2|=1,0.1 (top right and bottom respectively). The dashed lines show the asymptotic behavior along the directions (see text). the energy, φ2(x,t) (average kinetic energy per unit volume). The reference equilibrium state for this model is the ideal gas. The entropy per unit volume is given by the Sakkur-Tetrode expression which for dimension d=2is s(φ)=φ1ln 2πφ2 φ2 1+1(118) where we have considered h=1, kB=1, and m=1. Therefore, the mesoscopic reference equilibrium potential (119)is Veq[φ]= dxφ1(x)2lnφ1(x) φ1,eq −ln φ2(x) φ2,eq +φ1,eqφ2(x) φ1(x)φ2,eq +φ1,eq φ1(x)−2,(119) where (φ1,eq,φ 2,eq) is the macroscopic equilibrium state. In the BGK approximation [22] (where the usual interaction hardcore kernel is approximated by the local Maxwellian minus the one-particle distribution), we obtained for dimension d=2 the diffusion equations (3): g=0−1 2ωφ2 2/φ2 1−2(1 +ω)φ2/φ1(120) where ω∈[0,1]. In Ref. [16] two parameters appear: νthat is related with the collision part in the BGK approximation and αthat controls the frequency of the randomization mechanism for the particle velocities. In order to simplify computation we have assumed ω=2πα/(ν+2πα) and α=1/2π. Observe that the eigenvalues of gT(that we need for later computations) are λ(1,2) =−(1 +ω∓√1+ω2)φ2/φ1. The matrix Lcan be obtain by using Eq. (10) once we know Sfrom (118): S=−2/φ11/φ2 1/φ2−φ1/φ2 2⇒L=gS−1 =φ22φ2 2/φ1 2φ2 2/φ12(2 +ω)φ3 2/φ2 1.(121) At this point, we have all the ingredients to write down the closed equations for the static two-body correlations once we detail the boundary conditions and we find the corresponding deterministic stationary state [φ∗ 1(x),φ∗ 2(x)]x∈. Let us assume that our system is in a strip of width unity where we impose temperatures T0,T1at x1=0,Lrespectively. That permits only a flow of energy in the xdirection. We also assume that there is not a net flow of particles through the system and the average value of the density is a given constant. The stationary state is a solution of the equations (21) 024107-16
CORRELATIONS IN NONEQUILIBRIUM DIFFUSIVE SYSTEMS PHYSICAL REVIEW E 106, 024107 (2022) with constant currents: JD 1,i[φ∗;x]=0,JD 2,i[φ∗;x]=Jδi,1,i=1,2.(122) Both conditions imply φ∗ 1(x)=φ∗ 2 T(x1),T(x1)=T0−Tx1 L, φ∗ 2(x)=¯nT ln(T0/T1)≡φ∗ 2,(123) with T=T0−T1and ¯n=L−1L 0dxφ∗ 1(x) that are our system’s parameters. Let us note here that φ∗ 1(x) is the particle density and φ∗ 2(x) is the energy density. A bidimensional ideal gas at equilibrium has the equation of state φ2,eq =φ1,eqTeq where φ2,eq is the pressure. Then, the system’s stationary state has the local equilibrium property, that is, the local pressure is constant all over the system because there is no net flow of particles. We know from Eq. (36) that correlations are decomposed in the sum of two terms: the local equilibrium contribution, CLE αβ (x,y),and the correlation’s excess, ¯ Cαβ (x,y), that is solution of Eq. (37). We easily compute CLE: CLE αβ (x,y)=−(S−1)αβ (φ∗(x))δ(x−y),(124) where S−1(φ)=− φ1φ2 φ22φ2 2/φ1.(125) To compute ¯ C, we need to do an expansion around the equilibrium. For this model, we use T/Las the parameter in the expansions we defined in Sec. IV. We get all necessary items there: A(1,2) αβ ,g(1) αβ, by expanding the stationary state (123) for T≃0 to obtain h(1) 1(x)=L¯n T1x1 L−1 2, h(2) 1(x)=L2¯n T2 11−x1 L1 2−x1 L−1 12,(126) h(1) 2(x)=L¯n 2,h(2) 2(x)=−L2¯n 12T1 . ¯ C(1) αβ was extensively studied in Sec. V. However, we need their explicit expressions to study the next expansion order: ¯ C(2) αβ . One can check that Eq. (61) implies ˆ C(1) 11 (n,m;n⊥)=− 4¯n π2T1 (1 −ω)(1−δn,m) ×[1−(−1)n+m]ˆ D11(n,m;n⊥), ˆ C(1) 12 (n,m;n⊥)=−8¯n π2(1 −ω2)(1−δn,m)(127) ×(1−(−1)n+m)ˆ D12(n,m;n⊥), ˆ C(1) 22 (n,m;n⊥)=2ωT2 1ˆ C(1) 11 (n,m;n⊥), where ˆ D11(n,m;n⊥)=nm D(n,m;n⊥), ˆ D12(n,m;n⊥)=nm D(n,m;n⊥) n2+4n2 ⊥ n2−m2,(128) and D(n,m;n⊥)=ω(n2−m2)2+(1 +ω)2(n2+4n2 ⊥)(m2+4n2 ⊥).(129) After making the sinus Fourier transform of Eq. (47), we get a set of four linear equations with unknowns ˆ C(2) αβ . The solution is ˆ C(2) 11 (n,m;n⊥)=L2 2π2T2 1ω 1 D(n,m;n⊥)(n2+m2+8n2 ⊥)2 T1 (1 +ω)(m2+4n2 ⊥)(n2+4n2 ⊥)R2(n,m;n⊥) +(m2+4n2 ⊥)[2(1 +ω)2(n2+4n2 ⊥)−ω(n2−m2)]R1(n,m;n⊥) +(n2+4n2 ⊥)[2(1 +ω)2(m2+4n2 ⊥)+ω(n2−m2)]R1(m,n;n⊥), ˆ C(2) 12 (n,m;n⊥)=L2 π2T1 n2+4n2 ⊥ D(n,m;n⊥)(n2+m2+8n2 ⊥)−1 T1 (n2−m2)R2(n,m;n⊥) +(1 +ω)(m2+4n2 ⊥)R1(m,n;n⊥)−(1 +ω)(n2+4n2 ⊥)R1(n,m;n⊥), ˆ C(2) 21 (n,m;n⊥)=−m2+4n2 ⊥ n2+4n2 ⊥ ˆ C(2) 12 (n,m;n⊥), ˆ C(2) 22 (n,m;n⊥)=L2 π2 1 D(n,m;n⊥)(n2+m2+8n2 ⊥)2 T1 (1 +ω)(m2+4n2 ⊥)(n2+4n2 ⊥)R2(n,m;n⊥) +ω(n2−m2)[(m2+4n2 ⊥)R1(m,n;n⊥)−(n2+4n2 ⊥)R1(n,m;n⊥)],(130) 024107-17
P. L. GARRIDO PHYSICAL REVIEW E 106, 024107 (2022) where R1(n,m;n⊥)=4¯n L(1 −ω)nmω(1 −δn,m)(1−(−1)n+m)n2+m2+8n2 ⊥ D(n,m;n⊥) +16 π2ω[1+(−1)n+m](m2+4n2 ⊥) m=m [1 −(−1)m+m]m2 D(n,m;n⊥)(m2−m2)2 +8 π2(1 +ω)[1+(−1)n+m](n2+4n2 ⊥) m=m [1 −(−1)m+m]m2 D(n,m;n⊥)(m2−m2)2 ×1 n2−m2[(1 +3ω)(m2−m2)−2(1 +ω)(m2+4n2 ⊥)],(131) R2(n,m;n⊥)=4¯nT1 L(1 +2ω)δn,m−4¯nT1 Lω(1 −ω2)(1 −δn,m)[1−(−1)n+m]nm n2+m2+8n2 ⊥ D(n,m;n⊥) +64¯nT1 π2Lω(1 −ω2)[1+(−1)n+m]nm n=n [1 −(−1)n+n]n2(n2+4n2 ⊥) (n2−n2)(m2−n2) ×m2+4n2 ⊥ D(n,n;n⊥)(m2−n2)+n2+4n2 ⊥ D(n,m;n⊥)(n2−n2)+16¯nT1 π2Lω(1 +3ω)(1 −ω)[1+(−1)n+m]nm × n=n [1 −(−1)n+n]n21 D(n,n;n⊥)(m2−n2)+1 D(n,m;n⊥)(n2−n2)+32¯nT1 π2Lω(1 −ω2) ×[1+(−1)n+m]nm n=n [1 −(−1)n+n]n2(n2+4n2 ⊥)1 D(n,n;n⊥)(m2−n2)2+1 D(n,m;n⊥)(n2−n2)2. (132) We are interested in using these solution for ˆ C(2) αβ to compute the fluctuations of the field’s spatial average (74) at this order: neq,(2) αβ =1 π2L ∞ n=1 ∞ m=1 1 nm[1−(−1)n][1−(−1)m]ˆ C(2) αβ (n,m;n⊥=0).(133) Observe that to compute αβ all we need is ˆ Cαβ (n,m;n⊥)forn⊥=0 and n,modd values. Therefore, some sums appearing in Eqs. (131) and (132) just disappear and others should be done on even values. That permits us to obtain explicitly expressions for all of them: ˆ C(2) 11 (n,m;0) =4L¯n π2T2 1ω nm D(n,m;0)(n2+m2)(1 +ω)(1 +2ω)n2δn,m−16 π2(1 −ω)A11(n,m), ˆ C(2) 12 (n,m;0) =64L¯n π4T1 n3m D(n,m;0)(n2+m2)A12(n,m),(134) and A11(n,m)= 2 l=0α(l) 11 (n,m)Bl(n,m)+α(l) 11 (m,n)Bl(m,n), A12(n,m)= 2 l=0α(l) 12 (n,m)Bl(n,m)−α(l) 12 (m,n)Bl(m,n)(135) with α(0) 11 (n,m)=m4n2(2(1 +3ω+5ω2+5ω3+3ω4)n2+ω(1 +2ω+3ω2)m2), α(1) 11 (n,m)=2m2(ω(1 +ω+ω2)n2m2+ω2m4+(1 −ω2)(1 +ω+ω2)n4),(136) α(2) 11 (n,m)=ω(1 −ω2)n2m2, α(0) 12 =ω(1 −ω)(n2−m2), α(1) 12 =(1 −ω)(1 +ω+ω2)n4+ω(3 +ω)m4,(137) α(2) 12 =(1 +2ω+2ω2+3ω3)n2+ω(1 +3ω)m2, 024107-18
CORRELATIONS IN NONEQUILIBRIUM DIFFUSIVE SYSTEMS PHYSICAL REVIEW E 106, 024107 (2022) FIG. 10. Second-order fluctuations for the hydrodynamic twofield model computed using Eq. (133). Black dots: ¯ neq,(2) 11 (ω)= neq,(2) 11 (ω)T2 1/¯n. Red dots: ¯ neq,(2) 12 (ω)=neq,(2) 12 (ω)T1/¯n. The black dashed line is the analytic asymptotic behavior for small values of ω.¯ neq,(2) 11 (1) =0.061685, ¯ neq,(2) 12 (1) =0,and ¯ neq,(2) 12 (0) = −0.00210306. where the sums Bl(n,m)=∞ k=0 (2k)2+2l D(n,2k;0) 1 [(2k)2−m2)2(n2−(2k)2] (138) are explicitly done and given in Appendix D. With all these ingredients, we can compute neq,(2) αβ givenby(133). The analytic expressions are long functions of ωthat we do not explicitly write here. We have plotted their behavior in Fig. 10. We see that ¯ neq,(2) 11 (ω)=neq,(2) 11 (ω)T2 1/¯nbehaves effectively like a 1/ω function. In fact, its behavior for small values of ω is ≃0.06767/ω. This singularity for ω=0 is expected. When ω=0, the microscopic model loses the velocity randomization mechanisms. Its mesoscopic description dramatically changes because the momentum is locally conserved at such a limit; therefore, there are two more conserved fields and the model with two fields breaks down. On the other hand, ¯ neq,(2) 11 (1) =0.061685. This value differs by only by 0.006 from the extrapolation to one of the asymptotic expression around ω=0: ¯ neq,(2) 11 (ω)≃0.06767/ω. That is, the analytic complex and long formula for ¯ neq,(2) 11 (ω) only accounts for a tiny correction of simple asymptotic formula. Finally, we see that ¯ neq,(2) 12 (ω)=neq,(2) 12 (ω)T1/¯nis about 103times smaller than ¯ neq,(2) 11 (ω) and negative. It has a finite limit for ω=0 and a minimum near it. This fluctuation does not reflects the change on the mesoscopic description when ω→0. VIII. CONCLUSIONS We have studied the two-body equal time correlation functions for a diffusive sytem with a reference equilibrium state (DSe) with Mfields. We have derived the partial differential equations they follow and studied explicitly their solutions perturbatively around the equilibrium. We show the correlation’s complex and rich behavior as already observed in some particular nonequilibrium situations [15,21]: Generic power laws that depend on the path we follow when doing the longdistance limit. We show that the DSe correlations have two levels. The first level is the basic correlation function that is generic and it does not depend on the specific model, just on the eigenvalues’s ratio of the matrix gdefining the proportionality of fluxes and forces. It contains the basic spatial structure of the model. The second level is the linear combination of such basic correlation function to build the correlations. That combination strongly depends on the model’s details and on the form of the stationary state. Therefore, it seems interesting to define models that focus only on the basic correlation functions in order to study the main generic properties of those nonequilibrium systems. ACKNOWLEDGMENTS This work is part of the Project of I+D+i Ref. No. PID2020-113681GB-I00, financed by MICIN/AEI/10.13039/501100011033 and FEDER “A way to make Europe.” APPENDIX A: MESOSCOPIC PROBABILITY DISTRIBUTION OF EQUILIBRIUM Any system at equilibrium is completely determined by a small set of macroscopic variables. At mesoscopic level, such observables fluctuate around its equilibrium values. In order to obtain the corresponding probability distribution we can follow two equivalent strategies: the grand canonical ensemble or the Einstein fluctuation theory [17]. In this Appendix, we apply both paths to the case of a system with an equilibrium state defined by the number of particles, N,thevolume,, and the energy, E. 1. Grand canonical ensemble In the grand canonical ensemble, the equilibrium state is determined by (T,μ,) (temperature, chemical potential, and volume respectively). We know that these set of variables are related with (neq,eeq ) (particle density and energy per particle respectively) by β=∂s(n,e) ∂eeq , βμ =βeeq −s(neq,eeq )−neq ∂s(n,e) ∂neq ,(A1) where s(n,e) is the entropy per particle and β=1/T. The probability of finding the system with a given energy per particle, e=E/N, and density, n=N/, is given by P(n,e)=−1∞ N=1 eβμN dxNRNd dp Ne−βH(xN,pN) ×δe−H(xN,pN) Nδn−N .(A2) This relation simplifies when →∞, P(n,e)≃exp [−(f(neq,eeq)−f(n,e))](A3) 024107-19
P. L. GARRIDO PHYSICAL REVIEW E 106, 024107 (2022) with f(n,e)=n[βμ −βe+s(n,e)](A4) that can be written P(n,e)≃exp [−Veq (n,e)],(A5) Veq(n,e)=ns(neq,eeq )−s(n,e)−(eeq −e)∂s(n,e) ∂eeq +neq n(n−neq)∂s(n,e) ∂neq.(A6) Observe that Veq(neq,eeq)=0, ∂Veq(n,e)/∂n|eq=0,and ∂Veq(n,e)/∂e|eq=0. 2. Einstein fluctuation theory Boltzmann proposed that the entropy of a system at an equilibrium state defined by the macroscopic variables (A,B) is related with the number of compatible microstates, ω(A,B): S(A,B)=ln w(A,B),(A7) wherewehaveassumedkB=1. When we relax the constraint that fixed A, the system spontaneously evolves to a new equilibrium state defined by only the observable B. Then w(B)= A w(A,B),(A8) where the sum run over all the possible values of the observable A. Therefore, eS(B)= A eS(A,B).(A9) We consider that Sis a extensive variable proportional to a volume and therefore the sum is dominated by the term of the sum that maximizes S(A,B), that is, S(B)=S(A∗,B),∂S(A,B) ∂AA=A∗=0.(A10) Observe that w(B)≃w(A∗,B) and the state A∗is the one with larger number of microstates compared with any other value: w(A∗,B)⩾w(A,B). That is, the number of microstates compatible with the equilibrium state Bis overwhelmingly large compared to the number of microstates associated to any state (A,B), w(A,B) when Ais a macroscopic deviation from A∗. Einstein, following Boltzmann’s line of reasoning, proposed that the probability that a system at equilibrium in a state Bto be in a fluctuating macrostate Ashould be the ratio between the number of microstates compatible with Aand the total number of microstates compatible with B: P(A|B)=ω(A,B) ω(B)≃eS(A,B)−S(A∗,B).(A11) We can apply this idea to our example. The equilibrium state of a closed system Uis defined by the macrovariables (N,V,E) (number of particles, volume, and energy respectively). Let us divide the system in two disjoint subsystems 1 and 2. Let us assume that we have a set of constraints that can fix the equilibrium state at 1: (N1,V1,E1). That fixes the equilibrium state at 2: (N2,V2,E2)=(N−N1,V−V1,E−E1). Therefore, the total entropy of Uis just the sum of the entropies of both subsystems. Therefore, the U’s entropy per particle in this constrained system is s(n1,e1;n,e,α)=αn1 ns(n1,e1)+1−αn1 ns(n2,e2), (A12) where n=N/V,e=E/N,α=V1/V,n1=N1/V1,e1= E1/N1, and n2=n−αn1 1−α,e2=ne −αn1e1 n−αn1 .(A13) We have assumed that the entropy is extensive: S(N,V,E)= Ns(n,e)forNlarge enough. We see that the equilibrium state of the constrained system is defined by five variables: (n1,e1;n,e,α). Let us release the constraints over n1and e1 while keeping fixed n,e, and α. It can be checked that n∗ 1=n and e∗ 1=eare the values that make maximum the entropy (A12). We can apply now the Einstein theory of fluctuations. The probability to observe the subsystem 1 with values (n1,e1, while Uis at equilibrium state (n,e)is P(n1,e1;n,e,α)=exp [N(s(n1,e1;n,e,α)−s(n,e))]. (A14) It is straightforward to see that in the limit α→0 we get the same result we got using the grand canonical ensemble. APPENDIX B: STRUCTURE OF THE DETERMINISTIC STATIONARY FIELDS FOR DSe SYSTEMS The deterministic stationary fields for DSe systems, φ∗,are solutions of ∇x β gαβ (φ∗)∇xφ∗ β=0.(B1) We look for the conditions to have φ∗being solutions of jα= β gαβ (φ∗)∇xφ∗ β(B2) with j’s being constant vectors and of at least C2type. Of course, all φ∗solutions of Eq. (B2) are solutions of Eq. (B1), but that is not always true in the reverse case. In physics, jare the stationary currents and they contain implicit information of the form of the boundary conditions. We focus on asking that the cross derivatives of φ∗be equal (∂2 ijφ∗ α=∂2 jiφ∗ α)to guarantee continuity and C2differentiability. In order to study this property, we invert Eq. (B2), ∇xφ∗ β= σ (g−1)βσ jσ,(B3) and by doing the cross derivatives we get the differentiability condition α σ jα,jjσ,i γ ×∂(g−1)βσ ∂φγ (g−1)γα −∂(g−1)βα ∂φγ (g−1)γσ=0.(B4) We find several cases that accomplish Eq. (B4): 024107-20
CORRELATIONS IN NONEQUILIBRIUM DIFFUSIVE SYSTEMS PHYSICAL REVIEW E 106, 024107 (2022) (a) jα,jjσ,i=0∀i= j, or equivalently jα,i=qαδi,kfor agivenkdirection. That is, all the currents should follow the same vector direction. That is the case when, for instance, the boundaries are two hyperplanes of d−1 dimensions placed one in front of the other and with homogeneous values for the fields at the boundaries. (b) Conditions on the system: (i) Lis a constant matrix, (ii) L=˜s(φ)Awhere Ais a constant matrix and ˜sis the entropy of the reference equilibrium state, and (iii) M=1. We study in this paper systems with boundary conditions as described in the (a) case. Let us take x∈[0,L] as the axis perpendicular to the boundary hyper planes. Equation (B2)is then written jα= β gαβ (φ∗)dφ∗ β(x) dx (B5) and we take the boundary conditions φ∗ α(0) =φeq and φ∗ α(L)=φeq +φα,where φαare given constants. jαare constants that are determined by the boundary conditions. We now apply a perturbative expansion around the equilibrium solution: φ∗ α(x)=φeq +h1,α (x)+2h2,α (x)+···, jα=j1,α +2j2,α +···,(B6) where is related to the distance to the equilibrium. The original boundary conditions are translated to the h’s functions: h1,α (0) =h2,α (0) =0,h1,α (L)=φα,h2,α (L)=0. (B7) The expansion of Eq. (B5) gives the set of equations j1,α = β gαβ (φeq)dh1,β (x) dx , j2,α = βgαβ (φeq )dh2,β (x) dx + γ ∂gαβ ∂φγφ=φeq h1,γ (x)dh1,β (x) dx ... (B8) that can be solved order by order. The solutions for h1,α and h2,α (x)are h1,α (x)=φα Lx,j1,α = β gαβ (φeq)φβ L,(B9) 2h2,α (x)= β (g−1)αβ2j2,β x1−x L, 2j2,α =1 L σ γ ∂gασ ∂φγφ=φeq φσφγ,(B10) where φαare of order . APPENDIX C: CORRELATIONS FOR SYSTEMS WITH SOME STRICTLY CONSERVED FIELDS Let us assume that the system’s dynamics have a stochastic evolution that locally conserves the fields as our diffusive systems defined in the main text. Boundary conditions may or not break such a conservation law. For instance, open boundary conditions introduce fluctuations on the average field that periodic boundary condition does not. Moreover, for systems with M>1 fields, some of them may be strictly conserved while others are not. For example, think of a particle system where we permit open energy exchanges with the boundaries. Still, we fix the total number of particles, or the density field’s average is constant during the system’s evolution. This difference affects the form of the correlation functions. Let POB st [φ]≃exp[−V0[φ]] be the stationary distribution when 1 with a set of open-boundary conditions imposed on the system for all the Mfields (OB stands for open boundaries). That is, it is solution of the Fokker-Planck equation (23) in a such limit. It can be checked that the restricted distribution PSC st [φ]≃e−V0[φ] α∈˜ M δ dxφα(x)−||¯ φα(C1) is also a stationary solution of (23) compatible with the boundary conditions whenever ¯ φα=1 || dxφ∗ α(x),α∈˜ M,(C2) where φ∗ α(x) is the deterministic stationary solution for the αfield of the Langevin equation (SC stands for strictly conserved). In this case, we should think that at the microscopic level there is a strong constraint on the system that forces such strict conservation laws. Observe that in the limit →∞both systems, without or with constraints on some fields, have the same macroscopic representation. In this Appendix, we look for the relations between the two-body-correlations associated with the SC and the OB systems as we have defined them. Higher order correlations depend on other quasipotential’s perturbative terms that may be different for the OB and SC cases. We know that the two-body correlation for the OB case is related with the quasipotential’s second derivatives [8]: COB αβ (x,y)=(V−1)αβ (x,y),Vαβ (x,y)=δ2V0[φ] δφα(x)δφβ(y)φ∗ . (C3) Let us compute the correlations for the SC case. We define the functional generator: Z[B]=Dφexp {−F[φ,B]} × α∈˜ M δ dxφα(x)−||¯ φα,(C4) where F[φ,B]=V0[φ]− α∈M dxBα(x)φα(x).(C5) 024107-21
P. L. GARRIDO PHYSICAL REVIEW E 106, 024107 (2022) Observe that Mis the total number of fields and ˜ M⊆Mis the set of the strictly conserved fields. Then, the correlations are just derivatives of the functional generator (C4): CSC αβ (x,y)≡lim →∞[φα(x)φβ(y)−φα(x)φβ(y)] =lim →∞ δ δBα(x) δ δBβ(y) 1 ln Z[B]B=0 ,(C6) where ·=DφPSC st [φ]. We use now the Laplace representation of the Dirac’s delta function to get Z[B]≃ α∈˜ Mc+i∞ c−i∞ dsαDφexp{−G[φ,B,s]},(C7) where G[φ,B,s]=F[φ,B]+ α∈˜ M sα dx (φα(x)−¯ φα).(C8) We can get the dominant part of the integral (C4) when → ∞by expanding Garound the value (φ0[B],s0[B]) that make it a minimum. That is, Z[B]≃exp {−G[φ0[B],B,s0[B]]},→∞,(C9) where (φ0[B],s0[B]) are solution of the equations δG[φ,B,s] δφα(x)φ=φ0[B] s=s0[B]=0⇒δV0[φ] δφα(x)φ=φ0[B] =Bα(x)−δα∈˜ Ms0,α, ∂G[φ,B,s] ∂sαφ=φ0[B] s=s0[B]=0⇒ dxφ0,α (x) =||¯ φα∈˜ M.(C10) We see that for φ0(x;B=0) =φ∗(x) and s0[B=0] =0. Therefore, we can find the solution of Eqs. (C10) by doing a perturbative expansion around B=0. We get to first order: φ0,α (x)=φ∗ α(x)+φα(x),(C11) where φα(x)= β∈M dyCOB αβ (x,y)¯ Bβ(y), ¯ Bα(x)=Bα(x)−δα∈˜ Ms(1) 0,(C12) and s(1) 0is the first-order expansion in B.s0is the solution of dxφα(x)=0,α∈˜ M.(C13) Substituting the solution of φ0[B] into Eq. (C9), we get lim →∞ 1 ln Z[B]=−V0[φ∗]+ α∈Mλ dxBα(x)φ∗ α(x) +1 2 αβ∈M dx dyCOB αβ (x,y) ׯ Bα(x)¯ Bβ(y)+O(B3).(C14) Finally, in the form of Eq. (C6) we get the correlations for the strictly conservation case: CSC αβ (x,y)=COB αβ (x,y)− γδ∈˜ M (A−1)γδ dz1COB αγ (x,z1) × dz2COB βδ (y,z2),(C15) where Aαβ = dx dyCOB α,β (x,y),α,β∈˜ M.(C16) Please observe that dxCSC αβ (x,y)=0ifαand/or β∈˜ M,(C17) as we expected. For systems at equilibrium, we know that COB αβ (x,y)= ¯ Cαβδ(x−y). Therefore, CSC αβ (x,y)=¯ Cαβδ(x−y)−1 || γ¯γ∈˜ M ¯ Cαγ ¯ Cβ¯γ(A−1)γ¯γ, (C18) where Aαβ =¯ Cαβ,α,β∈˜ M,(C19) and they coincide in the thermodynamic limit. APPENDIX D: FOURIER TRANSFORMS AND SUMS .We use in this paper the sinus Fourier transform for the x-axis coordinates where the functions, f(x), are zero in the boundaries of the interval f(0) =f(L)=0: f(x)=∞ n=1 sin nπx Lˆ f(n),x∈[0,L].(D1) To use this transform, we need the properties 2 LL 0 dx sin nπx Lsin mπx L=δm,n,(D2) 1 LL 0 dx sin nπx Lcos mπx L =1 π[1−(−1)n+m]n n2−m2(n= m) =0(n=m).(D3) 024107-22
CORRELATIONS IN NONEQUILIBRIUM DIFFUSIVE SYSTEMS PHYSICAL REVIEW E 106, 024107 (2022) The normal Fourier transform is used for the x⊥∈D≡[0,L]d−1coordinates where the functions are periodic: g(x⊥)= n ei2π Ln·x⊥ˆg(n),n∈Zd−1,(D4) and we have the useful property 1 Ld−1D dx⊥ei2π Ln·x⊥=δn,0.(D5) We needed to derive in this work some Fourier sums: ∞ m=1 2msin(2mx) (2m)2+a2=π 4 sinh aπ 2−˜x sinh π 2a,0<˜x<π, =0,˜x=0, where ˜x=mod(x,π). ∞ m=1 (2m−1)sin ((2m−1)x) (2m−1)2+a2=sign(π−¯x)π 4 cosh aπ 2−˜x cosh π 2a,0<¯x<2π =0,¯x=0, with ¯x=mod(x,2π) and sign(0) =0. Taking a→ia,wealsoget ∞ m=1 2msin(2mx) (2m)2−a2=π 4 sin aπ 2−˜x sin π 2a,0<˜x<π, =0,˜x=0. In particular, if a=2n−1, n∈Z ∞ m=1 2msin(2mx) (2m)2−(2n−1)2=π 4cos ((2n−1)˜x),0<˜x<π, =0,˜x=0. ∞ m=1 (2m−1)sin((2m−1)x) (2m−1)2−a2=sign(π−¯x)π 4 cos aπ 2−˜x cos π 2a,0<¯x<2π, =0,¯x=0. In particular, if a=2n,n∈Z, ∞ m=1 (2m−1)sin ((2m−1)x) (2m−1)2−(2n)2=sign(π−¯x)π 4cos(2n˜x),0<¯x<2π, =0,¯x=0. In order to show these relations, we use some known result from Ref. [20]; for instance, in the first case we use Eq. (1.445.1): I(x,a)=∞ m=1 msin(mx) m2+a2=π 2 sinh (a(π−x)) sinh (πa).(D6) We separate the sum into even and odd terms: I(x,a)=Ie(x,a)+Io(x,a). But Ie(x,a)=I(2x,a/2)/2 and then we get the desired result: Io(x,a)=I(x,a)−I(2x,a/2)/2. Another relation that we use in the text is ∞ n=1 sin ((2n−1)x) (2n−1)2+a2=1 asin(x)∞ 0 dβe−βsin(βa)1+e−2β 1−2e−2βcos(2x)+e−4β. Finally, in Sec. VII, we need to solve sums of the form Bl(n,m)=∞ k=0 (2k)2+2l D(n,2k;0) 1 [(2k)2−m2)2(n2−(2k)2],(D7) 024107-23
P. L. GARRIDO PHYSICAL REVIEW E 106, 024107 (2022) where D(n,m;n⊥)isgivenbyEq.(129) and nand mare odd integers. These sums are done by breaking apart the denominators and then we use some of the above relations. After some trivial algebra, we get m= n: B0(n,m)=π 16ωπn2a0(ω) (n2−m2)(n2a0(ω)+m2)2(n2a1(ω)+m2) +1 [n2a1(ω)+m2]24√a1(ω) coth 1 2πn√a1(ω) n3[a1(ω)+1][a1(ω)−a0(ω)] −πm2[n2a1(ω)+m2] (m2−n2)[n2a0(ω)+m2]2 +4√a0(ω) coth 1 2πn√a0(ω) n3[a0(ω)+1][a0(ω)−a1(ω)][n2a0(ω)+m2]2,(D8) B1(n,m)=1 16ω−8 [a0(ω)+1][a1(ω)+1](n3−m2n)2+π2m2−8 (n2−m2)[n2a0(ω)+m2][n2a1(ω)+m2] −4a0(ω)πn√a0(ω) coth 1 2πn√a0(ω)−2 n2[a0(ω)+1][a0(ω)−a1(ω)][n2a0(ω)+m2]2+4a1(ω)πn√a1(ω) coth 1 2πn√a1(ω)−2 n2[a1(ω)+1][a0(ω)−a1(ω)][n2a1(ω)+m2]2 +8(n4{a0(ω)[m2−a1(ω)(m2−2n2])+m2a1(ω)}+m6) (m2−n2)2[n2a1(ω)+m2]2[n2a1(ω)+m2]2,(D9) B2(n,m)=1 16ω(π2m2−8)m2 (n2−m2)[n2a0(ω)+m2][n2a1(ω)+m2]−8 [a0(ω)+1][a1(ω)+1](m2−n2)2 +4a0(ω)2πn√a0(ω) coth 1 2πn√a0(ω)−2 [a0(ω)+1][a0(ω)−a1(ω)][n2a0(ω)+m2]2+4a1(ω)2πn√a1(ω) coth 1 2πn√a1(ω)−2 [a1(ω)+1][a1(ω)−a0(ω)][n2a1(ω)+m2]2 −8m2n2{a0(ω)[a1(ω)(2m2n2−3n4)+m4−2m2n2]+a1(ω)(m4−2m2n2)−m4} (m2−n2)2[n2a0(ω)+m2]2[n2a1(ω)+m2]2.(D10) n=m: B0(n,n)=1 64n7w[a0(ω)+1]3[a1(ω)+1]3[a0(ω)−a1(ω)] π(a1(ω)+1)πn(a0(ω)+1){a0(ω)[a1(ω)−3] −3a1(ω)−7}[a1(ω)−a0(ω)] +16a0(ω)[a1(ω)+1]2coth 1 2πna0(ω)−16π[a0(ω)+1]3a1(ω) coth 1 2πna1(ω),(D11) B1(n,n)=1 64n5w[a0(ω)+1]3[a1(ω)+1]3[a0(ω)−a1(ω)]16π[a0(ω)+1]3a1(ω)3/2coth 1 2πna1(ω) +π[a1(ω)+1]πn[a0(ω)+1][a1(ω)−a0(ω)][5a0(ω)a1(ω)+a0(ω)+a1(ω)−3] −16a0(ω)3/2[a1(ω)+1]2coth 1 2πna0(ω),(D12) B2(n,n)=1 64n3w[a0(ω)+1]3[a1(ω)+1]3[a0(ω)−a1(ω)] π[a1(ω)+1]πn[a0(ω)+1][a1(ω)−a0(ω)]{a0(ω)[9a1(ω)+5] +5a1(ω)+1} +16a0(ω)5/2[a1(ω)+1]2coth 1 2πna0(ω)−16π[a0(ω)+1]3a1(ω)5/2coth 1 2πna1(ω),(D13) where a0,1(ω)=1 ω(1 +ω+ω2±(1 +ω)1+ω2).(D14) Please note that these expressions only apply for nand mbeing odd integers. 024107-24
CORRELATIONS IN NONEQUILIBRIUM DIFFUSIVE SYSTEMS PHYSICAL REVIEW E 106, 024107 (2022) APPENDIX E: COMPUTATION OF THE BASIC CORRELATION FUNCTION F(x,z;σ, σ)FORd=1 The basic correlation function is defined by Eqs. (67) and (68). For dimension one, they reduce to F(x,z;σ,σ)=˜ F(x,z;σ,σ)−˜ F(z,x;σ,σ)(E1) with ˜ F(x,z;σ,σ)=− 8 π2λ(σ) ∞ n=1 sin π L(2n−1)x∞ m=1 sin 2πm Lz(2n−1)2m (2n−1)2−(2m)2 1 θ2 σσ(2n−1)2+(2m)2.(E2) Now we separate the fractions: 1 (2n−1)2−(2m)2 1 θ2 σσ(2n−1)2+(2m)2) =1 1+θ2 σσ 1 (2n−1)21 (2n−1)2−(2m)2+1 θ2 σσ(2n−1)2+(2m)2(E3) and we get ˜ F(x,z;σ,σ)=− 8 π2λ(σ) 1 1+θ2 σσ ∞ n=1 sin π L(2n−1)x 2n−1 ∞ m=1 sin 2πm Lz2m ×1 (2n−1)2−(2m)2+1 θ2 σσ(2n−1)2+(2m)2.(E4) At this point, we can use the formulas in Appendix Dto do explicitly the sum over m’s. We find ˜ F(x,z;σ,σ)=− 2 πλ(σ) 1 1+θ2 σσ ∞ n=1 sin π L(2n−1)x 2n−1−cos π L(2n−1)z+sinh π 2θσσ(2n−1)1−2z L sinh π 2θσσ(2n−1).(E5) The first sum can be done by converting the sinus cosinus product into a sum of sinus. Then, we use the Gradsteyn’s formula GR.1.442.1 [20] to get ∞ n=1 sin π L(2n−1)x 2n−1cos π L(2n−1)z=π 8[sgn(x−z)+sgn(L−(x+z))].(E6) The second sum in Eq. (E5) needs more work to get a simple version. First, we use Gradsteyn’s GR.3.743.1 that converts an hyperbolic sinus ratio into an integral: sinh(aβ) sinh(bβ)=2β π∞ 0 dy sin(ay) sin(by) 1 y2+β2;(E7) in our case we choose b=1, a=1−2z/L, and β=θσσ(2n−1)π/2. Therefore, we can write I=∞ n=1 sin π L(2n−1)x 2n−1 sinh π 2θσσ(2n−1)1−2z L sinh π 2θσσ(2n−1) =θσσ∞ 0 dy sin(y¯z) sin y ∞ n=1 sin π L(2n−1)x y2+π 2θσσ(2n−1)2.(E8) We can convert the last sum into another integral (see Appendix D) and we get I=−2 πcos π 2¯x∞ 0 dβe−β1+e−2β 1+2e−2βcos(π¯x)+e−4β∞ 0 dy sin(y¯z)sin(2 πθσσβy) ysin(y),(E9) where ¯x=2x/L−1 and ¯z=2z/L−1. We substitute the last integral with the relation that we derive in Appendix F: ∞ 0 dysin(ay)sin(by) ysin(y)=π 2sign(ab)∞ n=0 χ[2n+1−|a|<|b|<2n+1+|a|],|a|<1,(E10) where χ[condition] =1 whenever the condition holds and 0 otherwise. Therefore, we get I=−cos π 2¯x∞ n=0πθσσ(2n+1+¯z)/2 πθσσ(2n+1−¯z)/2 dβe−β1+e−2β 1+2e−2βcos(π¯x)+e−4β.(E11) 024107-25