The stochastic relativistic advection diffusion equation from the Metropolis algorithm

Nuclear collisions at high energy exhibit remarkable collective flows, which have been analyzed with considerable success using relativistic hydrodynamics [Heinz:2013th]. For large nuclei, ideal hydrodynamics provides a reasonable description of the observed flows. Viscous corrections are then incorporated by simulating (a version of) the relativistic Navier-Stokes equations, improving the description of the data and clarifying the theoretical consistency of the simulations. These simulations fit the shear viscosity to entropy ratio of QCD around the crossover temperature. Current Bayesian fits give $\eta/s\simeq 2\,\hbar/4\pi k_{B}$ [Bernhard:2019bmu, PhysRevLett.126.202301, PhysRevC.103.054904, Heffernan:2023kpm], which indicates that the medium is remarkably strongly coupled, with relaxation rates of the order $\sim k_{B}T/\hbar$.

The hydrodynamic description of heavy ion collisions can be tested by examining the collisions of light nuclei, such as ${\rm d+Au}$ and ${\rm He}^{3}+{\rm Au}$ at the Relativistic Heavy Ion Collider (RHIC) and proton-nucleus collision at the Large Hadron Collider (LHC). Remarkably, these events also exhibit correlations which are indicative of the collective flow [*[SeeforexampleSec.3.2.2andSect3.2.3in:][.]Arslandok:2023utm]. However, it must be emphasized, that the hydrodynamic description of these events is breaking down. This is in part because (a suitably defined) mean free path $\ell_{\rm mfp}$ has become comparable to the system size, and in part because the total number of particles produced in these events $N_{\rm ch}$ is becoming small, which leads to large fluctuations. One of the motivations for the current paper is to analyze thermal fluctuations in relativistic dissipative systems, with the ultimate goal of describing small colliding systems.

The current manuscript is also motivated by two classical $2^{\rm nd}$ order phase transitions in QCD, which may be an observable in heavy ion collisions. In both of these transitions incorporating thermal fluctuations into the hydrodynamic description is essential to describing the underlying physics in the critical region. The first phase transition is the $O(4)$ chiral transition of QCD, which is a $2^{\rm nd}$ order phase transition in the limit of two massless quark flavors [Pisarski:1983ms, Rajagopal:1992qz]. There is a strong motivation from lattice QCD to look for signatures of the $O(4)$ transition in the heavy ion data at the LHC [HotQCD:2019xnw, Kaczmarek:2020sif]. At lower temperatures and higher baryon density, strong theoretical arguments suggest that there should be a Ising critical point in the $(T,\mu_{B})$ plane [Stephanov:2004wx], with $T$ and $\mu_{B}$ being the temperature and the baryon chemical potential, respectively. Currently at RHIC, there is an ongoing search for the Ising critical point where the beam energy is scanned in an effort to tune the baryon chemical potential to the critical region of the phase diagram [Du:2024wjm].

We are also motivated to study thermal fluctuations in relativistic fluids by algorithmic developments and the mathematical structure of stationary stochastic processes. In statistical mechanics, the Metropolis algorithm is used to generate field configurations of a field $\phi$ from a known probability distribution, $P[\phi]\propto e^{S[\phi]}$, where $S[\phi]$ is the entropy [Landau_Binder_2014]. In the algorithm a proposal is made for a change in the fields, $\phi\rightarrow\phi+\Delta\phi$. This proposal is accepted or rejected according to the magnitude and sign of the change in the entropy, $\Delta S\equiv S[\phi+\delta\phi]-S[\phi]$. The Metropolis steps are guaranteed to converge to the required equilibrium distribution. Recently, in the context of simulating the $O(4)$ critical point, we simulated the stochastic diffusion of a conserved charge coupled to the order parameter using a variant of the Metropolis algorithm [Florio:2021jlx]. A proposal is made for the transfer of charge between the fluid cells, and this proposal is then accepted or rejected based on the change in the entropy of the system. For small enough time steps the Metropolis updates naturally reproduce the Langevin dynamics of the diffusion equation. The equivalence of the Metropolis and Langevin dynamics for small time steps $\Delta t$ has been repeatedly observed over the years [OldMCMC, Moore:1998zk].

