Non-equilibrium steady states in quantum critical systems with Lifshitz scaling

We study out-of-equilibrium energy transport in a quantum critical fluid with Lifshitz scaling symmetry following a local quench between two semi-infinite fluid reservoirs. The late time energy flow is universal and is accommodated via a steady state occupying an expanding central region between outgoing shock and rarefaction waves. We consider the admissibility and entropy conditions for the formation of such a non-equilibrium steady state for a general dynamical critical exponent z in arbitrary dimensions and solve the associated Riemann problem. The Lifshitz fluid with z=2 can be obtained from a Galilean boost invariant field theory and the non-equilibrium steady state is identified as a boosted thermal state. A Lifshitz fluid with generic z is scale invariant but without boost symmetry and in this case the non-equilibrium steady state is genuinely non-thermal.


Introduction
Fluid theory is one of the oldest effective descriptions in physics. 1 It is based on general symmetry principles and applies in the limit of long wavelength and low frequency compared to characteristic microscopic length and time scales of the system in question. A fluid description can thus stand on its own and be useful even when no microscopic description, based on particles or quasiparticles, is available. There has been considerable recent interest in extending fluid theory to systems with unconventional symmetries, including Lifshitz scale symmetry, with potential applications to quantum critical systems [2][3][4][5]. Motivated by these developments, we will consider a problem involving out of equilibrium energy transport in fluids with Lifshitz symmetry.
It remains an open problem to develop a general fluid dynamics formalism for systems that are far from thermal equilibrium, but there has been interesting recent progress in this direction involving relativistic fluids. Investigating out of equilibrium energy transport between two relativistic quantum critical heat baths led to the discovery of the emergence of a universal Non-Equilibrium Steady State (NESS) between the two heat baths, described by a Lorentz boosted thermal state [6][7][8][9][10]. In the present paper, we extend this analysis to more general quantum critical fluids, in particular to non-relativistic fluids with Lifshitz scale symmetry (referred to as Lifshitz fluids in the following), and find that a NESS emerges here as well. For the special case of a Lifshitz fluid with dynamical critical exponent z = 2, the resulting NESS can be viewed as a Galilean boost of a thermal state. For Lifshitz fluids with z = 2, there is no underlying boost symmetry [11]. It turns out there is still an emergent NESS at generic z, but in this case it cannot be obtained as a boosted thermal state.
In order to gain further insight into emergent hydrodynamic behaviour, we adapt the local quench construction of [6,7] to the case of a non-relativistic fluid with Lifshitz scaling symmetry and study the subsequent time evolution for different values of the dynamical critical exponent. We begin in Section 2, where we introduce general properties of such fluids and continue in Section 3 by describing the setup involving a pair of quantum critical heat baths that are brought into contact at t = 0. In Section 4 we briefly review the theory of shock and rarefaction waves that can appear in this context and associated stability conditions. In Sections 5 and 6 we apply the general theory to our specific system, first for the case of a z = 2 scale invariant fluid with Galilean boost invariance and then for a general z = 2 fluid without boost symmetry. Finally, we discuss some open questions and possible future directions in Section 7.

Perfect fluids with Lifshitz symmetry
For simplicity, below we will focus on the special case of perfect fluids. These are idealised fluids, that are without shear, strain or bulk viscosity and do not conduct heat. We begin by introducing the symmetries we will be assuming and the definition of the dynamical critical exponent z.

Symmetries of relativistic and non-relativistic critical fluids
Symmetries play a central role in any fluid description. The most basic symmetries are time translations, spatial translations and spatial rotations, generated by the operators g = {Ĥ,P i ,Ĵ ij }, respectively, whose commutators form the so-called Aristotelian algebra.
A relativistic fluid is not only invariant under these symmetries, but also under Lorentz boostsL i relating observers moving with respect to each other with constant velocity, where v describes the relative velocity between the two observers and γ = 1/ 1 − v 2 /c 2 . At low velocities v c, the Lorentz boost reduces to the Galilean boostĜ i , The Aristotelian algebra is extended to the Poincaré algebra or the Galilei algebra, depending on which of these boost generators is added to g. Furthermore, the Galilei algebra allows for a central extension, known as the Bargmann algebra [12], by the inclusion of an additional symmetry generatorM, such that the non-vanishing Galilean boost commutators are given by (2. 3) The chargeM corresponds to the non-relativistic kinetic mass [13] and needs to be included when describing a fluid with mass density. In a theory with Galilean boost symmetry, the kinetic mass is a measure of the amount of matter in the system and does not vary between inertial frames. It is a conserved quantity in an isolated system. On top of this, in a relativistic critical fluid there is an additional symmetry under dilations of the form Invariance under this symmetry implies that physical processes happen in the same way, at all distance scales or, alternatively, energy scales. For relativistic fluids, the scale symmetry is compatible with Lorentz symmetry and together they place powerful constraints on the allowed dynamics of the fluid. A non-relativistic critical fluid can be scale invariant too, but in this case dilationsD take the more general form of a Lifshitz symmetry, where z ≥ 1 is referred to as the dynamical critical exponent. In the absence of boost symmetries, a closed algebra exists for any z consisting of the generators g z = {Ĥ,P i ,Ĵ ij ,D}.
A key observation, however, is that Lifshitz symmetry with generic z > 1 is in general not compatible with boost symmetry. Indeed, Lorentzian boost symmetry is only compatible with z = 1, which gives the scaling (2.4) and the no-go result of [11] implies that the Galilean boost symmetry (2.2) is only compatible with z = 2 Lifshitz scaling. In the special case of z = 2 the Bargmann algebra can be be further extended to the Schrödinger algebra [14] involving the set {Ĥ,P i ,Ĵ ij ,Ĝ i ,D (z=2) }. For this reason, when discussing the out of equilibrium dynamics of non-relativistic fluids, we will consider separately the cases z = 2 and z = 2, leading to different conclusions about the nature of the emergent steady state.

