Linearized regime of the generalized hydrodynamics with diffusion

We consider the generalized hydrodynamics including the recently introduced diffusion term for an initially inhomogeneous state in the Lieb-Liniger model. We construct a general solution to the linearized hydrodynamics equation in terms of the eigenstates of the evolution operator and study two prototypical classes of initial states: delocalized and localized spatially. We exhibit some general features of the resulting dynamics, among them, we highlight the difference between the ballistic and diffusive evolution. The first one governs a spatial scrambling, the second, a scrambling of the quasi-particles content. We also go one step beyond the linear regime and discuss the evolution of the zero momentum mode that does not evolve in the linear regime.

In this work, we consider the resulting dynamics for the Lieb-Liniger model, a gas of 1d bosons with ultra-local interactions. The GHD forms a complicated, infinite, set of non-linear equations. Our aim is to analyze solutions to those equations and display some of their qualitative features, especially those relevant to the question of equilibration. Specifically, we focus on the linearized regime of the GHD in which the resulting equations take a form of an infinite set of integro-differential equations.
The manuscript is organized in the following way. After the introductory part of the next section, follows three sections with results. In Section 3 we display a general solution to the linear GHD in terms of the eigenstates of the diffusion operator. We also study numerically the properties of this operator and conjecture their effect on the dynamics. In the following Section 4 we solve numerically the GHD equations for two prototypical initial states: the localized and the delocalized state. In Section 5 we consider a next order correction to the linear GHD and study the dynamics of the k = 0 spatial mode which does not evolve in the linear approximation. Finally, the Appendix A is devoted to set the notation for integral operators.
2 Generalized hydrodynamics, diffusion and the Lieb-Liniger model We start this section by presenting the Generalized Hydrodynamics following the original works [41,43]. Then, we recall the Lieb-Liniger model and introduce the ingredients necessary for its GHD description. We will consider the Lieb-Liniger model with repulsive interactions, which contains quasi-particles of only one type. Below, we present the GHD appropriate for this case. Generalization to situations with multi-species quasi-particles is possible [43].