The advantages of a Metropolis based approach is that detailed balance and the Fluctuation-Dissipation-Theorem (FDT) are automatically preserved, independently of $\Delta t$, which guarantees that Markov-chain will equilibrate to a specific action. In non-linear theories this simplifies the renormalization of the theory and clarifies discretization ambiguities that arise in non-linear Langevin equations [Arnold:1999uza, Gao:2018bxz]. For fluids at rest, the Metropolis algorithm has been used to implement the non-linear Langevin dynamics in a number of challenging applications, leading to the study of sphaleron transitions of hot QCD [Moore:1998zk], the real time dynamics of $O(4)$ critical point in QCD [Florio:2021jlx, Florio:2023kmy], and the dynamics of Model B and Model H [Chattopadhyay:2023jfm, Chattopadhyay:2024jlh] in the Halperin-Hohenberg classification of dynamical critical phenomena.

Our principal task in this and a companion paper is to generalize the Metropolis approach to relativistic fluids in general coordinates, where the (somewhat complicated) form of the dissipative stress should arise naturally from the accept/reject steps of the Metropolis algorithm. As a first step, in this paper we will consider the diffusion of charge in a relativistic fluid.

To understand the issues that arise with relativity, consider the relativistic advection-diffusion equation in flat spacetime in the Landau-Lifshitz frame, momentarily neglecting the stochastic noise for simplicity. We are considering a fluid moving with three velocity $v^{i}$ in the lab frame and following the diffusion of a dilute conserved charge within the fluid, $\partial_{\mu}J^{\mu}=0$. The four velocity is $u^{\mu}=(\gamma,\gamma{\bf v})$ and the local charge density in the rest frame of the fluid is $n_{\rm\scriptscriptstyle LF}=-u_{\mu}J^{\mu}$. Landau and Lifshitz define a hydrodynamic frame where the chemical potential is given by the value of $n$ and ideal equation of state to all orders in the gradient expansions, i.e. $n_{\rm\scriptscriptstyle LF}=\chi\,\mu_{\rm\scriptscriptstyle LF}$ where $\chi$ is the charge susceptibility [landau2013fluid].

The problem with the covariant Landau-Lifshitz approach is that the diffusive current, which was spatial in the rest frame of the fluid $j_{D}^{i}=-D\partial^{i}n$, involves time derivatives in an arbitrary frame. This leads to equations which are second order in time, which in turn leads to runaway solutions and other pathological behavior [Hiscock:1985zz, Gavassino:2021owo]. One way to correct this pathology is to promote the diffusive current to an additional dynamical field which relaxes on collisional timescale to the expected form. This procedure results in Maxwell-Catteneo or Israel Stewart type equations.

There is merit to the Israel-Stewart approach – it provides an effective way to realize the dynamics of the relativistic diffusion equation and the relativistic Navier Stokes equations more generally. Indeed, almost all large scale simulations of flow in heavy ion collisions are based on this approach. However, the Israel-Stewart formulation involves fast variables whose physical significance should be questioned. Indeed, there have been many reformulations of viscous hydrodynamics in the relativistic domain [Israel:1976tn, Israel:1979wp, Geroch:1990bw, OTTINGER1998433, Bemfica:2017wps, Kovtun:2019hdm], and each of these reformulations involve some additional variables (or non-hydrodynamic modes) which relax quickly to the form constrained by “first-order” hydrodynamics [Lindblom:1995gp]. A kind of theorem has emerged, which states that it is impossible to construct a causal and stable relativistic theory of hydrodynamics without incorporating non-hydrodynamic modes [Heller:2022ejw, GAVASSINO2023137854, Gavassino:2023mad].

In this paper we will investigate an alternative to the Israel-Stewart approach developed to describe fluids without an underlying boost symmetry [Novak:2019wqg, deBoer:2020xlc, Armas:2020mpr]. In particular, we found the “density frame” discussed by Armas and Jain a clarifying formalism when implementing Metropolis updates [Armas:2020mpr]. The density frame has no non-hydrodynamic modes and no additional parameters compared to the Landau theory of first order hydrodynamics, at the price of not being fully boost invariant. In hydrodynamics without boosts the constitutive relations are written down for setups where the underlying interactions are not Lorentz invariant and thus there is a preferred “lab” frame¹¹1For a physical example, consider a fluid flowing over a table and study the diffusion of charge in this background fluid flow. The table sets a preferred lab frame.. If the microscopic interactions happen to be Lorentz invariant, the additional boost symmetry imposes relations between the coefficients of the gradient expansion, which is an expansion in lab-frame spatial derivatives. Indeed, the equations of motion in the density frame follow from the Landau ones if the ideal equations are used to rewrite lab-frame time derivatives appearing in the dissipative strains as spatial derivatives. With this rewrite the equations are strictly first order in time and are stable.