Thermodynamics and stress-energy tensor
Based on the considerations above, we will consider a fluid whose description is invariant under time and space translations as well as rotations. In addition, we will also assume a global U (1) symmetry whose corresponding conserved charge is N . This is realized by the basic set of generators {Ĥ,P i ,Ĵ ij ,M}. Additional symmetries under boosts and rescaling will be considered below.
Global quantities in this fluid include the energy E, momentum P, entropy S and charge N . Locally, we have the energy density E = E/V , momentum density P i = P i /V , entropy density s = S/V and charge density n = N/V . Assuming a configuration where these can be uniformly defined, the fundamental thermodynamic relations relating the change of the internal energy to the changes in the rest of the thermodynamic state functions are 6) or, in terms of the associated densities, The thermodynamic forces associated to these parameters are the temperature T , the pressure P , the fluid velocity v and the chemical potential µ. As argued in [11], assuming a fluid with uniform velocity v in the presence of rotational symmetry, the momentum density must be proportional to the only directed quantity in the fluid, i.e. the velocity, 8) and the above thermodynamic relation becomes The quantity ρ is referred to as the kinetic mass density. In a theory with Galilean boost symmetry it is proportional to the charge density n but in the absence of boost symmetry the relation between n and ρ is more complicated. The dynamical variables enter into the stress-energy tensor of the fluid T µ ν and the current J µ , whose conservation equations read 2 Classically the symmetry generators are realised by which provides direct interpretation for various components of the stress-energy tensor and current. In particular, the energy density is E = −T 0 0 , the momentum density is P i = T 0 i , and the charge density is n = J 0 in any frame.
For a perfect fluid there exists a reference frame, the rest frame, in which there is no momentum density. The charge current then reduces to just the charge density and the stress-energy tensor involves only two parameters, the energy density and pressure. Explicitly, in this frame we have In any other frame of reference the description will also depend on the velocity v and in the absence of boost symmetry the v dependence can be non-trivial. If the perfect fluid has Lorentz boost symmetry, the stress-energy tensor and current in the moving frame are related to those in the rest frame by a Lorentz boost transformation (2.1). In Section 5, we will be interested in describing a non-relativistic perfect fluid with Galilean boost symmetry under (2.2). In this case the stress-energy tensor and current in the moving frame are obtained from the following transformation rules [15], where we define Γ µ = 1 2 |v| 2 , v . Note that this version of T µ ν does not follow the usual tensor transformation properties, because it does not have tensorial status in the context of Galilean relativity. However, it is possible to combine T µ ν and J µ into an d × (d + 1) dimensional objectT = (T, J) which acts as a tensor. 3 The conservation equations (2.10) are merged into one, and spacetime is embedded into a higher-dimensional construction of Bargmannian coordinates where a tensorial description arises naturally. For an overview of this description in the context of Bargmann theory, see [15] and [16].
Applying (2.13) to a perfect fluid which is flowing at constant velocity v, and described in the rest frame by (2.12), we obtain the following stress-energy tensor and current components [17], (2.14) where E = E 0 + 1 2 n v 2 adds kinetic energy to the energy density. From the off-diagonal components we read off the momentum density P i = T 0 i = n v i , which fixes the coefficient in (2.8) to be ρ = n.
This last observation can also be obtained from the Ward identity corresponding to Galilean boost symmetry. The boost generator can be written asĜ i = t ∂ i = G µ i ∂ µ . Due to the non-vanishing Poisson bracket [P i ,Ĝ j ] in (2.3), the boost current is b µ i = t T µ i − x i J µ and the associated Ward identity gives T 0 i = J i [18], from which ρ = n follows. The physical interpretation is that the flow of matter gives rise to momentum density and the inhomogeneous term in the transformation of the stress-energy in (2.13) accounts for the addition of momentum density under Galilean boosts.
In Section 6, where we consider critical fluids with generic z, we do not assume any boost invariance and the kinetic mass density ρ and the particle number density n are no longer identified with each other. Instead, we adopt an ansatz where they appear separately in the stress-energy tensor and the current [11], and then study out of equilibrium evolution. The Lifshitz scaling relation (2.5) with z = 1 implies that space and time coordinates have different scaling behavior and this affects how dimensional analysis is carried out. The energy E is a conserved quantity associated to time translations, so it must scale as the inverse of time, and thus the energy density scales as E = Λ −(d+z) E. On the other hand, the individual terms in the thermodynamic relation (2.9) must all have the same scaling and from there one can infer the scaling behavior of the various thermodynamic variables of the Lifshitz fluid: Note that it is only for z = 2 that the kinetic mass density scales in the same way as the charge density.
The symmetry under Lifshitz scaling (2.5) leads to the Ward Identity, z T 0 0 + T i i = 0, which in turn implies the equation of state where d is the number of spatial dimensions. For the particular case of z = 2, the equation of state reduces to d P = 2 E − n v 2 and it is easy to see that a Galilean boost of the form (2.13) to the rest frame gives the equation of state for a fluid at rest d P = 2 E. However, as mentioned above, scale invariance with generic dynamical critical exponent z is incompatible with Galilean boost invariance and we will see this explicitly in Section 6 when we study non-equilibrium steady states of a quantum critical fluid with z = 2. In this case, the state variables of a uniformly moving fluid are not equivalent to those of an equilibrium configuration viewed in a moving reference frame.
3 Local quench between semi-infinite heat baths The specific system we consider consists of two semi-infinite heat reservoirs in d spatial dimensions, which are brought into contact at time t = 0 across a flat interface orthogonal to the x-coordinate axis. An equilibrium state of a charged quantum critical fluid is characterized by two energy scales, often taken to be the temperature and the chemical potential (due to scale invariance it is only the ratio T /µ that is physically relevant). In the case at hand, we find it convenient to instead use the pressure P L,R and charge density n L,R of the two reservoirs to describe the initial state, and our solution to the resulting fluid dynamical problem will be expressed in terms of the scale invariant ratios P L /P R and n L /n R . In what follows, we will consider P L /P R > 1 without loss of generality, and arbitrary charge ratio, 0 < n L /n R < ∞.
A local quench of this type, with sharp jump functions θ(x), can serve as a first step towards studying out of equilibrium dynamics in a fluid. The pressure difference between the two reservoirs drives a fluid flow between them. One might intuitively expect the sharp initial gradient to be steadily smoothed out with the system approaching local equilibrium in the central region, but at the level of leading order hydrodynamics this is not the case. Instead, as time evolves, a non-equilibrium steady state (NESS) occupies a growing region between the two heat baths, characterised by the presence of a non-zero, constant energy flow, as was discussed in [6,7]. The properties of the NESS are constrained by the equation of state of the heat baths and the conservation of the stress energy tensor and the charge current across the wavefronts, which emanate from the contact region (see Figure 1).
An initial value problem in hydrodynamics with piecewise constant initial data, where two fluids at equilibrium are joined across a discontinuity, is an example of a so-called Riemann problem [19] in the theory of partial differential equations. A solution, which generically involves shock and rarefaction waves propagating outwards from the initial discontinuity, can be found via the techniques described in Section 4, allowing the fluid variables that characterise the resulting non-equilibrium steady state to be determined in terms of the relevant input data. A Riemann problem for a relativistic quantum critical fluid in general dimensions was studied in [7]. Initially, both outgoing wavefronts were assumed to be shockwaves but it was later realized [8,9] that above two spacetime dimensions, a solution with one shockwave and one rarefaction wave is preferred, based on entropy arguments and backed by numerical analysis. The existence and universality of the steady state for higher dimensional CFTs was studied in [20].

