Relaxation and fluctuations of a mass- and dipole-conserving stochastic lattice gas

Han et al. [Phys. Rev. Lett. \textbf{132}, 137102 (2024)] have recently introduced a classical stochastic lattice gas model which, in addition to particle conservation, also conserves the particles' dipole moment. Because of its intrinsic nonlinearity this model exhibits unusual macroscopic scaling behaviors, different from those of lattice gases that conserve only the number of particles. Here we investigate some basic relaxation and fluctuation properties of this model at large scales and at long times. These properties crucially depend on whether the total number of particles is infinite or finite. We find similarity solutions, describing relaxation of the dipole-conserving gas (DCG) in several standard settings. A major part of our effort is an extension to this model of the macroscopic fluctuation theory (MFT), previously developed for lattice gases where only the number of particles is conserved. We apply the MFT to the calculation of the variance of nonequilibrium fluctuations of the excess number of particles on the positive semi-axis when starting from an (either deterministic, or random) constant density at $t=0$. Using the MFT, we also identify the equilibrium Boltzmann-Gibbs distribution for the DCG. Finally, based on these results, we determine the probability distribution of, and the most probable density history leading to, a large deviation in the form of a macroscopic void of a given size in an initially uniform DCG at equilibrium.


I. INTRODUCTION
Lattice gases, composed of particles undergoing symmetric random hopping, exhibit diffusive behavior at large scales and long times [1][2][3][4].As these models conserve the number of particles locally, the gas density ρ(x, t) at large scales and at long times obeys the continuity equation with the diffusion current For a limited class of models -generically for the so called gradient lattice gases [1,3,5] -the diffusion coefficient D(ρ) can be calculated from the microscopic model exactly.In simple cases, which include noninteracting random walkers, D(ρ) is independent of the density, and Eqs.(1) and (2) reduce to the simple (linear) diffusion equation where ∆ is the Laplacian operator.How is the simple diffusion model modified if the lattice gas is constrained to satisfy additional local conservation laws (besides the conservation of the number of particles)?Systems with additional conservation laws have attracted much interest in recent years.Many studies concerned quantum systems with infinitely many conservation laws, specifically integrable systems in one spatial dimension.The corresponding generalized hydrodynamics is an active research area, see Refs.[6,7] for reviews.* meerson@mail.huji.ac.il The generalized hydrodynamics was also used for the description of the hydrodynamic behavior in classical integrable systems in one spatial dimension [8][9][10].Quantum higher-moment conserving models attracted much attention in the context of fractons [11][12][13][14][15][16][17]; see Refs.[18][19][20][21][22] for a sample of articles about the emerging 'fracton hydrodynamics'.Lattice gases with additional conservation laws provide a convenient and flexible theoretical platform for a systematic derivation (rather than a postulation) of a 'fracton hydrodynamics'.
Han et al. [23] have recently introduced an interesting stochastic lattice gas model which manifestly conserves, in addition to the number of particles, their dipole moment.The microscopic model of Han et al. involves a continuous-time lattice gas, where a randomly chosen pair of neighboring particles randomly hop in opposite directions in pairs, so that their center of mass is conserved.Han et al. obtained the continuum limit of this model via a standard derivative expansion of the deterministic rate equation for the particle densities on each site.In one spatial dimension, the resulting large-scale deterministic description is given by a nonlinear partial differential equation (PDE) of fourth order [23]: As one can see, here the simple first-order continuity equation (1) gives way to the second-order equation [12] with the current Equation ( 4) can be generalized to an arbitrary spatial dimension [23]: arXiv:2406.03462v2 [cond-mat.stat-mech]13 Aug 2024 where ∆ ≡ ∇ 2 is the Laplace's operator.The transport coefficient D comes from the microscopic model, and it has the units of length d+4 /time, where d is the dimension of space [24].
Because of an interplay of the nonlinearity and the fourth-order spatial derivative, Eqs. ( 4) and ( 7) exhibit unusual macroscopic scaling behaviors [23].Here we will further explore these scaling behaviors by deriving some similarity solutions to these equations, which describe relaxation of the dipole-conserving gas (DCG) in several standard settings.As we will see, the scaling properties of the relaxation dynamics crucially depend on whether the total number of particles is infinite or finite.
The main focus of this work, however, is on large-scale fluctuations in the DCG.A study of fluctuations obviously requires going beyond the deterministic limit, described by Eqs. ( 4) and (7).Han et al. [23] have already made this important step by deriving, from a microscopic lattice gas model, a Langevin equation for this system, see Eq. ( 19) below.In addition to the terms present in the deterministic equation ( 4), this stochastic PDE also includes a noise term.A similar-in-spirit Langevin description of the mass-only conserving lattice gases is known by the name of 'fluctuational hydrodynamics' [1][2][3][4].Starting from the Langevin equation ( 19), here we develop a macroscopic fluctuation theory (MFT), which is suitable for studying large deviations of different fluctuating quantities in the DCG.In the mass-only conserving lattice gases the corresponding MFT was developed by Jona-Lasinio et al., see Ref. [25] for a review, and it has been employed and further developed in numerous subsequent works.
Here we use the MFT to establish the form of the Boltzmann-Gibbs distribution for the DCG at equilibrium.We also apply the MFT to the calculation of the variance of nonequilibrium fluctuations of the excess number of particles on the positive semi-axis when starting from a (either deterministic, or random) constant density at t = 0. Finally, we determine the probability distribution of, and the most probable density history leading to, a large deviation in the form of void of a given size in an initially uniform DCG at equilibrium.
Here is a plan of the remainder of this paper.In Sec.II we present some similarity solutions of the deterministic Eqs. ( 4) and (7), which involve infinite and finite mass, and discuss their properties.Sections III and IV deal with fluctuations in the DCG.Starting from fluctuational hydrodynamics, as described by the Langevin equation (19), we introduce in Sec.III the problem of full statistics of the excess number of particles on the positive semi-axis.Using this setting as an example, we formulate the MFT of large deviations in the DCG and calculate the variance of the excess number of particles.Section IV is devoted to the MFT at equilibrium.Here we introduce the free energy density of the DCG at equilibrium and determine the probability distribution of, and the most probable density history leading to, the formation of a void in an initially uniform gas.Section V presents a brief summary and discussion of our main results.Some technical details of the derivation of the MFT equations and boundary conditions are relegated to Appendix A. In Appendix B we present an independent calculation of the variance of the particle excess directly from the Langevin equation (19).