The procedure amounts to a non-covariant choice of hydrodynamic frame where the relation between the chemical potential and the lab frame charge per volume $J^{0}$ is given by ideal hydrodynamics at all orders in the gradient expansion, i.e. $J^{0}=\chi\mu u^{0}$. The frame choice and the resulting equations of motion in the density frame are not invariant under Lorentz transformations; but, they are invariant under Lorentz transformations followed by a change of hydrodynamic frame, which reparametrizes the hydrodynamic fields. The results obtained in different Lorentz frames will vary, but the variation is beyond the accuracy of the diffusion equation. Different choices of Lorentz invariant hydrodynamic frames, such as the Landau, Eckart or BDNK²²2The acronym is short for Bemfica, Disconzi, and Noronha [Bemfica:2017wps] and Kovtun [Kovtun:2019hdm]. These authors studied a class of Lorentz invariant frames where the temperature is related to the Landau choice up to derivatives. choices, will also give different (if Lorentz invariant) results to this accuracy. The density frame approach was known to “work” in simple cases, but the work on hydrodynamics without boosts formalized the procedure. Finally, Armas and Jain made an important connection of this approach to modern treatments of hydrodynamics [Crossley:2015evo, Glorioso:2017fpd] where the conserved charge and the corresponding canonical conjugates (for instance a $U(1)$ charge $Q$ and the associated phase $\varphi$) play a dual role [Armas:2020mpr]. The numerical utility of the symplectic structure inherent in this duality remains to be fully exploited.

The structure of the current paper is as follows. Section II discusses the relativistic advection-diffusion equation in the Landau and density frames, and derives the density frame constitutive relation from covariant kinetic theory. Section III compares the density frame relativistic advection-diffusion equation with relativistic kinetics. As discussed above, the density frame is not Lorentz invariant, but it is invariant under Lorentz transformations followed by a reparametrization of the hydrodynamic variables. In a specific test problem discussed in Section III, we show that, in the regime of validity of hydrodynamics, the deviations of the density frame from a underlying covariant microscopic theory are controllably small, even for highly boosted fluids. We also study the convergence of the gradient expansion of the density frame, making connections with Lorentz covariant approaches. Stochastic dynamics in the density frame is studied in Section IV. Since the hydrodynamics equations are defined using a given foliation of space-time, the Metropolis updates discussed above are very natural and easy to implement. One simply makes random transfers of charge in between ideal advective steps. These charge transfers are then accepted or rejected using the entropy defined on the spatial slice as the statistical weight. We present a first study of numerical correlation functions from the Metropolis algorithm for an equilibrated boosted fluid in Section IV. Although we have used the Metropolis algorithm for the relativistic hydrodynamics in the density frame, it should be useful for other approaches to stochastic relativistic hydrodynamics, e.g. approaches based on BDNK [Mullins:2023tjg, Gavassino:2024vyu] or Israel-Stewart [Singh:2018dpk]. Finally, in Section V we conclude with a short discussion of the next steps.

We are considering the advection and diffusion of a charge in fluid moving at relativistic speeds in flat spacetime, $\eta_{\mu\nu}=(-,+,+,+)$. The charge density is low and the temperature and flow velocity $u^{\mu}=(\gamma,\gamma{\bf v})$ may be considered fixed. (The generalization of this problem to a time dependent background fluid is discussed in App. A.) In first order hydrodynamics in the Landau frame the conserved current obeys [landau2013fluid]

where the first term in $J^{\mu}$ is ideal advection and the second term is the diffusive correction, expressed as

Here $T$ is the temperature, $\sigma$ is the conductivity, $\hat{\mu}_{{\rm\scriptscriptstyle LF}}\equiv\mu/T$ is the scaled chemical potential, thermodynamically conjugate to charge density $n_{\rm\scriptscriptstyle LF}$. $\Delta^{\mu\nu}$ is the spatial projector

and satisfies $\Delta^{\mu\rho}\Delta_{\rho\nu}=\Delta^{\mu}_{\phantom{\nu}\nu}$. Since the density is low, $n_{{\rm\scriptscriptstyle LF}}=\chi\mu_{{\rm\scriptscriptstyle LF}}$ where $\chi$ is the (temperature dependent) susceptibility.

In the next subsections we will briefly review the density frame pointing out the differences with the Landau frame. To keep the presentation self contained and pedagogical, sections II.2 and II.3 review a small portion of [Armas:2020mpr] with a focus on the diffusion equation. Section II.4 describes how the density frame constitutive relation arises in relativistic kinetic theory.

Consider a portion of a fluid in perfect global equilibrium within a lab frame measurement volume, $V_{0}\equiv\int{\rm d}\Sigma_{0}$. The entropy, energy-momentum, and charge on this slice are