Formulation of the Riemann problem
In the present Riemann problem, the heat reservoirs are brought into contact across a planar surface, that we can take to be orthogonal to the x-axis. Following [6,7], we look for a solution with wave fronts, traveling in the x-direction, that separate space into regions.

1.
A region on the left, with the fluid at rest and stress-energy tensor as in (2.12) with E L , P L and n L .
2. Steady state region (or regions) in the middle, with the fluid flowing at a constant flow velocity v, and stress-energy tensor as in (2.14) with E s , P s and n s .

3.
A region on the right, with the fluid at rest and stress-energy tensor as in (2.12) with E R , P R and n R .
Drawing from the expressions presented in (2.14), in each region the conservation equations (2.10) take the following form: These equations are supplemented with the equation of state (2.17) that relates E and P in a way that reflects the scaling symmetry of the fluid system. Thus, the dynamics is governed by a set of hyperbolic conservation laws of the form where φ and f are functions of the same fluid variables and f (t, x) represents the flux of the conserved quantity φ(t, x). In our non-relativistic quantum critical fluid, the conserved quantities are charge, momentum and energy densities, and the resulting conservation equations (3.2) may be written as Let us now discuss briefly the possible wave solutions that will emerge in this system.

Wave analysis
Generically, let us consider a conservation law of the form mentioned above, for a field φ(t, x), together with a piecewise constant initial condition: This problem was first considered by Riemann in the 19th century [19]. Note that for any given solution of this problem φ sol (t, x), the rescaled function φ θ (t, x) = φ sol (θt, θx) is also a solution for any θ > 0. In fact, the initial condition (4.2) selects, out of all possible solutions of the conservation equations, those which are invariant under such a scaling transformation. These solutions are constant along rays emanating from the origin (t = 0, x = 0) due to the scaling, and they can generically be understood in terms of waves.