II. DETERMINISTIC RELAXATION A. Infinite-mass scaling
To start with, let us study expansion of the DCG into vacuum.Suppose that the initial gas density has the form of a step-function: The relaxation of this system is described by the following similarity soluton of Eq. ( 4): In this case the dynamical exponent 4 is the same as in the linear fourth-order equation x u originally studied in the context of surface diffusion [26].
The dimensionless scaling function R(ξ) obeys an ordinary differential equation (ODE), The boundary conditions are can be obtained by solving Eqs.(10) with these boundary conditions numerically.Alternatively, we can solve numerically the full time-dependent PDE (4) after bringing it to a dimensionless form by rescaling ρ 0 x → x, ρ 5 0 Dt → t and ρ/ρ 0 → ρ. Figure 1 gives an example of such a time-dependent solution for the rescaled initial condition ρ(x, t = 0) = 1 − tanh(15 x).The top panel shows this initial condition and the resulting density profiles at rescaled times t = 5, 10 and 15.The bottom panel shows the same three density profiles, but plotted against the similarity coordinate ξ.As one can see, the profiles collapse into a single curve, which describes the scaling function R(ξ).Salient features of this similarity solution are its oscillatory decay at x → −∞ and its semi-compact support: the solution is defined for −∞ < ξ < ξ * ≃ 2.5.The asymptotic of the solution near the edge is R(ξ) ≃ (ξ * /48)(ξ * − ξ) 3 , so that the first and second derivatives of R vanish at ξ = ξ * alongside with R.