Generalized hydrodynamics
The integrable theories are characterized by stable quasi-particles. The stability of the quasiparticles is related to the presence of an extensive number of local (or quasi-local) conserved densitiesq i (x, t). The hydrodynamic picture relies on the assumption that the state of the system can be locally described by the averagesq i (x, t) = q i (x, t) . These averages can be conveniently parametrized introducing local distribution of the quasi-particles ρ p (θ, t, x) where h i (θ) is the one-particle eigenvalue of the conserved charge with density q i (x, t) and rapidity θ parametrizes the quasi-particle. At the Euler scale ρ p (θ, x, t) corresponds to the density of a local Generalized Gibbs Ensemble [44][45][46]. Beyond that scale, including the next term in the hydrodynamic derivative expansion, the average values of generic observables depend not only on ρ p (θ, x, t) but also on its spatial derivatives.
At the Navier-Stokes level the dynamics of the GHD can be understood through the following qualitative picture. When a quasi-particle travels through an inhomogeneous system it experiences two effects. First, the changing distribution of the quasi-particles leads to a difference in the propagation velocity. Second, the scattering processes with other particles lead to diffusion, see Fig. 2. These two effects are combined in Navier-Stokes-like equations of GHD with the effective velocities v eff n and diffusion operator D n depending on the local (in the space-time) distribution of quasi-particles. Solving (2) is difficult even in the absence of the diffusion term, however, in some cases the problem can be rephrased into more trackable integral equations [31]. The difficulty lies in the fact that each mode θ propagates with a velocity v eff n (θ) depending non-linearly on the local density function ρ p (θ, t, x) of quasiparticles. Therefore, already at the ballistic level, the dynamics is complicated and non-linear [47] and displays the equilibration phenomena. The known results describe, for example, the self-similar solution [27] or numerical solution [38] for the nonequilibrium steady states.
Inclusion of the diffusion term in (2) further complicates the situation, e.g. it leads to mixing between densities of different rapidities θ. Our aim is to understand better the role of the diffusion term in the dynamics governed by (2).
Let us specify now in some details the ingredients entering Eq. (2). The effective velocity The functions (E ) dr (θ) and (p ) dr (θ) are dressed derivatives of the effective energy and momentum of a particle with rapidity θ in the presence of the other particles specified by filling function n(θ) (see Appendix A for the definition of the dressing procedure and notation).
The diffusion operator is given by the integral kernel where T is the differential scattering kernel (explicitly defined for the Lieb-Liniger model in Function ρ s (θ) specifies the maximal allowed density of particles in the range [θ, θ + ∆θ] and is related to the particles density ρ p (θ) through the filling function n(θ): ρ s (θ) = n(θ)ρ p (θ).
The dependence of all the above functions on space-time variables has been skipped for the compactness of the notation.

Lieb-Liniger model
The Lieb-Liniger model [48,49] describes a gas of bosonic particles in one spatial dimension with ultra-local interactions. The Hamiltonian for N particles is where we set = 2m = 1. We focus on repulsive interactions, c > 0, and assume periodic boundary conditions for a system of length L. The Lieb-Liniger model is exactly solvable. At the center of the exact solution is the scattering kernel The thermodynamics in the limit of finite particles density N/L (which we assume from now on to be 1) is known [45,[50][51][52] and takes a universal form where (θ) solves an Thermodynamic Bethe Ansatz (TBA) equation. At the thermal equilib- where T is the temperature and µ is the chemical potential. The generalized TBA [45] describes in the similar framework the filling functions of non-thermal states. An example of such a state is a state reached after the interaction quench from the BEC state, the ground state of the non-interacting gas c = 0. In that case the distribution of quasi-particles is known explicitly to be [53] n quench (θ) = a(θ/c) a(θ/c) + 1 , a(x) = 2π cx sinh(2πx) with I j (x) the modified Bessel functions of the first kind.
The bare momentum and energy are p(θ) = θ, E(θ) = θ 2 . The physical observables in the Lieb-Liniger model depend on the particles density ρ p (θ) which follows from the filling function n(θ) through the same dressing procedure where we understand that 1/(2π) is a constant function.
In this work, we focus our attention on the filling function n(θ) instead of a particle density ρ p (θ). The reason is two-fold. First, the linearized hydrodynamics is most easily formulated Figure 1: The particle density ρ p (θ) and filling function n(θ) for a thermal gas at T = 1 and with c = 1.
for the n(θ) function. Therefore to understand the resulting dynamics it is most natural to look at n(θ) instead of ρ p (θ). Second, beside the non-linear relation the quantitative features of ρ p (θ) are visible already at the level of n(θ). The main effect in the transformation (13) from n(θ) to ρ p (θ) is in smoothing the shape of the function and, for small values of c, making the distribution more compact, see Fig. 1.
This has two applications. First, when one indeed considers perturbing a space-time independent equilibrium with a small disturbance. 1 The final state of the GHD evolution is then the original equilibrium state.
The second case is in the large-time hydrodynamics of an initially strongly inhomogeneous state. We expect that at late times the GHD dynamics smooth out the quantum fluid to a state which can be described by the linearized equations around a homogeneous state. One must, however, remember that the hydrodynamical modes can only describe systems which are in local thermodynamical equilibrium [26].
The Eq. (2) of the GHD was expressed in terms of the particles' density ρ p (θ, t, x). To consider its linear regime it is convenient to write it as an equation for the filling function and its linearized form is Here we have explicitly written the action of the kernel operator on the filling function. The effective velocities v eff and the diffusion kernelD(θ, α) are determined from the equilibrium data n as in Eq.(14) 2 .
The time evolution of linearized hydrodynamics (16) is driven by two phenomena. First, the velocity with which particles propagate depends on their rapidity. This leads to spatial spreading of an initially localized profile δn(θ, t, x = 0) even in the absence of the diffusion. On top of this the diffusion processes lead to the redistribution of the particle content δn(θ, t, x) through scattering processes with the background, see Fig. 2 The effect of the GHD dynamics is the reorganization of the spatial and rapidity dependence of the filling function. The linear dynamics is conservative in the sense that the total profile of rapidities dx δn(θ, t, x), does not vary in time.

Solution to the linearized GHD
Equation (16) is an integro-differential equation for δn(θ, t, x). The variables are separated and we can solve it in a standard manner. For a single momentum mode where the operator explicitly depends on the Fourier mode k. Hereafter, we shall consider k = 0 modes only: k = 0 does not evolve in time at the linear level and should be part of the background equilibrium state. At the quadratic level δn 0 acquires non trivial dynamics as it shall be discussed in Sec. 5.
The operator D k being an integral, linear operator on non-compact quasi-momenta space has an infinite number of eigenstates f k,ω (θ) and corresponding eigenvalues z k,ω ∈ C, where ω enumerates the eigenvalues, κ ω (k) and λ ω (k) are both real and control the ballistic and diffusive propagation respectively. In writing (21) we have anticipated that, once for each k a suitable ordering of ω's is made, the z k,ω 's are well-defined functions of k. From the numerical diagonalization, we observe that the spectrum of D k is non-degenerate and can be ordered with κ ω . An example of such spectrum is given in Fig. 3. From (20) it follows that κ ω (k) = −κ ω (−k) and λ ω (k) = λ ω (−k). Moreover, we shall see that all λ ω are positive and display the equilibration rate of the mode f k,ω . The slowest equilibration is given by the smallest λ ω .
A general solution to (19) in terms of the eigenstates of D k is where coefficients c k,ω are specified by the initial state The integration is performed over ω parametrizing the spectrum of D k operator. In practice, when we perform numerical computations, operator D k becomes a finite matrix and its spectrum contains finite number of eigenvectors which turns the integration into a sum. The integration measure dω is irrelevant, after discretization it can be absorbed into the coefficients c k,ω . Note that D k operator is not normal and therefore its eigenstates do not form on orthonormal basis. However, they do form a complete basis and (23) has a unique solution.
To quantify the time evolution we focus on the dynamics of momentum modes, and on the moments of the filling function δn(θ, t, x), We will also look at the diffusion coefficient describing the linear in time growth of the spatial variance of the particles' distribution. To this end, we define the diffusion coefficient D j associated to the j-th moment as where Withx j we denote the average position with respect to the j-th moment, defined as The relation (26) can hold only approximately. The time evolution, even in the linearized regime, is quite involved, and can not be simply described by a number of diffusion coefficients.
However, we will see that for intermediate times the relation (26) holds and helps in clarifying the quantum fluid dynamics. At larger times it breaks due to the finite length L of the system.

The diffusion operators
We will now display some properties of the diffusion operatorD (5) and a related operator D k (20). The eigenvalues ofD are non-negative real numbers as follows from its construction [43]. Examples of spectra ofD obtained through numerical diagonalization for thermal and BEC-quench distributions are shown in Fig. 3 . We observe that all eigenvalues ofD but one depend continuously on the eigenvalue index. The special zero eigenvalue is related to the Markov property of the diffusion matrix [43]. The corresponding left eigenfunction is, in the Lieb-Liniger model, the constant function. This can be easily seen from the definition of the diffusion kernel (4).
The Lieb-Liniger model in the c → 0 limit is a theory of free bosons, whereas for c → ∞ maps to free fermions (Tonks-Girardeau gas). In both cases, there is no diffusion because there is no interaction between quasi-particles. This is supported by the large c analysis of the diffusion operator presented in [43]. In the inset of Fig. 3 we plot the dependence of the trace ofD for the BEC-quench saddle point state as a function of the interaction strength c which shows the expected behaviour.
We turn now to the analysis of the spectrum of the D k operator, see also Fig. 3. The diagonal operator with values v eff reigning the ballistic evolution andD do not commute with each other. This rebuilds the spectrum of the resulting operator D k with respect to that ofD.
One of the consequences is that D k has no zero modes. Numerical computations reveal that the eigenvalues are smooth functions of k and with great accuracy follow a scaling with k, as illustrated in Fig. 4. This simple scaling implies that modes with larger k diffuse much faster than modes with smaller k. The eigenvalues κ ω (k) describe the angular frequency of ballistic propagation of a mode, whereas eigenvalues λ ω (k) are responsible for its diffusive decay. The dimensionless ratio of the numbers defines two regimes of the propagation where either ballistic or diffusive propagation dominates. Because of the observed scaling, we can introduce an effective momentum k * Modes with k > k * ω propagate diffusively, modes with k < k * ω propagate ballistically. The effective momentum depends on the eigenstate and changes over few orders of magnitude as also illustrated in Fig. 4. The scaling of the eigenvalues is due to the fact that the diffusion is

Examples
In this section we consider the dynamics governed by the linearized GHD over a thermal state. We have checked that the dynamics over the BEC-quench saddle point state (12) yields qualitatively similar results. We consider two types of the initial state. The first one has a well-defined momentum and therefore its time evolution is especially simple.

Single Fourier mode initial state
We start with the single mode perturbation δn of a well-defined momentum k. The initial perturbation is of the product form 3 where n is the GGE background filling fraction, see Spatial time evolution is displayed in Fig. 6 where we plot the 0-th moment of δn(θ, t, x).
Because of the wave-like form of the initial state the time evolution is periodic in space.
Therefore, we focus on x = 0, which corresponds to one of the maxima of the initial state, and consider local moments of δn(θ, 0, t), that is µ j (0, t). During the time evolution, the local 3 In the following we shall suppress an overall small parameter because it scales out from the linear equation. 4 The Fermi momentum kF = π for the unit density of the background state conglomerates of particles at x multiplicities of k/π get smeared over the whole system. Particles with higher rapidities travel faster and therefore the higher moments of the distribution vanish faster, see also Fig. 6.
Finally, we look into the rapidities distribution at a specific position in space. We choose again x = 0. The evolution of δn(θ, 0, t) is shown in Fig. 7 for linearized GHD with and without the diffusion. We observe that whereas the averaged distribution in both cases vanishes similarly the mechanism is different. The presence of diffusion causes a steady decay of the perturbation whereas purely dispersive evolution leads to a superposition of decoherent for the entropy production (33). The discrepancy between the GHD predictions and observed entropy originates from linearization of the full dynamics.
fluid [43] as we will now present.
In purely non-diffusive GHD, with full non-linear dynamics, the entropy is conserved. Including the diffusion leads to the production of the entropy until the final state of the system settles in, as shown by the solid line in the right panel of fig. 7. When we linearize the GHD dynamics we introduce an error of the second order in the perturbation. This leads to the entropy production already at the non-diffusive level and to entropy production at the diffusive level larger then the one predicted by the full non-linear GHD with diffusion. These expectations are confirmed in the right panel of Fig. 7. The results show that scrambling in the space of rapidities is related to the entropy production.
We can summarize these findings concluding that the evolution of the perturbation is governed mainly by the dispersion suppression mechanisms. Particles with different rapidities propagate with different velocities which quickly (at the order of a single frequency period κ ω (k)) leads to equilibration of the initial spatially delocalized profile. However, looking into the details of local filling functions δn(θ, t, x) reveals the role of the diffusion in the equilibration. Whereas the ballistic propagation leads just to a decoherent superposition of waves, it is the diffusion that redistributes the rapidities and leads to a uniform, for a chosen initial state, long-time distribution.

Localized initial state
We consider now a bit more involved example of a perturbation localized initially in space.
The initial configuration is over the same thermal equilibrium background n(θ). The time evolution and ultimately the equilibration, of the first and second moments µ j (t, x), as defined in (25), is shown in Fig. 9.
The final state of the system is given by the thermal background over which there is a uniform in space and time reminiscent of the initial perturbation given by (see (22)) It is worth pointing out that there are no dissipation processes and that the total density of quasi-particles is conserved dx δn(θ, t, x) = 2π dω c 0,ω f 0,ω (θ).
At large times the density of quasi-particles gets spread uniformly over the space as (35) indicates.
To quantify the process of equilibration we consider the evolution of modes δn k (t) defined in (24). The results are presented in Fig. 10   the decay and that the diffusion slightly speeds up this process. We also estimate the diffusion coefficient for different moments µ j (t, x), see also Fig. 9. The higher moments diffuse faster with the density (0-th mode) decaying the slowest.
Finally, we look at the "occupation numbers" c k,ω (θ, t) quantifying the number of ballistic and diffusive modes. We define the "ballistic occupation" and "diffusive occupation" by with k * ω defined in (30). We expect the diffusive part to decay much faster than the ballistic part. These intuitions are confirmed on Fig. 11 where we show the time evolution of both quantities for the initial perturbation given by the Gaussian packet (34).
Concluding, the time evolution of the localized perturbation is again driven by the ballistic propagation with diffusion playing a secondary role. Moreover, the modes with diffusive dynamics decay faster. Despite this, there is a region in time in which the dynamics is diffusive, in the sense that the spatial variance of expectation values of local conserved charges grows linearly in time. This is caused by an interplay of the ballistic propagation, which mixes quasi-particles in the real space, and diffusive propagation, which mixes modes in the rapidities space. and c diffusive (t) from equations (37) and (38).

Quadratic corrections
In this section we go one step beyond the linear approximation discussed until now and consider quadratic GHD. We focus on the k=0 mode as this mode does not evolve at the linear level. We also dismiss the quadratic part of the diffusion operator as it gives only subleading corrections to k = 0 modes.
The full GHD equation (15) can be easily expanded to the second order in δn(θ, t, x) keeping small diffusion term at linear order only: The leading second order correction comes from the effective velocity term because the latter depends on the full state n(θ, t, x). Directly from the definition of the dressing procedure in Expanding δn and δv in spatial Fourier modes leads to the following equation for the k = 0 Here, δn q is the linear solution to the GHD and determines δv p according to (40). The equation (41) can be now solved numerically with the techniques introduced in the previous sections.
In Fig. 12 we plot the dynamics of the k = 0 mode resulting from (41). In the left panel of Fig. 12 we consider the delocalized initial state from Eq. (31). Such a state has a well-defined momentum k and we consider its three different values. The slower the momentum the slower the resulting evolution of the k = 0 mode to its late-time value.
We consider also the dynamics of the k = 0 mode with the localized initial state (34). To simplify the problem we approximate the full dynamics given by (41) by including only the lowest contributing modes. This provides a lower bound on the equilibration process. In the right panel of Fig. 12, we plot the moments of the resulting distribution of quasi-particles, δn j 0 (t) = dθ θ j δn 0 (θ, t).
The higher moments, which depend more strongly on large θ part of the distribution δn 0 (θ, t) equilibrate faster. This is simply caused by their larger effective velocity.
It is worth to stress that the k = 0 mode tends to a constant non-zero value as times goes to infinity. That means that quadratic fluctuations modify the equilibrium state. Whereas linear dynamics describes decaying fluctuations around the initial profile, inclusion of the quadratic terms leads to the change in the homogeneous profile itself. This is expected from general properties of GHD which has now quite solid ground based on various numerical checks (see e.g. recent works [54,55]). Here we can see it in a much simpler setup.

Conclusions
In this work, we have considered numerical solutions to the Generalized Hydrodynamics in the linear regime. This required turning an initial integro-differential equation into an eigen-problem of the diffusion operator. The latter being exactly known in terms of the quantities of the Thermodynamic Bethe Ansatz. We have performed numerical diagonalization of this operator. Questions about the analytic construction of its eigenstates remain open.
The general solution to the linearized GHD presented in Section 3 was then applied to two concrete examples of initially delocalized and localized states. In Section 4 we focused on the time evolution over the thermal background, however, we have checked that time evolution over the non-equilibrium saddle-point state like the one obtained after the BEC-quench leads to qualitatively similar results. The results show the complementary role of the ballistic and diffusion effect on the dynamics. Whereas the ballistic propagation dominates the spatial homogenization of the initial state, it is the diffusive dynamics that reorganizes the quasiparticle content and ultimately equilibrates distribution also in the rapidities space. It also increases the total entropy of the fluid but the linear approximation applied in this paper overestimates the exact result for the entropy production.
On the other hand, we have observed that the effect of the diffusion is rather weak and potentially can be treated in a perturbative way. A first step in developing such an approach, in a context of the bi-partite quench, was performed in [43]. It would be interesting to pursue this approach further.
Finally, we have also considered the effect of the quadratic terms. These are most important for the evolution of the k = 0 spatial mode which does not evolve in the linear regime.
Quadratic corrections do not show any sign of instabilities which could appear due to positive interference of linear modes. Instead, they tend to non-zero values at late times modifying the homogeneous GGE state. This behavior supports the GHD as a valid approach to non-equilibrium physics of the integrable systems.
In this work we focused on the simplest and most easily accessible observables. It would be interesting to consider the dynamics of other observables, e.g. the local n-body correlations functions for which exact expressions in the homogeneous Lieb-Liniger model were recently obtained [56]. Under the hydrodynamic assumption, these can be extended to the inhomogeneous case and evaluated for the linearized or full GDH dynamics. We leave this problem for future work.