Linear problem
In the problem we will be considering, φ is a vector whose components are the energy density, pressure and fluid velocity, but, for the present discussion, we simply take it to be a generic vector of k components. A simple special case is obtained when ∂ that is, for the strictly hyperbolic system where A is a matrix of constant coefficients. In this case, any solution can be written as a superposition of traveling waves. A generic initial condition φ(0, x) =φ(x) defines a wave profile that is shifted to the left and right as it evolves in time, in such a way that the height of the evolved profile at a given point is the sum (superposition) of heights at different points of the original profile. The explicit solution takes the form where λ i are the eigenvalues of the matrix A, that determine the speed of propagation of each component of φ, while the coefficients of the superposition, r i , are the components of the corresponding eigenvectors of A, and they determine the direction of the rays along which the wave travels. By diagonalising the matrix, the problem is decomposed into k scalar Cauchy problems that can be solved separately. Figure 1: Propagation of shock, contact discontinuity and rarefaction waves for P L > P R .

Non-linear problem
More generally, the Jacobian in (4.1) is a function of φ itself, This adds non-linearity to the problem. The solution can still be written in terms of waves, but the waves can interact with each other, producing additional waves. This is because the eigenvectors r i are generalised into functions which depend on φ. The eigenvalues λ i also depend on φ, and so the shape of the various components of the solution will vary in time, leading to wave dispersion and compression.
In [21], Lax provided a classification of the waves that can arise in non-linear wave problems with initial conditions of the form (4.2). To do so, he introduced a simplifying assumption: that each λ i (φ), that is, the i th eigenvalue of the Jacobian matrix (4.5), corresponds to either a genuinely non-linear wave, such that ∇λ i (φ) · r i (φ) = 0 for all φ, or to a linearly degenerate wave, such that ∇λ i (φ) · r i (φ) = 0 for all φ. The quantity ∇λ i · r i can be understood as the directional derivative of λ i (φ) in the direction of the vector r i .
As we will see below, this assumption holds in our Riemann problem for Lifshitz fluids and the resulting solutions have a simple structure consisting of different kinds of waves or discontinuities, which can be classified as follows: • The linearly degenerate case ∇λ i · r i = 0, for which λ i is constant along each integral curve of the corresponding field of eigenvectors r i . In this case the profile of the solution does not change in time, generating a so-called contact discontinuity.
• The genuinely non-linear case with ∇λ i · r i > 0 such that the i th eigenvalue λ i is strictly increasing along the integral curve of the corresponding field of eigenvectors r i . This leads to a rarefaction wave, displaying a smooth profile that widens and decays over time.
• The genuinely non-linear case with ∇λ i · r i < 0. This leads to a shock wave, displaying a compression which makes it become steeper over time.
When the simplifying assumption described above is valid, a set of stability conditions can be formulated which guarantee uniqueness and a continuous dependence on the initial data [22]. The one relevant for our analysis is Lax's shock wave admissibility condition [23], which can be easily visualised for the Riemann problem, where the initial configuration of φ(0, x) jumps from a left state φ L to a right state φ R at some value of x. The information contained in the piecewise initial condition propagates forward at speeds given by λ i (φ L ) on the left and λ i (φ R ) on the right. In order to prevent new characteristics spawning away from the shock interface, which would amount to non-uniqueness for our Cauchy problem, one must impose λ i (φ L ) ≥ λ i (φ R ). Furthermore, a shock wave connecting the states φ L , φ R moving at speed λ = u s , must satisfy Lax's admissibility condition applies to shock waves but not to rarefaction waves. For a rarefaction wave, the solution's admissibility is determined by requiring λ i (φ) to increase smoothly along the profile.
5 Rarefaction and shock waves for a z = 2 Lifshitz fluid As already mentioned in Section 2, a non-relativistic Lifshitz fluid with scaling exponent z = 2 is special. This is due to a number of reasons. First of all, a Galilean boost invariant field theory describing such a fluid has been explicitly constructed [3,24]. In addition, for z = 2, the Schrödinger group (consisting of the Bargmann group, enhanced by the addition of the dilation operatorD), can have an additional generator,Ĉ, corresponding to special conformal transformations. Finally, as shown in [11] and [25], it is only for this particular value of z that one can have a Galilean boost invariant fluid with Lifshitz scaling symmetry with a discrete Hamiltonian and number operator spectrum.
In view of this, we first consider a z = 2 Lifshitz fluid in d spatial dimensions taken to be invariant under Galilean boosts in addition to the scaling symmetry. In this case, we have the relation ρ = n by virtue of a Ward identity, so the momentum density (2.8) is Then, the conservation equations (3.4) become where the combination has been introduced and the right hand side has been expressed solely as a function of the conserved variables. This has the form of a Riemann problem (4.1) with φ = (E, q, n). The flux vector f (φ) can be read off from the right hand side and the Jacobian matrix is easily evaluated, One of the eigenvalues of the Jacobian, along with the corresponding eigenvector, is This is linearly degenerate, ∇λ 1 (φ) · r 1 (φ) = 0, and corresponds to a contact discontinuity. The remaining eigenvalues and eigenvectors are and It follows that the λ 2 (λ 3 ) eigenvalue corresponds to a left-moving (right-moving) wave. Lax's admissibility condition for a shock wave turns out to be satisfied if and only if the pressure in the region behind the wave front exceeds the pressure in the region ahead of it. In our problem, where we assume that P L > P R , this is the i = 3 right-moving wave. The left-moving i = 2 wave, on the other hand, advances into a region of higher pressure and is therefore a rarefaction wave, whose profile widens over time. 4 Figure 2 shows a snapshot of the wave profile for a particular choice of initial data in (3.1), with a rarefaction wave on the left, a shock wave on the right, and a contact discontinuity in between. The shape is similar to the solution of the corresponding Riemann problem for a relativistic critical n  n R P  P R Figure 2: Snapshot of wave profiles for the pressure and charge density at t = t 0 > 0 for P L > P R and n L = n R . The NESS region, bordered by the left-moving rarefaction wave and the right-moving shock wave, contains a contact discontinuity in the charge density.
fluid considered in [8,9]. In particular, as we'll see below, the pressure remains constant across the contact discontinuity in the NESS region while the charge density jumps. In the relativistic case, the charge density decouples from the equations that determine the pressure but this is not the case here. For a non-relativistic Lifshitz fluid, the pressure still remains constant across the contact discontinuity but its value in the NESS region is nevertheless influenced by the initial values for the charge density of the two reservoirs (see e.g. (5.33) below).