We will notate the charge density with $N\equiv J^{0}$, and the temporal components of the energy momentum tensor with $(E,M^{i})\equiv(T^{00},T^{0i})$, reserving the calligraphic symbols $\mathcal{S}$, $\mathcal{P}_{\mu}$ and $\mathcal{N}$ for the charges, as opposed to the charge densities, $S,E,M,N$. When a boost symmetry is ultimately adopted, the entropy will take the form $S=su^{0}$ where $s$ is a Lorentz scalar.

Given the conserved charges $\mathcal{P}_{\mu}$ and $\mathcal{N}$, the temperature, velocity and chemical potential can be determined from the micro-canonical equation of state $S(\mathcal{P},\mathcal{N},V_{0})$

where $\hat{\mu}\equiv\mu/T$ and $v^{i}\equiv\beta^{i}/\beta^{0}$ is the fluid velocity. The Gibbs-Duhem relation follows from the extensivity of the system

Then for small charge densities the entropy density as a function of $N$ takes the form

So far we have not used boost symmetry, and the parameters, such as $\beta^{0}$, $\hat{\mu}$ and $\chi^{00}$ are functions of $E$ and $M$. After imposing boost symmetry in the next paragraph, $\beta^{0}$ will be a Lorentz scalar $\beta$ times $u^{0}$, and $\chi^{00}$ will be a Lorentz scalar $\chi(\beta)$ times $u^{0}u^{0}$

justifying the notation a posteriori.

When the fluid has an underlying Lorentz symmetry the dependence on the velocity is determined by the symmetry. Before imposing the symmetry, the conservation laws take the form

At zeroth order in the gradient expansion (ideal hydrodynamics) the three current $J^{i}$, the energy flux $T^{i0}$, and the spatial stress tensor $T^{ij}$ are algebraically related to the conserved charges $E$, $M^{i}$, and $N$. Using the conservation laws, the thermodynamic relations (eqs. (5) and (7)), the symmetry of the stress tensor imposed by Lorentz invariance $T^{i0}=M^{i}$, and requiring that $\partial_{t}S+\partial_{i}S^{i}$ be non-negative, fixes the form of fluxes $J^{i}$, $M^{i}$, $T^{ij}$ and entropy current $S^{i}$ to the form of ideal hydrodynamics [Armas:2020mpr]. In particular, the charges and entropy are parametrized by $\beta^{\mu}\equiv\beta u^{\mu}$ and the chemical potential $\hat{\mu}$

Here $e$, $p$, $n$ and $s$ are scalar functions of $\beta\equiv\sqrt{-\beta^{\mu}\beta_{\mu}}$ and $\hat{\mu}$ as given by the equilibrium equation of state. At small chemical potentials, the local charge density and entropy takes the form

where $\chi(\beta)$ is a function of temperature. Comparing the definitions in eq. (8) and eq. (12), we find $\chi^{00}/\beta^{0}=T\chi u^{0}$ as claimed above.

In the density frame the familiar relations of ideal hydrodynamics, eqs. (11), serve to define the temperature, chemical potential, and flow velocity in terms of the lab frame charges, $E$, $M^{i}$ and $N$, at every order in the gradient expansion, i.e. the relation between the charges $E,M^{i},N$ and their conjugates do not receive viscous corrections. In particular, the chemical potential in the density frame is defined at all orders in the gradient expansion as

This definition should be contrasted with the chemical potential in the Lorentz invariant Landau frame where

With the density frame definition a fiducial observer needs to count the charge in a given measurement volume in order to determine the chemical potential. With the Landau frame definition, the three current $J^{i}$ also needs to be measured. Thus, the Landau frame involves counting the charges at different times in order to define the chemical potential.

Here we will derive the density frame equations of motion by first considering fluids without a boost symmetry, and then specializing the equations to Lorentz covariant fluids. An example of a two dimensional non-Lorentz invariant fluid is a fluid flowing over a flat surface at relativistic speeds. The diffusion of a charge in this fluid depends on the speed of the fluid relative to the surface.

The advection-diffusion equation in the density frame consists of the conservation law

together with a constitutive relation for the diffusive current $J_{D}^{i}$

The diffusive current is expanded in spatial gradients of the conserved charge, or its thermodynamic conjugate $\hat{\mu}$. The most general form of $J_{D}^{i}$ at first order in gradients of $\hat{\mu}$ is

where $\hat{\mu}\equiv\partial S/\partial N=\beta^{0}N/\chi^{00}$ and $\hat{v}^{i}=v^{i}/|v|$ is a flow unit vector. The first and second terms on the right hand side of eq. (17) capture the diffusion parallel and perpendicular to the fluid motion respectively. Demanding that entropy production be positive leads to the requirement that $\sigma_{\parallel}(\beta^{0},v)>0$ and $\sigma_{\perp}(\beta^{0},v)>0$, but no further constraints can be derived on general grounds.