B. Finite-mass scaling
The long-time evolution of a system with a finite number N of particles is described by the similarity solution which exhibits a different dynamical exponent 5 [23].The dimensionless scaling function R(ξ) obeys the normalization condition following from the conservation of the total number of particles.The form (11) is exact if the system starts from the initial condition Otherwise, Eq. ( 11) describes a long-time asymptotic of the solution [27].For the scaling function R(ξ) we obtain an ODE 1 5 Integrating once, we obtain while the integration constant must be zero.Since R(ξ) is an even function, we can solve Eq. ( 15) on the half-line ξ > 0 with the boundary conditions The solution must be nonnegative, and the a priori unknown constant a is to be determined from the normalization condition (12).The nonnegativity of the solution and the boundary condition R(ξ → ∞) = 0 demand that the solution have a compact support.At the edges of support both R(ξ) and R ′ (ξ) must vanish, thus providing continuity of the flux, see Eq. ( 15).As one can check, the ODE (15) remains invariant under rescaling ξ → C −1/4 ξ and R → C −1 R, where C > 0. Therefore, once we have found the solution R 1 (ξ) of the problem ( 15) and ( 16) for a = 1, we can find the solution R a (ξ) for arbitrary a by the rescaling transformation Using Eqs. ( 12) and ( 17), we obtain so what remains is to find R 1 (ξ).This can be achieved numerically by the shooting method.We set a = 1 and trade the boundary condition at infinity R 1 (∞) = 0 for the condition R ′′ (0) = γ, where γ < 0 serves as the shooting parameter.Having found R 1 (ξ) and employing Eqs. ( 17) and ( 18), we obtained the numerical solution shown in Fig. 2.Here a = 0.356 . . ., while the edges of support are at |ξ| = ξ * = 3.140 . . . .(The proximity of the latter number to π raises curiosity but most likely coincidental.)The asymptotic of R(ξ) near the edges is R(ξ) ≃ (ξ * /60)(ξ * − |ξ|) 3 , so that the first and second derivatives of R vanish alongside with R at |ξ| = ξ * similarly to the step-like solution of the previous subsection.

III. MACROSCOPIC FLUCTUATION THEORY. FLUCTUATIONS OF EXCESS NUMBER OF PARTICLES
Fluctuation hydrodynamics of the DCG is described by the Langevin equation which has been recently derived in Ref. [23]: where η(x, t) is a white Gaussian noise, ⟨η(x 1 , t 1 )η(x 2 , t 2 )⟩ = δ(x 1 − x 2 )δ(t 1 − t 2 ), and we confine ourselves to one spatial dimension.
As we will show in Sec.IV, the equilibrium state of the DCG can be described by the Boltzmann-Gibbs distribution with a well-defined free-energy density F (ρ).That is, when the gas is in equilibrium at density ρ 0 , the probability density of observing an arbitrary density profile ρ * (x) is given by − ln P[ρ * (x)] ≃ S eq , where Remarkably, the free energy density of this gas, coincides with that of the lattice gas of noninteracting random walkers (RWs).The latter is described by the more familiar second-order Langevin equation [1] Plugging Eq. ( 21) for F (ρ) into Eq.( 20), we obtain which coincides with the corresponding expression for the gas of noninteracting RWs.Large deviations in the DCG can be described by the MFT [25].The MFT is a weak-noise theory which is based on the Langevin equation ( 19) and relies on a problem-specific small parameter, that we identify below.When such a parameter is present, the probability distribution of interest can be approximately determined via a saddle-point evaluation of the exact path integral corresponding to Eq. ( 19) and the problem-specific constraints.
We will introduce the MFT for the DCG on the example of the statistics of the excess number of particles K on the positive semi-axis at specified time T .One can consider two different settings.In the first of them we start from a constant gas density ρ 0 and condition the process on the value of the integral at time T .This setting corresponds to infinite total mass of the system.In the second setting (a finite total mass) there is a large but finite number of particles, N ≫ 1, in the system.Here the particle excess condition at t = T is where −N/2 ≤ K ≤ N/2.As we will see shortly, the scaling behavior of the particle excess statistics in the two settings is quite different.In the infinite-mass setting the result also depends on whether the initial condition is quenched: that is, deterministically prepared, or annealed: that is, randomly sampled from the equilibrium distribution of the gas at density ρ 0 .