Rarefaction wave profile
Let us start by analysing the i = 2 rarefaction wave. For this it is convenient to introduce the concept of Riemann invariants. A function R (i) (φ) that is constant along the integral curves of the eigenvector r i , is called an i-Riemann invariant. A system with k eigenvalues has k−1 linearly independent i-Riemann invariants and they provide a convenient way to construct elementary wave solutions that are the building blocks of a full solution to the Riemann problem [19]. In the case at hand, we have two independent Riemann invariants per family of solutions, satisfying ∂R For the first family, λ 1 = q n = v is itself a Riemann invariant, which means that the speed of the fluid is the same on both sides of the contact discontinuity, and additionally that the discontinuity itself moves at the same constant speed. In fact, this wave is called a contact discontinuity precisely because it moves at the fluid flow speed. A second Riemann invariant for the first family is given by the pressure, P = 2E d − q 2 d n , so this quantity remains constant across the discontinuity as well.
For the two genuinely non-linear families, we find the following pairs of Riemann invariants: where γ ≡ d+2 d . In order to facilitate their interpretation, these expressions can be rewritten using the equation of state, where c was defined in (5.9) and we have dropped a multiplicative constant from R 1 and R 1 . We note that c and γ are, respectively, the speed of sound and the ratio of specific heats at fixed pressure and volume in an ideal gas of z = 2 Lifshitz particles in d spatial dimensions [11].
The first Riemann invariant is the same for both the second and third families and involves a combination of pressure and particle density, P n −γ , which remains constant during an isentropic process in an ideal gas. In other words, the conservation of R (i) 1 amounts to the conservation of specific entropy, i.e. the entropy per particle, along integral curves of r i . To see this, write the first law of thermodynamics in the form T ds = de − P n 2 dn, (5.14) where s and e are, respectively, the specific entropy and specific internal energy. When expressed in terms of the specific internal energy, the equation of state (2.17) becomes which reduces to d P = 2ne for z = 2. This implies Inserting (5.16) into (5.14) and applying the ideal gas law, one obtains Thus, the first Riemann invariant in (5.12) may be interpreted in terms of entropy and we note that the second one has the expected form of a Riemann invariant obtained for a compressible Eulerian fluid [26].
For the i = 2 rarefaction wave, the conservation equations (3.4) are solved implicitly by the requirement that both Riemann invariants remain constant along the wave profile, 2 (E, q, n) = R 2 (E L , q L , n L ) .

(5.18)
The left reservoir values E L , q L , n L are realised at the leading edge of the rarefaction wave profile and can therefore be taken as a reference. The above requirement translates into the following two relations: where v s1 denotes the fluid flow velocity to the right of the rarefaction wave (see Figure 1) and v L is the fluid flow velocity in the heat bath on the left (v L = 0 in a heat bath at rest). Equivalently, the first relation in (5.19) can be used to express the flow velocity in terms of pressure rather than charge density, (5.20) The phase velocity of the wave is given by the eigenvalue λ 2 , as seen in (4.4), which in the present case is given by λ 2 = v − c (with c > 0). Taking the wave profile to be parametrised by n, the condition for a valid rarefaction wave solution is On the curve we have and the rarefaction condition holds provided the charge density is higher in the region ahead of the wave front than behind the wave. This is indeed the case when P L > P R . Note that since the wave has a smooth profile with spatial dependence n(x, t), the phase velocity of the rarefaction wave also acquires a profile, λ 2 (x, t). On the leading left wavefront, where n = n L , it evaluates to λ 2 = −c L , that is, to the speed of sound in the heat bath on the left. Similar considerations apply when P L < P R , except in this case the rarefaction wave belongs to the i = 3 family and moves to the right.