Lorentz invariant fluids can be treated as a special case of eq. (17). Indeed, the boost symmetry determines the dependence of $\sigma_{\parallel}$ and $\sigma_{\perp}$ on the velocity. The easiest way to derive this relation is to return momentarily to the Landau frame. Comparison to the density frame form gives

and

We now use the lowest order equation of motion³³3When the fluid velocity is not constant this expression receives additional corrections proportional to derivatives of velocity. The generalized constitutive relation is given in App. A.,

to approximate the Landau frame expression for the diffusive current

Substituting eq. (21) into eq. (19) gives the current in the density frame

where we have defined a frequently occurring matrix

A generalization of the constitutive relation in (21) for fluids depending on space and time is given App. A.

Comparison with the general form in eq. (17) shows that

In summary, in the density frame the equation of motion is

and when $\hat{\mu}$ is written in terms of the charge $N=\chi\mu u^{0}$, we arrive at an advection-diffusion equation

with a diffusion matrix

Here $D=\sigma/\chi$ is the scalar diffusion coefficient of the Landau frame. Apart from the tensor structure the equation is numerically similar to the non-relativistic advection-diffusion equation and can be solved by familiar numerical methods⁴⁴4In our numerical tests without noise we used an IMEX scheme with a standard advection step and a Crank-Nicholson like implicit step using the PETSc library [petsc-web-page]. .

The $\gamma$ factors in the diffusion matrix can be easily understood physically. The diffusion coefficient has units of distance squared per time. The rate of transverse diffusion is suppressed relative to a fluid at rest by one factor of $\gamma$ due to time dilation. The rate of longitudinal diffusion is suppressed by three factors of $\gamma$ due to time dilation and length contraction, i.e. each spatial step in the random walk is length contracted by $\gamma$ and the steps add in square.

In this subsection we will briefly describe how the density frame constitutive relation arises naturally in relativistic kinetic theory. Specifically, we will show how the conductivity matrix $T\sigma^{ij}$ follows from covariant kinetics and find how this matrix is determined by the particles through the first viscous correction $\delta f$ to the phase-space distribution function. For simplicity, we will assume a relaxation time approximation and consider single species of classical relativistic particles, which carry the charge of the system

Here ${\mathcal{C}}_{p}$ is a momentum dependent parameter controlling the collision rate in the rest frame of the medium.

In global equilibrium the phase space distribution function is characterized by a constant chemical potential, temperature, and flow velocity. If the density of the charged particles depends slowly on space and time then the parameter $\hat{\mu}(t,{\bf x})$ is no longer a constant but reflects this dependence

In the density frame $\hat{\mu}(t,{\bf x})$ is adjusted to reproduce the charge density in the lab frame $J^{0}$, while in the Landau frame $\hat{\mu}(t,{\bf x})$ is adjusted to reproduce the charge density in the rest frame, $n(t,{\bf x})=-u_{\mu}J^{\mu}$. These two definitions of the chemical potential agree when gradients are neglected, and in this case $f_{\rm eq}(t,{\bf x},{\bf p})$ is a solution to the Boltzmann equation. $\hat{\mu}(t,{\bf x})$ obeys the equations of ideal advection equation at lowest order

Following a standard procedure to find the first viscous correction [Teaney:2009qa], we parameterize $f=f_{0}+\delta f(t,{\bf x},{\bf p})$ and solve for $\delta f$ order by order in the gradients

In the Landau frame one decomposes the gradient into its temporal and spatial components as

We neglect the temporal term in Eq. (33) exploiting the lowest order equations of motion, Eq. (31). Then we substitute into Eq. (32), which leads to the familiar form of the first viscous correction in the Landau frame

where $\nabla_{\alpha}=\Delta_{\alpha}^{\phantom{\nu}\mu}\partial_{\mu}$. Evaluating the diffusive current

yields the expected form of the current in the Landau frame

The conductivity in this expression is defined from the transport integrals

In the density frame one proceeds similarly, but uses the lowest order equations in the lab frame

which yields the form of the first viscous correction

The diffusive current $J_{D}^{i}$ in the density frame is the difference between the current and the ideal advection $J^{0}v^{i}$

Substituting the approximate distribution function $f_{0}+\delta f$ leads to an appealing positive definite symmetric matrix which determines the mean current⁵⁵5As discussed below, the matrix $\mathcal{K}^{ij}$ also determines functional form of the noise in the density frame. In Ref. [Armas:2020mpr] the matrix $\mathcal{K}$ is written as $D_{jj}^{\rho\sigma}$. This section shows how this form arises in a microscopic theory.