A. Infinite-mass scaling
Let us rescale the variables: t/T → t, x/(ρ 0 DT ) 1/4 → x, and ρ/ρ 0 → ρ.In the new variables Eqs. ( 19) and ( 24) become respectively.At long time the noise becomes effectively weak.The presence of a large parameter (ρ 5 0 DT ) 1/4 ≫ 1, which is a typical number of particles within a region with the length scale (ρ 0 DT ) 1/4 , makes it possible to develop the MFT, that is to perform a saddle-point evaluation of the exact path integral for Eq. ( 26) with account of the constraint (27).The calculation boils down to minimization of the action functional (see Appendix A for details) and leads to the following Hamilton's equations that describe the optimal (that is, the most likely) path of the system conditioned on Eq. ( 27): The boundary condition at the observation time where λ is a Lagrange multiplier, introduced when accommodating the constraint (27), and δ ′ (x) is the xderivative of the delta function.The quenched initial condition is The annealed initial condition is introduced shortly.
Once the optimal path is found, the probability distribution is given, up to a pre-exponent, by the action along the optimal path: − ln P(K, T, ρ 0 ) ≃ S, where where The rescaled excess of the number of particles j is defined in Eq. ( 27).
For the annealed initial condition with average rescaled density 1 the full action includes the cost of creating the optimal initial condition ρ(x, 0): S annealed = S + S 0 , where S is the dynamical action, described by Eqs. ( 32) and ( 33), and As a result, Eq. ( 31) gives way to a different condition [28], which describes a relation between the a priori unknown ρ(x, 0) and v(x, 0):

B. Finite-mass scaling
In this case the parameter ρ 0 is absent, and the rescaling of variables is different: Notice that the infinite-mass dynamical exponent 4 gives way to the finite-mass exponent 5.The rescaled Eqs. ( 19) and ( 25) become and where |j| ≤ 1/2.Here the saddle-point expansion, leading to the MFT, relies on the large parameter N ≫ 1.
The rescaled MFT equations coincide with Eqs. ( 28) and (29).The boundary condition ( 30) is also the same, but Eq. ( 31) gives way to the condition ρ(x, 0) = δ(x) . (39) Here the probability distribution of the number of particles, transferred to the right, is independent of the observation time T : − ln P(K, N ) ≃ S, where and s(j) is again described by Eq. ( 33).

C. Variance of particle excess
The statistics of typical (that is small) fluctuations of the excess number of particles K is Gaussian, and its variance scales as a characteristic number of particles involved.In the infinite-mass case this is the typical number of particles over the dynamical length scale (ρ 0 DT ) 1/4 , that is var K ∼ (ρ 5 0 DT ) 1/4 .In the finite-mass case it is simply var K ∼ N .Essentially, the role of theory (a first-order perturbation theory in |λ| ≪ 1, developed in Ref. [29]) is to provide the numerical coefficients O(1) in these expressions.
For the infinite-mass case the calculations are straightforward.Upon linearization with respect to λ, the MFT equations ( 28) and (29) become We can solve Eq. ( 42) backward in time with the "initial condition" (30).The solution can be obtained by differentiating with respect to x the previously known solution for the initial condition in the form of a deltafunction (see e.g.Ref. [30] and Eq.(B6) in Appendix B below).The result has the similarity form The scaling function V (z) can be expressed via the hypergeometric function 0 F 2 : Now we can evaluate the rescaled dynamical action (33) in terms of λ.Within the linear theory in λ, we should replace the density ρ(x, t) in Eq. ( 33) by 1.We obtain The integral over t gives 4. The integral over z can be evaluated numerically:

Quenched initial condition
For the quenched setting, the formula s = 4αλ 2 suffices for expressing the action in terms of j.Indeed, using the "shortcut relation" ds/dj = λ (see, e.g.Ref. [31]), we obtain which gives λ = j/(8α).As a result, Back to the dimensional variables, we obtain the variance of the typical fluctuations of K: Exactly the same result (48) for the quenched initial condition can be obtained directly from the linearized version of the Langevin equation (26).This calculation is presented in Appendix B. The MFT, however, also enables one to calculate the variance for the annealed initial condition, where a direct calculation with the linearized Langevin equation does not seem to be available.

Annealed initial condition
In the annealed setting one should also take into account (the small-λ expansions of) Eqs. ( 34) and (35), which are the following: respectively, where δρ(x, t = 0) is the a priori unknown small perturbation on the background of the constant density ρ = ρ 0 , and v(x, t = 0) is also small.Evaluating Eq. ( 43) at t = 0, plugging the result into Eq.( 50) and integrating the latter equation twice over x, we determine the optimal initial density field ρ(x, 0) ≃ 1 + δρ(x, 0), where and θ(x) is the step-function.The double integral in this expression can be evaluated with "Mathematica" analytically, but the result is too cumbersome to present it here.
Plugging it into Eq.( 49) and evaluating the resulting integral numerically, we obtain where β = 0.1581 . . . .The total annealed action is, therefore, where α = 0.05798 . . ., as obtained above, and we have again used the shortcut relation ds/dj = λ to express the action through j.At given j, the annealed action (53) is smaller than the quenched action (47).Back in the original variables, the variance of the typical fluctuations of K in the annealed case, is larger than that in the quenched case [see Eq. ( 48)], as to be expected on physical grounds.