Jump conditions and shock wave
Riemann invariants are useful when the wave profile is smooth but other methods are needed for dealing with the sharp transitions that occur across a shock wave. A solution can be found by imposing so-called Rankine-Hugoniot jump conditions [23,27], which express the conservation laws across the wavefront and relate variables in adjacent regions. For the problem (3.3), the jump conditions can be stated as where u s is the speed of the wave front in question and [x] denotes the change in the variable x across the wave front, as described above. Writing w = v − u s and ν = n w, these conditions can be expressed as where we have used the equation of state (5.2) and the definition c 2 = γ P n . A trivial and immediate solution is ν = [P ] = 0, which is the contact discontinuity described by the linearly degenerate i = 1 family of the previous subsection. As discussed below (5.11), the pressure and fluid speed are the same on both sides of the contact discontinuity, P s1 = P s2 ≡ P s and v s1 = v s2 ≡ v s , but in general the energy and particle densities will be discontinuous across the wave front.
A right-moving i = 3 wave presents a non-trivial solution to the jump conditions. Assuming that ν = 0, we introduce dimensionless variables: where the right-most equality follows from the first jump condition (5.25). The remaining jump conditions (5.26) and (5.27) can be re-expressed as respectively. Combining these conditions and solving for y or Π s gives Substituting y back into (5.29), and choosing the branch of the square root that corresponds to a wave moving to the right, leads to the following expression for the shock speed, Here v R is the fluid speed in the heat bath on the right (v R = 0 for a heat bath at rest). With this choice of sign, Lax's admissibility conditions (4.6) are satisfied for the shock wave. Indeed, with λ 3,R = v R + c R = c R , the requirement is u s > c R , i.e. that the speed of the wave front exceeds the speed of sound in the medium that the shock wave expands into. This, in turn, amounts to the condition P s > P R .
Finally, we can use the relation y = w R /w s from (5.28) to obtain the fluid speed v s2 in the region between the shock wave and the contact discontinuity in Figure 1,

NESS variables and Galilean boost symmetry
Earlier we observed that pressure and fluid flow speed are the same on both sides of a contact discontinuity and the discontinuity itself propagates at the same speed. Demanding equality of the expressions for v s1 in (5.20) and v s2 in (5.32) gives us the following scale invariant condition on the pressure in the NESS region between the rarefaction and shock waves, The initial data of the two reservoirs enters through the ratios Π L = P L /P R and η = n L /n R . The above condition is non-linear but can be solved numerically and one finds a unique value of Π s for given Π L and η. The full solution to the Riemann problem can then be mapped out by evaluating the following expressions for the remaining NESS variables in  Left panel: Steady state pressure P s , charge densities n s1,s2 , and energy densities E s1,s2 . Right panel: Flow speed v s , shock speed u R , and wave speed u L across rarefaction profile. Left panel: Steady state pressure P s , charge densities n s1,s2 , and energy densities E s1,s2 . Right panel: Flow speed v s , shock speed u R , and wave speed u L across rarefaction profile.
terms of the pressure, and evaluating (5.31) for the speed of the right-moving shock wave. The speed of the fluid flow in the NESS region can be obtained by evaluating either (5.20) or (5.32). Solutions for d = 3 spatial dimensions are presented in Figure 3 as a function of P R /P L for fixed n L /n R and in Figure 4 as a function of n L /n R for fixed P R /P L . In the solution of the corresponding Riemann problem for a relativistic quantum critical fluid [8,9] the NESS was described by a Lorentz boosted thermal state with a contact discontinuity in the charge density in the fluid rest frame. The behaviour of a z = 2 non-relativistic critical fluid is analogous, although in this case the boost symmetry is Galilean rather than Lorentzian. The fluid variables in the NESS region of the z = 2 flow have a stress-energy tensor and current of the form (2.14). The pressure and fluid speed are the same on both sides of the contact discontinuity but the energy density and the charge density take different values on the two sides. Nevertheless, if we perform a Galilean boost with velocity −v s to the NESS rest frame following the rule (2.13), we obtain a stress-energy tensor of the form (2.12) with P = P s and a uniform energy density, Furthermore, the fluid variables in the rest-frame satisfy the equation of state of z = 2 fluid at rest, E 0 = d 2 P . Since n does not transform under a Galilean boost, there is still a contact discontinuity in the charge density. Indeed, in the NESS rest frame the two fluids are at rest in hydrostatic equilibrium but the charge density is discontinuous across the contact surface. The charge density remains unchanged with time as there is no fluid flow across the boundary and therefore no charge transport. This kind of a sharp charge discontinuity is allowed when we restrict ourselves to leading-order hydrodynamics but is presumably smoothed out by higher-order corrections, which we do not consider here. We note that analogous behaviour was seen in the NESS rest frame of a relativistic fluid in [8,9].
6 Rarefaction and shock waves at general z In this section we turn our attention to a perfect Lifshitz fluid with a general dynamical critical exponent z > 1. This is motivated by the existence of quantum critical condensed matter systems with a general dynamical critical exponent z = 2, such as the heavy fermion metals discussed in [28] and [29]. For generic values of z such a system is without boost symmetry and it is interesting to see how this affects the solution to the fluid Riemann problem that we have been considering. The first thing to note is that the kinetic mass density ρ can no longer be proportional to the charge density n when z = 2. If we assume that ρ can still be expressed as a function of n alone, then the scaling relations (2.16) imply a relationship of the form ρ = m n α , (6.1) and m a constant of proportionality. In principle, one could allow for more general behaviour, for instance by letting ρ depend explicitly on the velocity v as well as on the charge density, but we will not pursue this here. A scaling ansatz of the form (6.1) provides an example of a Lifshitz fluid without boost symmetry and this is sufficient for our present purposes. In what follows, we will take m = 1 for simplicity.
With the above ansatz the thermodynamic relation (2.9) takes the form The dv 2 terms can be absorbed by defining an internal energy and a shifted chemical and then the familiar form of the first law of thermodynamics is recovered, dÊ = T ds +μ dn .
with q = n α v. The analysis of the Riemann problem proceeds along the same lines as before. The equations are more involved when z = 2, and we have to rely on numerical evaluation to a greater extent, but the NESS variables can still be solved for. The Jacobian matrix, df (φ) for general z is and its eigenvalues and eigenvectors can readily be evaluated. They correspond to a linearly degenerate wave, which is a contact discontinuity, together with two genuinely non-linear waves, and where we've introduced the shorthand notation, It is easily checked that the corresponding expressions in Section 5 are recovered when we insert z = 2 in (6.6) -(6.11). Furthermore, by using the equation of state (6.5) one obtains The eigenvalues corresponding to genuinely non-linear waves can then be written, As before, we find that λ 2 (λ 3 ) corresponds to a left-moving (right-moving) wave, and that the leading wavefront of a rarefaction wave will advance at the speed of sound in a heat bath at rest.