Noting that

we find that $\mathcal{K}^{ij}$ has the expected density frame form

where we used the integrals defined in eq. (37). To summarize this section, we have shown how the form of the dissipative current in the density frame arises in covariant kinetic theory (eqs. (41) and (43)). This form is not covariant although the underlying kinetic theory is covariant. This arises because we are trying to approximate the full results of kinetic theory in a specific frame.

In this and the next subsections we will study an analytically tractable covariant kinetic model in $1+1$ dimensions and investigate how the current in this model approaches the density frame constitutive relation. The model consists of massless “particles” moving with the speed of light along a one-dimensional line. The particles experience Poissonian random kicks which changes the direction of their velocities. The dynamical evolution of this system can be mapped to the famous telegraph equation from which an all-orders gradient expansion can be derived. We will then compare the results of this all-order gradient expansion with the predictions from the density frame formalism.

In this subsection we will analyze this model in a static background case and then in the next subsection consider a fluid background moving with constant velocity. The results in this subsection are not new and can be found in many places, but they will serve to define the terms for the boosted case.

Let us denote the respective densities of the left and right moving particles as $n_{-}$ and $n_{+}$ (see Fig. 1).

Due to the random kicks, the particles change direction with the rate $\Gamma\equiv 1/2\tau_{R}$ where $\tau_{R}$ is the relaxation time. The density of the left/right movers obey the following kinetic equation

Adding eqs. (44) and (45) and rearranging leads to the conservation equation

where $n=n_{+}+n_{-}$ is the density and $j=c(n_{+}-n_{-})$ is the current. Similarly, subtracting eqs. (44) and (45) and rearranging gives the relaxation equation for the current

As usual, this set of first order equations can be expressed as a second order equation for $n$

which is known as the Telegraph equation. The exact Green functions associated with this system are known analytically and can be expressed in terms of modified Bessel functions which are given in the Appendix B for convenience. For simplicity, we will set $c=1$ for the remainder of this section.

Consider some initial state at $t=0$ (not necessarily in equilibrium) that is specified by two independent functions $n(t=0,x)=n_{0}(x)$ and $j(t=0,x)=j_{0}(x)$. Let us further assume that the initial current can be expressed as a gradient expansion of the initial density, with some coefficients $j_{0}=-\sum_{\ell=0}^{\infty}\tau_{R}^{2\ell+1}b_{2\ell+1}\partial_{x}^{2\ell+1}n_{0}(x)$. Note that the coefficients $b_{2\ell+1}$ partially characterize the initial state and they are dimensionless by definition.

By using the exact Green functions, it can be shown that at late times $t\gtrsim\tau_{R}$ the system obeys a universal dispersion relation

where

with $C(\ell)=(2\ell)!/(\ell!(\ell+1)!)$ being the $\ell^{\rm th}$ Catalan number and $P_{\ell,m}(t)$ is a polynomial in $t$ of order $\ell-m+1$ whose coefficients depend on $b_{i\leq\ell}$. It is clear from this expansion that the dependence on the initial conditions decays exponentially for $t\gtrsim\tau_{R}$, meaning that the system thermalizes and obeys a universal constitutive relation. The constitutive relation can be expressed as an all-order gradient expansion whose coefficients are given by $(-1)^{\ell}C(\ell)$. The Catalan numbers are nothing but the Taylor coefficients of the dispersion curve $\omega(k)$ of the diffusion mode:

associated with the linear system expanded around $k=0$. The dispersion relation has a branch singularity at $k_{*}=\pm 1/(2\tau_{R})$. Furthermore the large $k$ expansion is consistent with the stability and causality conditions given in Refs. [Heller:2022ejw, GAVASSINO2023137854, Gavassino:2023mad, Hoult:2023clg] where the Telegraph equation was previously analyzed in these terms.

The additional eigenfrequency from (51) is gapped and approaches $\omega(k)=-i/\tau_{R}+\mathcal{O}((\tau_{R}k)^{2})$ for $k\rightarrow 0$. The gapped modes are responsible for the exponential decay in eq. (50). As we will see in the next section, when the fluid is moving with velocity $v$, the density frame constitutive relation is approached on a short timescale of order $\tau_{R}/\gamma$, which again reflects the dynamics of gapped non-hydrodynamic modes.

Let us now consider the same kinetic model of the previous section in the background of a fluid moving uniformly with velocity $v$ in the lab frame, as shown in fig. 2.