IV. MFT OF LARGE DEVIATIONS AT EQUILIBRIUM
Now we suppose that the DCG is at equilibrium and, using the MFT, evaluate the probability density of observing a specified density profile ρ * (x).Importantly, the Hamiltonian MFT equations are still Eqs.( 28) and ( 29) (where we return to the dimensional variables except setting D = 1).The boundary conditions in time, however, become ρ(x, t = −∞) = ρ 0 and ρ(x, t = 0) = ρ * (x) [29,32].
At the microscopic level, the DCG obeys detailed balance.At the macroscopic level the detailed balance manifests itself as the Onsager-Machlup reversibility principle [33].In particular, the optimal activation path ρ(x, t), leading from ρ(x, t = −∞) = ρ 0 to ρ(x, t = 0) = ρ * (x), must coincide with the time-reversed relaxation path from ρ(x, t = −∞) = ρ * (x) to ρ(x, t = ∞) = ρ 0 .That is, the optimal path must obey the equation Combining this equation with Eq. ( 28), we obtain the important relation which describes the equilibrium manifold of this system.The relation (56) brings about two important consequences.Firstly, as one can check by a direct substitution, Eq. ( 29) is now obeyed automatically.Secondly, the mechanical action can be calculated as follows: After two integrations by part over x, this expression becomes By virtue of Eq. ( 55), this expression can be recast as Performing the integration over time, we arrive at the announced Boltzmann-Gibbs relation (23).As a simple but instructive example, let us evaluate the probability distribution P void of observing a void of size 2L in a uniform gas at equilibrium: Using Eq. ( 23), we obtain P void ∼ exp(−2ρ 0 L), as in the gas of noninteracting RWs [32].The optimal path of the system toward the void formation, ρ(x, t) is obtained by the time reversal of the relaxation dynamics of the void, see Eq. ( 55).We computed this optimal path numerically, and the results are shown in Fig. 3. Noticeable are spatial oscillations of the density, which are absent in the optimal path of the void formation in the gas of noninteracting RWs, see Ref. [32].

V. SUMMARY AND DISCUSSION
We have considered some basic macroscopic relaxation and fluctuation properties of the mass-and dipoleconserving stochastic lattice gas.We have extended the MFT approach to this system.Using some carefully selected examples, we have demonstrated how one can apply the MFT both to typical fluctuations of the DPG, and to its large deviations.We hope that the MFT formalism will find additional applications to particular settings relevant to experiments.
We would like to conclude this work with a fascinating observation.Let us briefly recall one basic property of diffusive lattice gases that conserve only the number of particles.The fluctuational hydrodynamics of such gases is described by the Langevin equation [1] which generalizes Eq. ( 22) to a broad class of diffusive lattice gases.The free energy density F (ρ) of this class of gases is determined by the diffusion coefficient D(ρ) and the mobility σ(ρ) via the Einstein relation [1] F ′′ (ρ) = 2D(ρ) σ(ρ) . (62) In the particular case of noninteracting RWs, see Eq. ( 22), one has D(ρ) = D = const, and σ(ρ) = 2Dρ.Then Eq. (62) yields the free energy density F (ρ) = ρ ln ρ − ρ.The remarkable coincidence of this expression with the free energy density (21) for the DCG may suggest the validity of the Einstein relation (62) for dipole-conserving gases at the level of the Langevin equation (26).Indeed, if, by analogy with Eq. (61), we identify the functions D(ρ) = Dρ (which enters the highest-derivative term) and σ(ρ) = 2Dρ 2 , then Eq. ( 62) correctly reproduces the free energy density (21).This observation hints at possible additional surprises in the studies of the dipole-conserving lattice gases.

FIG. 2 .
FIG.2.The scaling function R(ξ) of the similarity solution(11) found by numerically solving the problem (15) and (16) (solid line) and the full time-dependent PDE (4) with a localized initial condition (dashed line).Only the positive-ξ region is shown.