Rarefaction wave profile
Now consider initial data of form (3.1) for a Lifshitz fluid with general z and assume that P L > P R . In parallel with the z = 2 case considered in Section 5, this results in a left-moving rarefaction wave, a right-moving shock wave, and a central NESS region with constant flow velocity and a contact discontinuity moving with the fluid. The key difference compared to the z = 2 case is that now there is no boost symmetry and the steady state flow in the central region will no longer be a boosted thermal state. We use Riemann invariants to analyse the i = 1 contact discontinuity and the i = 2 rarefaction wave. The Riemann invariants for the first family of wave solutions are again given by the pressure P and the velocity v, which coincides with the eigenvalue λ 1 = q n α . Therefore, the contact discontinuity will still propagate at the same speed as the velocity of its surrounding fluid regions on the left and right.
For the genuinely non-linear families, we find generalisations of the pairs of Riemann invariants, which took the form (5.13) for z = 2, but are now given by , and K may be read off from (6.12). As before, we require that both Riemann invariants are constant along the characteristic curves of the left-moving rarefaction wave. From R with K expressed as a function of P , v, and n through the relations in (6.12). These conditions are non-linear and do not allow for analytic solution for generic values of d and z. In order to facilitate their numerical solution, we find it convenient to first eliminate the pressure between them by inserting (6.15) into (6.16). This results, after some algebraic manipulations, in the following equation, relating the scale invariant variables n/n L and v/c L , . A numerical solution for v/c L in terms of n/n L can then be inserted into (6.15) to determine P/P L . In order to check the validity of the rarefaction wave solution so obtained, we have evaluated the characteristic speed λ 2 along the integral curve for specific initial data. Numerical results for several different values of z are shown in Figure 5 and in each case the rarefaction condition (5.21) is indeed satisfied.

Shock wave
A shock wave solution for a Lifshitz fluid at general z satisfies the following Rankine-Hugoniot jump conditions, where we have used the equation of state (6.5). The contact discontinuity corresponds to the trivial solution w = [P ] = 0.
To find a right-moving shock wave solution corresponding to the i = 3 eigenvalue family, we again introduce dimensionless variables, where the right-most equality follows from the jump condition [n w] = 0. For a shock wave propagating into a fluid at rest, the other two jump conditions can be re-expressed as .
The two equations can now be combined and solved either for y or Π s , By substituting Π s into (6.21), the speed of the shock wave can be written in terms of the dimensionless variable y as, The shock wave admissibility conditions are satisfied when the shock front moves faster than the speed of sound in the medium the wave is expanding into, i.e. when u s > c R . It is easily checked that this holds for all values of y that correspond to P s > P R . The fluid velocity in the region between the shock wave and the contact discontinuity can also be expressed in terms of y via the relation, v s2 = (y − 1) y u s . (6.24)