In the lab frame, the density of left and right movers obey

where transition rates are

The appearance of the kinematical factors $\kappa_{\pm}=\sqrt{(1\pm v/c)/(1\mp v/c)}$ can be understood in the following way: the right movers in the local rest frame, denoted by the subscript $R$, follow the trajectory $x_{R}=ct_{R}$. A Lorentz boost to the lab frame coordinates, $(t,x)$, leads to the time dilation factor $t=\gamma(1+v/c)t_{R}=\kappa_{+}t_{R}$ where $\gamma=(1-v^{2}/c^{2})^{-1/2}$ is the Lorentz factor. Since the mean free path time in the local rest frame is $\Gamma^{-1}=2\tau_{R}$, the mean free path time in the lab frame is therefore $\Gamma^{-1}_{+}=\kappa_{+}2\tau_{R}$.

Similarly to the static case, by adding and rearranging equations (52) and (53), we get

where $N=N_{+}+N_{-}$ and $J=c\,(N_{+}-N_{-})$ denote the density and the total current respectively. Following the density frame approach, we write the current as a sum of the advective and diffusive parts,

The diffusive part satisfies the relaxation equation

For clarity we will set $c=1$ for the rest of this section.

Let us now consider an initial state given at $t=0$ in the lab frame. Based on our findings in the static case, we expect the system to lose information about the initial conditions for $t\gtrsim\tau_{R}/\gamma$ and obey a universal gradient expansion

Note that due to the nonzero velocity that explicitly breaks parity, we expect both even and odd terms in the derivative expansion as opposed to the static case, which only contains odd terms. We factored out an overall factor of $\gamma^{-2}$ and the characteristic time scale $\tau_{R}/\gamma$ explicitly to simplify the expressions for $c_{n}$. As in the static case, the coefficients of the gradient expansion, $c_{n}$, can be calculated from the dispersion relation

The branch singularity in this case is in the complex plane, $k_{*}=(\pm 1-iv\gamma)/(2\tau_{R})$. Note that the radius of convergence, $|k_{*}|=\gamma/(2\tau_{R})$, grows with $\gamma$ so that the hydrodynamic description applies to modes of wavenumbers $k\lesssim\gamma/\ell_{\rm mfp}$ in the lab frame. Here and below we define the rest-frame mean-free-path:

The other branch in the dispersion relation is gapped, $\omega(k)=-i/(\tau_{R}/\gamma)$ for $k\rightarrow 0$. This branch describes a non-hydrodynamic mode decaying exponentially on a timescale of $\tau_{R}/\gamma$. Notably, the decay time is not time-dilated, but instead time-contracted by a factor of $\gamma$ relative to the static case.

The $c_{n}$ follow from Taylor expanding $\omega(k)$ in eq. (61) around $k=0$. The even and odd terms are found as

For reference we write down the first few terms below:

Just like the static case, the large $k$ expansion of the dispersion relation given in eq. (61) necessarily satisfy the causality and stability conditions of Refs. [Heller:2022ejw, GAVASSINO2023137854, Gavassino:2023mad, Hoult:2023clg].

Putting everything together, we find the gradient expansion in the lab frame as

Let us discuss this result. The leading term in eq. (66) agrees with the prediction from the advection-diffusion equation in the density frame. In some sense, the higher order corrections encode the underlying microscopic dynamics and they restore Lorentz causality as can be seen from the large $k$ expansion, for instance. Of course Lorentz causality is violated if the gradient expansion is truncated at any finite order.

One might ask when the physics of the covariant kinetic model is captured by the density frame diffusion framework. We first require that the system reaches local thermal equilibrium and therefore can be described by the universal constitutive relation. This happens in a timescale of order $t\sim\tau_{R}/\gamma$, which is set by the inverse frequency of the gapped modes. This is illustrated in Fig. 3 where at $t=0$ we initialize a

Gaussian drop of charge, $N_{0}(x)=N_{\rm max}\exp(-x^{2}/2\sigma^{2})$, and set the initial diffusive current to zero in the kinetic model, $J_{D}=0$. For the test case shown we took a fluid with $\gamma=10$ and set $\sigma=L/\gamma$ with $L/\ell_{\rm mfp}=50$. The left figure shows the time evolution of the charge density in the lab frame. In the short period of time considered in the figure $\Delta t\sim\tau_{R}/\gamma$, the charge does not have enough time to diffuse and it is simply advected by the background flow. The right figure shows the evolution of the lab frame diffusive current over the same time period. Clearly the diffusive current in the kinetic model relaxes on a time of order $\tau_{R}/\gamma$ to a steady state, and at $t=5\tau_{R}/\gamma$ the steady state agrees with the leading derivative term of density frame gradient expansion given in eq. (66) to a very good precision.

The relaxation timescale $\tau_{R}/\gamma$ is easily understood. Indeed, the mean collision times of right and left movers are of the order $\gamma\tau_{R}$ and $\tau_{R}/\gamma$ respectively. In equilibrium, where $J=Nv$, the number of right and left movers is

and thus for $v\rightarrow 1$ there are almost $N$ right movers and of order $\sim N/\gamma^{2}$ left movers in the equilibrium sample. If all the particles are initially left movers, then in a time of order $\tau_{R}/\gamma$, these initial particles will scatter with probability one, reaching approximate equilibrium with $N_{+}\simeq N$. Similarly, if all the particles are initially right movers, then in a time of order $\tau_{R}/\gamma$, $N/\gamma^{2}$ of these initial particles will scatter, generating a yield of left movers commensurate with equilibrium, $N_{-}\sim N/\gamma^{2}$. Thus, the typical equilibration time is of order $\tau_{R}/\gamma$.

Of course, the applicability of the gradient expansion depends on the size of the system as well, and we expect it to break down for smaller systems. In Fig. 4 we show the diffusive current for the same setup as above with

$\gamma=10$ at $t=5\tau_{R}/\gamma$, but with different initial Gaussian sizes ranging from $L=(1\ldots 8)\,\ell_{\rm mfp}$ in the rest frame. From the figure we see that for a small system when $L$ is one mean-free-path, the gradient expansion does not converge at all. As the system size gets larger, the gradient expansion becomes more and more accurate. Note that the leading term in the gradient remains qualitatively well behaved even for very small systems. Remarkably, the convergence starts already when the system size is only $3$ or $4$ mean-free-paths long, where the second and third order terms in the gradient expansion capture the physics quite accurately. Notably in this intermediate region, the effect of the second order term is clearly seen by comparing the red (exact) and blue (first order) curves in the lower left figure. This skewness is due to the background motion, and is absent in the static case where parity is unbroken.

The convergence of the gradient expansion can be analyzed more rigorously. Let us consider the time interval $\tau_{R}/\gamma\ll t\ll\tau_{R}\gamma^{3}$ where the system is in local thermal equilibrium but the initial wave packet has not diffused yet. In this regime $N(t,x)\approx N_{0}(x)$, therefore, the gradient expansion, eq. (59), becomes

where $H_{n}(x)$ denotes the $n^{th}$ Hermite polynomial that follows from the gradient expansion of the Gaussian wave-packet with the width $\sigma$. We do not keep track of the overall magnitude of the initial density as it is inconsequential for our argument. It is straightforward to observe that the coefficients, $c_{n}$, grow exponentially in $n$ ⁶⁶6This can be seen from the dispersion relation, eq. (61), whose Taylor coefficients around $k=0$ are $c_{n}$.. Given that asymptotically, the Hermite polynomials grow factorially as $H_{n}\sim\Gamma(n/2)$ (apart from the finite number of points where they vanish), the gradient expansion is an asymptotic expansion. The effective coupling constant of this asymptotic series is $\tau_{R}/(\sigma\gamma)$. A well known property of an asymptotic series that grow as $g^{n}\Gamma(n/2)$ is that it starts to diverge at order $n^{*}\sim 2/g^{2}$. According to the optimal truncation procedure a la Poincaré, the gradient expansion can be directly summed up to $n^{*}$ terms before the series start to diverge. This result implies that when the system size is comparable to mean-free-path, namely $\sigma=c\tau_{R}/\gamma$ for some constant $c$, we get $n^{*}=1$ and the optimal truncation breaks down.

Furthermore in the ultra-relativistic limit, $v\rightarrow 1$, the gradient expansion coefficients significantly simplify; $c_{n}=2^{n-1}+{\cal O}(\gamma^{-1})$ and we find the asymptotic behavior of the gradient expansion to be

Therefore the breakdown of the optimal truncation occurs when $\sigma=2\tau_{R}/\gamma$, confirming our heuristic argument earlier.

In this section we will add noise to the density frame diffusion equation and study the stochastic dynamics of a boosted fluid. The dissipative current in the density frame take the form

where $\xi^{i}$ is the noise. For the stochastic process to equilibrate to the probability distribution determined by the entropy of the system,

the noise must respect the fluctuation-dissipation theorem

In App. A we describe how the noises added to the system should be generalized when the fluid velocity depends on space and time.

The form of the noise matrix in the density frame can also be found by algebraically manipulating the current in the Landau frame. In the Landau frame, we have