NESS variables
We now have everything in place to construct the full solution to our Riemann problem for a Lifshitz fluid with general z in d spatial dimensions, with initial data given by P L , P R , n L , and n R (with P L ≥ P R ). Once again, there will be a growing NESS region between a leftmoving rarefaction wave and a the right-moving shock wave, with a contact discontinuity in between, as depicted in Fig. 1. The solution can be constructed in a number of ways but the key observation is that pressure and fluid flow speed remain constant across the entire NESS region, while the charge density is piecewise constant and makes a jump at the contact discontinuity. We will proceed by first solving for the pressure and flow speed in terms of the charge density on either side of the contact discontinuity. We then require that the results are the same on both sides and this, in turn, fixes the charge densities in terms of the initial data. On the one hand, the NESS pressure and flow speed are expressed in terms of the dimensionless variable y = n s2 /n R in (6.22) and (6.24), respectively. These relations follow directly from the Rankine-Hugoniot jump conditions across the shock wave front.
On the other hand, we can obtain the same quantities in terms of another dimensionless variable x = n s1 /n L by considering the trailing end of the rarefaction wave profile, where n = n s1 . In this case, (6.17) reduces to . This can be solved numerically for v s as a function of x and the result is then inserted into (6.15) to obtain the NESS pressure, The requirement that v s and P s take the same values on both sides of the contact discontinuity gives rise to two independent relations between the variables x and y, which is sufficient to determine their values for given initial data for the reservoirs. 5 The remaining NESS variables are easily obtained once the dimensionless charge densities x and y have been solved for numerically. For instance, the NESS pressure is obtained by inserting y into the equation on the right in (6.22), while the shock wave speed and the fluid speed in the NESS region are given by (6.23) and (6.24), respectively. Solutions for d = 3 spatial dimensions and z = 3 are presented in Figure 6 as a function of P R /P L . Figure 7 shows how the solution changes with z for a particular choice of P L /P R and n L /n R .
As stated above, the NESS for z = 2 cannot be recognized as a boosted thermal fluid. The equation of state (6.5) is incompatible with the Galilean boost transformations (2.14), 5 As in the z = 2 case, the initial data only enters through the ratios PL/PR and nL/nR.  which leave invariant P and n while shifting E → E + 1 2 n v 2 . Furthermore, the momentum density P = ρ v does not match the one obtained from a Galilean boost. Therefore no temperature can be associated to the solution obtained here. It is genuinely a non-thermal out-of-equilibrium state in a theory without boost symmetry.
It is interesting to compare the NESS variables we find at z = 1 to the solution of the corresponding Riemann problem for a relativistic fluid presented in [8,9] in the limit of low flow velocity. The steady state flow is slow when P R /P L is close to 1, i.e. when the pressure difference between the two reservoirs is small. Figure 8 shows the NESS variables n s1 , n s2 , Π s , and v s at different values of P R /P L for d = 3, z = 1, and n L = 2n R . The corresponding variables in a relativistic fluid (taken from [8]) are indicated by red dashed curves in the figure. We see a close match for all the NESS variables as P R /P L → 1.
In the relativistic case, the charge density decouples from the equations that determine the steady state pressure and flow speed but in general this is not the case for our nonrelativistic Lifshitz fluids. The decoupling of the charge density is, however, recovered in the limit of small pressure difference in the z = 1 Lifshitz case. To see this, one carries out an expansion in powers of small ∆ = Π L − 1 in (6.25) and (6.26) that determine Π s and v s at z = 1 and observes that η = n L /n R indeed decouples from the equations to leading order in ∆. For large values of ∆ the steady state flow speed is no longer small and there is no reason to expect a match between a relativistic fluid and a z = 1 Lifshitz fluid.

Discussion
The above study of the Riemann problem for Lifshitz fluids had a twofold purpose. On the one hand, it extends to a non-relativistic setting some recent work on the out-of-equilibrium flow of relativistic quantum critical fluids [6][7][8][9][10], and, on the other hand, it provides an  Figure 8: NESS variables for a Lifshitz fluid at d = 3, z = 1 and n L /n R = 2 as a function of P R /P L . For comparison, the corresponding variables for a relativistic fluid considered in [8,9] are shown by the red dashed curves. The solutions are well matched as P R /P L → 1.
application to a concrete physical setup of a recently developed general formalism for perfect fluids without boost symmetry [11].
We have established that a non-equilibrium steady state, of the type seen previously in a relativistic scale invariant fluid, will also develop in a non-relativistic critical fluid when two reservoirs are brought into contact across a hypersurface. Consistent with the Lax entropy conditions, the non-relativistic NESS is bounded on one side by an outgoing shock wave and on the other side by a rarefaction wave propagating in the opposite direction. Inside the NESS there is a contact discontinuity where the charge density jumps but the pressure stays unchanged.
In the special case of a z = 2 Lifshitz fluid the NESS is a Galilean boost of a thermal equilibrium state, in direct analogy with the Lorentz boosted thermal state seen in the corresponding relativistic problem. Using a simple scaling ansatz for the kinetic mass density of a Lifshitz fluid at generic z, we found that the fluid variables in the central region can be solved for and a NESS forms in this case as well, but the solution is genuinely non-thermal.
There are several future directions to be explored. In this study, we have concentrated on perfect fluids without impurities or lattice effects which break translational invariance. Proceeding along the lines of [8], where this has been done for a conformal fluid, one could allow for diffusion and momentum relaxation in the hydrodynamics equations, to obtain the time scale up to which the non-relativistic NESS persists.
Another interesting direction is to analyse a dual gravitational description of nonequilibrium steady states of Lifshitz fluids. In this context, it would especially be interesting to identify a gravitational dual of a z = 2 Lifshitz fluid flow without boost symmetry.