Hydrodynamics and instabilities of relativistic superfluids at finite superflow

We study the linear response of relativistic superfluids with a non-zero superfluid velocity. For sufficiently large superflow, an instability develops via the crossing of a pole of the retarded Green's functions to the upper half complex frequency plane. We show that this is caused by a local thermodynamic instability, i.e. when an eigenvalue of the static susceptibility matrix (the second derivatives of the free energy) diverges and changes sign. The onset of the instability occurs when $\partial_{\zeta}(n_s\zeta)=0$, with $\zeta$ the norm of the superfluid velocity and $n_s$ the superfluid density. The Landau instability for non-relativistic superfluids such as Helium 4 also coincides with the non-relativistic version of this criterion. We then turn to gauge/gravity duality and show that this thermodynamic instability criterion applies equally well to strongly-coupled superfluids. In passing, we compute holographically a number of transport coefficients parametrizing deviations out-of-equilibrium in the hydrodynamic regime and demonstrate that the gapless quasinormal modes of the dual planar black hole match those predicted by superfluid hydrodynamics.


Introduction
Superfluids are quantum phases of matter whose most remarkable experimental manifestation is frictionless flow, [1][2][3].They arise when a global, continuous U(1) symmetry is spontaneously broken: in other words, an operator charged under the symmetry acquires a nonzero vacuum expectation value.The phase of this complex order parameter is related to the appearance of a new gapless mode in the spectrum, the Goldstone boson φ, directly responsible for frictionless flow.
The Goldstone mode transforms nonlinearly under U(1) gauge transformations, and so gauge invariance demands that the theory only depends on gradients of the Goldstone ∇ µ φ.This makes this degree of freedom topological: its winding can only be relaxed by topological defects, the vortices.A modern understanding of superfluidity is in terms of an anomalous emergent higher-form symmetry, reflecting the conservation of winding of the Goldstone field, [4,5].The mixed 't Hooft anomaly with the original U(1) current sources a nonzero overlap of the charge current with a conserved quantity (the higher-form current) and in turn causes the conductivity to be infinite in the zero frequency limit, for much the same reason it is infinite in a phase with unbroken spatial translations.This is the origin of frictionless flow.
A finite density of the higher-form current corresponds to a non-trivial background winding and a nonzero supercurrent.As it turns out, increasing the background winding eventually triggers an instability and superfluidity is destroyed.Landau gave a famous argument to determine the onset of the instability in Helium 4, [2].Besides the superfluid component, at nonzero temperature the current is also transported by a normal component made up of thermally excited 'phonons', in other words quasiparticles with bosonic statistics and a dispersion relation ϵ k linear at small wavenumbers k in the superfluid rest frame.Assuming that the system enjoys Galilean invariance and boosting back to the normal fluid rest frame, the phonon energy becomes ϵ ′ k = ϵ k + k • v s , where v s is the superfluid velocity transporting superflow.Superfluidity is lost when it becomes energetically favorable to excite phonons, i.e. when ϵ ′ k < 0. The critical velocity at which ϵ ′ k becomes zero depends on the angle between the vectors v s and k and is lowest when they are antiparallel.This gives Landau's criterion for the critical velocity In a previous work [6], a subset of us showed that the instability at large superflow can be recast as a local thermodynamic instability, expressed simply as where p is the pressure (equal to minus the free energy density), n s the superfluid density and v s the norm of the superfluid velocity, as we will define more precisely in Section 2. Computing the pressure for a non-relativistic system of bosonic thermal excitations with dispersion relation ϵ k allows one to reproduce Landau's criterion (1.1).
The advantage of (1.2) is that it does not assume any particular boost symmetry,1 or that the ground state is weakly-coupled.It also gives a formulation independent of the microscopic details of the system under consideration.
Gauge/gravity duality [8,9] allows one to construct strongly-coupled relativistic superfluids, [10][11][12].Holographic superfluids at finite superflow were studied in [13][14][15][16][17][18][19][20][21].[20,21] observed that for a critical value of the superfluid velocity, a pole of the retarded Green's function crosses to the upper half plane, leading to perturbations exponentially growing with time.Here, we show that the critical velocity when this happens is precisely predicted by (1.2).This paper is structured as follows.In Section 2, we review relativistic superfluid hydrodynamics including first-order derivative corrections (earlier work on relativistic superfluid hydrodynamics includes [17,[22][23][24][25][26][27]).In Section 3, working in the limit when the wavenumber and the superfluid velocity are collinear for simplicity, we describe the spectrum of linear, collective excitations and demonstrate (1.2).In Section 4, we construct a holographic superfluid phase at finite superflow, and check that it becomes dynamically unstable when (1.2) holds.In passing, we compute all the thermodynamics of the state as well as the collinear transport coefficients, as well as check that the hydrodynamic dispersion relation for the gapless modes matches the location of the gapless quasinormal modes of the dual planar black hole, as expected. 2We then conclude with some outlook.In Appendix A, we give details on the numerical calculations, while in Appendix B, we analyze the so-called probe limit where temperature and normal velocity fluctuations are frozen.

Review of superfluid hydrodynamics
We follow the derivation in [29] (see also [30,31]), although we will choose a different frame which ultimately proves more convenient for our purposes.
We consider a (d + 1)-dimensional manifold equipped with a metric g µν and a U(1) gauge field A µ , and impose that the theory is invariant under diffeomorphisms χ µ and U(1) gauge transformations Λ.Under these, the background fields transform as (2.1) where the set of gauge parameters is X = {χ µ , Λ χ } and L χ represents the Lie derivative along χ.
In thermal equilibrium, we can use the existence of a time-like Killing vector field β µ to define our notion of temperature and fluid velocity which transform appropriately under diffeomorphisms as a scalar and a vector, respectively.T 0 is a normalization constant which does not play any role, and with our definition, the fluid velocity is timelike and u µ u µ = −1.We can also define a gauge and diffeomorphism invariant chemical potential as The set of gauge parameters X is of course arbitrary, but it is advantageous to choose them to be B = {β µ , Λ B }.Then, the variations of the background fields under B are ) and vanish in thermal equilibrium.
In a superfluid phase, we also define the gauge-covariant superfluid velocity ξ ν , which is built out of the Goldstone field and the external gauge field: together with ζ µ , the component of the superfluid velocity normal to the fluid velocity and the associated projector ∆ µν = g µν + u µ u ν .
The stress tensor T µν and the current J µ of the theory can be computed by varying the background fields g µν and A µ and computing the corresponding variation of the Lagrangian density: Here we have added a variation of φ, as φ transforms under gauge transformation as Thus in superfluid phases, there is another 'current' K, which as we will see obtains a constitutive relation similarly to T µν and J µ .In thermal equilibrium, Λ B = µ/T − β µ A µ , so that we recover the familiar Josephson relation for the Goldstone field Requiring that the variation (2.6) does not depend on gauge parameters yields the conservation equations of the theory in the presence of background sources With these equations in place, we now turn to the constitutive relations.They can be determined by requiring that the divergence of the entropy current is nonnegative: where we have decomposed the entropy current into an ideal part where the normal fluid transports the entropy density s, and a non-ideal part sµ .This decomposition is based on the physical assumption that the superfluid component does not transport entropy at ideal level, but turns out to be correct once we impose that the ideal terms in the constitutive relations produce no entropy.
The pressure depends on all possible diffeomorphism and gauge invariant scalars, p(T, µ, ζ 2 ), where we have chosen to work with ζ 2 rather than ξ 2 , which allows us to work in a convenient frame once we include derivative corrections.The pressure does not depend explicitly on the fluid velocity, as a state with a nonzero background fluid velocity can always be generated by a Lorentz boost from the state where the normal fluid is at rest.However, the norm ζ ≡ ζ 2 of the superfluid velocity relative to the normal fluid velocity is a physical parameter.The variation of the pressure is where s and n correctly turn out to be the entropy and charge densities as expected, while n s will be interpreted as the superfluid density, the fraction of the charge carried by the superfluid component.n, s and n s all depend on T , µ, ζ 2 .
To determine the constitutive relations, it is more convenient to compute which can be shown to be non-negative off-shell as well [29], i.e. without needing to impose the equations of motion (2.9).After some algebra, and with a little bit of hindsight, the constitutive relations are found to be3 (2.17) In passing we have also defined the energy density Restricting for a moment to ordinary fluids, (2.17) shows that no entropy is produced when T µν , Jµ are set to zero, i.e. when the stress tensor and charge current take their ideal forms.The structure of (2.17) and the requirement that ∆ ≥ 0 implies that derivative corrections to T µν , J µ should either drop out from ∆, in which case they are called 'non-dissipative', or contribute positively to ∆, in which case they are deemed 'dissipative', [30,31].Further, these derivative corrections can be separated in those that are terms proportional to (derivatives of) δ B g µν and δ B A µ , which vanish in thermal equilibrium and are then called 'non-hydrostatic', and terms which are not proportional to any of the δ B (•) and so do not vanish in thermal equilibrium, which are called 'hydrostatic'.The hydrostatic corrections capture deviations of the equilibrium state from homogeneity, and can be derived from writing derivative corrections to the pressure in the generating functional, [34][35][36].By varying the external sources, one can then determine how they contribute to the constitutive relations of currents.Instead, the non-hydrostatic terms (either dissipative or not), are true non-equilibrium corrections.
There is a subtlety with this derivative expansion related to the presence of a superfluid.Indeed, the constitutive relation for K involves derivative mixing, coming from our choice of φ as a dynamical field, [29,37].This choice forces us to take φ ∼ O(∇ −1 ) so that gradients of φ contribute to the first law or constitutive relations at order O(∇ 0 ).Then the first term in (2.15) is actually an O(∇) correction, and positivity of entropy production implies that This type of derivative mixing is typical of theories with spontaneous symmetry breaking and can be avoided by trading φ for a higher-form current and the associated anomalous conservation law stemming from the conservation of winding in the absence of topological defects, [4,38].Then all currents are treated on the same footing.
Here we will stick to the more familiar formulation of superfluid hydrodynamics keeping the Goldstone field as the dynamical field, and pay the corresponding price of having to contend with derivative mixing in the constitutive relations.Putting together (2.15) and (2.19), which shows that the Josephson relation follows from the left-hand side at leading order in derivative, while first-order corrections are (possibly partially) captured by the right-hand side.
The upshot of writing (2.20) is that instead of parametrizing non-hydrostatic corrections to K in terms of (derivatives of) δ B φ, we can equivalently parametrize non-hydrostatic derivative corrections to δ B φ in terms of the combination K − T ∆ µν ∇ µ ns µT ζ ν , which also must vanish in thermal equilibrium.The advantage of doing so is that it avoids the proliferation of time derivatives at first order in gradients, which greatly facilitates the linear response analysis.This is also the reason why we chose to write the first law (2.11) in terms of ζ 2 instead of ξ 2 , as otherwise the combination in the right-hand side of (2.20) would have involved ∇ µ [n s ξ µ /(µT )], which generates time derivatives.This partial choice of frame corresponds to the "modified phase frame" of [17].
Before moving on to determining the non-ideal terms T µν , Jµ , we need to finish fixing our frame, i.e. make a specific choice for our definition of temperature, chemical potential and fluid velocity out of equilibrium.This requires a set of d+2 constraints.In this work, we choose to work with the Landau frame: which constrains the non-ideal corrections to T µν and J µ to be orthogonal to the fluid velocity.In the Landau frame, where we have decomposed T µν = T µν T L + ∆ µν T /d into its traceless and tracefull pieces and denoted −μ the derivative corrections to T δ B φ. σ µν is the shear tensor, which is the transverse, traceless combination of derivatives of the fluid velocity For relativistic superfluids and in the absence of parity violation, there are no hydrostatic corrections to the homogeneous pressure at first order in derivatives, [36].We then parametrize the non-hydrostatic corrections in the following way This matrix must be positive (semi)definite in order to achieve positivity of entropy production.In particular, its diagonal elements must be nonnegative.Its elements can be parametrized as follows: where ζµ ≡ ζ µ / √ ζ ν ζ ν , P µν ≡ ζµ ζν −∆ µν /d is a traceless, transverse tensor and P µν ≡ ∆ µν − ζµ ζν is the projector onto the subspace orthogonal to both u µ and ζ µ .There are a total of 21 non-hydrostatic parameters.The diagonal elements of this matrix are all dissipative, together with the symmetric part of the off-diagonal elements, which is a total of 14 non-hydrostatic dissipative parameters.The 7 parameters left in the anti-symmetric parts of the off-diagonal elements are on the other hand non-dissipative, [29].
Onsager reciprocity requires the matrix appearing in (2.24) to be symmetric, which sets all non-dissipative parameters to zero: We thus obtain a total of 14 dissipative parameters, which is consistent with [18].
Note that in the limit of vanishing background superfluid velocity ζ µ → 0 the unit norm vector ζµ is not well-defined [17], so that one must set β A = α C = α G = α L = 0, while other transport coefficients multiplying terms non-linear in ζ µ (such as eg α B ) also drop out.After imposing Onsager invariance, this means that the only nonvanishing dissipative contributions to the constitutive relations will be parametrized by a total of five coefficients: α A , β L , α F , α H and α S .
Imposing conformal invariance (as we will be interested in doing for holographic superfluids) requires the stress tensor to be traceless, which in turn sets leaving us with a total of 10 dissipative coefficients after imposing Onsager reciprocity, [18].In the limit of zero background superfluid velocity this set is further reduced to three coefficients: α A , β L and α F .

Linear response
We now turn to the study of linear response in the presence of a nonzero background superfluid velocity.That is, we perturb around a state where the spatial component of the superfluid velocity has a nonzero background value The temperature and the chemical potential are perturbed similarly, while the spatial component of the normal fluid velocity is u i = δu i e −iωt+ik j x j .For simplicity, we work in the collinear limit ζ//u//k//ê x , where all these vectors are aligned along the same direction êx .We will henceforth drop the subscript on ζ x = ζ when discussing the magnitude of the superfluid velocity.
We choose to work in the grand-canonical ensemble, holding fixed the temperature and chemical potential, which couple to the energy and charge current densities T tt and J t , respectively.The other sources are u i the spatial component of the normal velocity, which couples to the momentum density T ti .By construction, the Goldstone field couples to K. As we saw in the previous Section, in thermal equilibrium, K = T ∆ νλ ∇ ν [(n s ζ λ )/(µT )].It will be more convenient to consider δ∂ x φ as a fluctuation, which as a consequence couples to T δh x ≡ T δ[(n s ζ x )/(µT )] after integrating by parts the previous relation.This defines a vector of sources δs A = {T δ(µ/T ), δT /T, δu x , T δh x }, which is conjugate to the vector δO A = (δJ t , δT tt , δT tx , δ∂ x φ).
The linear equations of motion with external sources turned off take the form where Here we have labelled the bulk viscosity, ζ b , so as not to confuse it with the spatial superfluid velocity ζ.Interestingly, at nonzero superfluid velocity, the bulk viscosity is non-vanishing even when conformal invariance is imposed.
The dissipative coefficients can be extracted from the retarded Green's functions through the following Kubo relations: To compute the retarded Green's functions and the location of the hydrodynamic poles, we need to work out the matrix of static susceptibilities which relates the O A to the s for which we find and where Here we have often written the individual matrix elements in terms of the ζ derivatives rather than the h derivatives, which will prove more direct to evaluate in holographic states.More explicitly, Both the χ and M matrices are Onsager symmetric, recalling that the eigenvalues under the time reversal symmetry of the O A are respectively {1, 1, −1, −1} and that ξ, ζ should be flipped under time reversal.
With this in hand, the retarded Green's functions can be computed through It is also helpful to trade the energy density for the entropy density, which is conserved at linear level in fluctuations.This can be done using the first law to combine the equations of motion (3.1) and obtain the matrices M and χ in the new basis δO A = (δJ t , δs, δT tx , δ∂ x φ) with sources δs A = {δµ, δT, δu x , δh x }.One can check that the row of the matrix M corresponding to the δs equation matches the linearized canonical entropy current given in (2.16).
The longitudinal sector contains four sound modes, ω = vk + iΓk 2 /2 + O(k 3 ), which are the zeroes of the denominator of retarded Green's functions computed with (3.10).These modes reduce in the zero superfluid velocity case to the so-called first and second sound modes [39][40][41].In general, with a non-vanishing background superfluid velocity, these modes take complicated expressions which are not very helpful to write explicitly.
Clearly, χ ξξ in (3.8) diverges when As we showed in [6], we expect a dynamical instability when this happens.Indeed, one can check that it is at this critical value that one of the sound modes becomes unstable (its imaginary part crosses onto the upper half plane while its real part changes sign).
Around the critical value (3.12), one of the sound modes takes a particularly simple form.Setting ζ = ζ c + δζ + O(δζ 2 ) as before, we find a sound mode ω = vk + iΓk 2 /2 which is linear in δζ: The attenuation constant Γ can be written as: with and Most importantly, the matrix M AB is positive (semi)definite.This follows from positivity of entropy production, applied to matrix (2.24) in the collinear limit.Thus we see that for δζ > 0 (i.e.ζ > ζ c ), the imaginary part of this mode becomes positive, thereby signaling a dynamical instability.

Holographic superfluids
To verify our hydrodynamic predictions, we turn to a class of numerically solvable holographic superfluids [11,12,42].Holographic superfluids are strongly interacting d + 1dimensional quantum systems with a large N number of colors that spontaneously break a global U (1) symmetry.Via the gauge/gravity duality [8,9], at leading order in N , these systems can be described by a d + 2-dimensional classical gravitational action minimally coupled to a charged complex scalar field.A generic action for the system is given by Here AdS ) is the cosmological constant in units of the AdS radius, L AdS , and M = 0, ..., d + 1 runs over all bulk spacetime indices, whereas µ = 0, ..., d runs over only the boundary spacetime indices.
is the field strength for an abelian gauge field under which the complex scalar, Ψ, is charged, In this work, we will only consider the simplest example of holographic superfluids where the scalar potential consists of just a mass term, 3) The first term in the action (4.1) is generically UV divergent.These divergences are regulated by introducing a boundary, ∂Σ, at large r and adding local counterterms via The phase diagram for the holographic superfluid we study (also reported in [6]) as obtained from hydrodynamics and quasinormal modes.The curves on the upper left are slices through the phase diagram for which we report detailed hydrodynamic results.Here, we choose Q = 1 where it is known that there is only one value of the condensate at fixed ζ/µ and T /µ [15].Finally, the star marks the two-stream instability shown in Figure 2.
S bdy .S = S bulk + S bdy is then regular as r → ∞.The procedure leads to a well-behaved variational principle [43,44].The choice of boundary action is dimension dependent and also depends on the choice of m 2 .In what follows, we will restrict to d = 2 and m 2 = −2.
The boundary action in this case is Here, γ M N is the induced metric on the surface ∂Σ, K M N is the extrinsic curvature of the surface, K is its trace, and R (γ) is the Ricci tensor associated with γ M N .Without loss of generality, we can take Ψ = Ψ = ψ.We furthermore set G = 1/16π and L AdS = 1.
The equations of motion which follow from varying (4.1) are The phase diagram of these holographic superfluids was constructed in [6] and is shown in Figure 1.There one can see that the superfluid phase becomes dynamically unstable at the critical superfluid velocity (3.12).Moreover, in some regions of the phase diagram if one increases the superfluid velocity beyond the critical value, a new instability can appear where two of the velocities of the sound modes pick up equal and opposite imaginary components and has been termed the "two-stream" instability [45][46][47].We have illustrated this behavior in Figure 2.  (3.12)where one of the velocities changes sign and is accompanied by a change in sign of the attenuation (here, we just show the velocity but in Figure 6 we also show the attenuations).The second, marked by a star, is the two-stream instability where the velocities become complex and one has a positive imaginary component.In general, this occurs when the first instability is already present, as demonstrated in the phase diagram of Figure 1.

★ ★
In the remainder of this section we will probe that phase diagram and carefully check that the hydrodynamic theory constructed in previous sections describes these holographic superfluids

Thermodynamics and linear response from holography
The equilibrium thermodynamics of holographic superfluids obeys the Landau-Tisza model of relativistic hydrodynamics described in Section 2, as was shown in [15].Here, we ex-tend those results to include the leading order out of equilibrium transport by considering linearized fluctuations of the holographic solutions.
We start by considering the equilibrium thermodynamics.Summarizing [15], variations of the on-shell action give access to the holographic stress tensor as well as the conserved U(1) current.This allows one to match to the non-tilded quantities in (2.13) and (2.14) and obtain the quantities ϵ, p, n, s, n s , µ, ζ and T .By considering a family of solutions parametrized by µ, ζ, and T , we have access to the entire matrix of static susceptibilities.A subset of these equilibrium quantities is illustrated in Figures 3 and 4.
We then consider linearized fluctuations of the gravitational and matter fields with plane wave dependence e −iωt+ik j x j which, on-shell, give access to the retarded Green's functions [48].Having obtained the retarded Green's functions via the holographic dictionary, we then use the Kubo formulae in (3.5) to obtain the dissipative transport coefficients.A subset of the dissipative transport coefficients is illustrated in Figure 5.
In the absence of sources, the equations of motion for the linearized fluctuations only have solutions for specific values of ω and k.These solutions are called "quasinormal modes" and the values of ω(k) for which the solutions exist are the quasinormal mode spectrum.At the same time, via the holographic dictionary, source-free solutions to the linearized equations of motion correspond to poles in the retarded Green's functions.Therefore, the gapless quasinormal modes, which have lim k→0 ω(k) = 0, can be identified with the spectrum of hydrodynamic modes [49].
Hence, to verify that holographic superfluids satisfy the theory of relativistic superfluid hydrodynamics including finite superflow and dissipation, as described in Section 2, we compare the numerically obtained spectrum of quasinormal modes to the prediction for the hydrodynamic spectrum using values for transport coefficients obtained via the holographic dictionary.A representative comparison is shown in Figures 6a, 6b, 7a, 7b, where it is observed that the agreement is excellent. 4n the main text, we consider the fully backreacted holographic superfluid.However, one can also consider decoupling the energy-momentum sector from the charge sector in a superfluid, which is equivalent to considering the case of a probe holographic superfluid.A convenient aspect of probe superfluids is that, for a d = 3 dimensional superfluid, there exists an analytic solution [14,17].This case provides an instructive introduction to the analysis of the main text, so we have included a discussion in Appendix B.  1 that are used to calculate the hydrodynamic dispersion relations.For convenience, in the last plot, we have include the values of µ that we use to construct numerical solutions when r h = 1.

Equilibrium thermodynamics
In this Section, we discuss how to obtain the on-shell constitutive relations for the currents (again, working in d = 2 with m 2 = −2).Choosing a superfluid velocity along the x direction, we write the background metric as and the U(1) gauge field and charged scalar as where r is the AdS radial vector and we choose r → ∞ to be the UV AdS boundary.The choice Ψ = Ψ fixes the phase of the superfluid φ = 0 in equilibrium.Given this ansatz for the coordinate dependence of the fields, we find that there are 5 second order and 2 first order equations for the 7 fields, so that to completely specify a solution we must impose 12 boundary conditions.The boundary conditions are specified as follows.
Finite temperature states of the superfluid are described in terms of black hole spacetimes.In Schwarzschild-like coordinates, this requires that there exists a black hole horizon at some radial coordinate r h where D(r h ) = C xt (r h ) = B −1 (r h ) = 0. Near this surface, the metric reads The matter fields are required to be well-behaved near the horizon-in particular, we require that A t (r h ) = 0. Near r = r h , then, the matter fields have the form, Higher order terms in r − r h are determined in terms of the constants above.The vanishing of D, B −1 , C xt and A t at the black hole horizon constitutes 4 boundary conditions.
At large r, we require that the metric on a fixed r hypersurface is conformal to Minkowski space, which fixes five boundary conditions.After fixing these boundary conditions, the equations of motion lead to an asymptotic form for the metric and matter fields Note that the gauge invariant combination ζ = ∂ x φ − A x fixes the sign of the leading falloff for A x when φ = 0.Here we have chosen a suggestive form for the coefficients of the subleading powers of r.Below, we will vary the on-shell action and observe that ψ s acts as a source for ⟨O ψ ⟩.One can also check that the remaining coefficients indeed correspond to the boundary theory quantities associated with the equilibrium hydrodynamics presented earlier.Interpreting ⟨O ψ ⟩ as the superfluid condensate, solutions with the boundary condition ψ s = 0 describe spontaneous symmetry breaking.The final two boundary conditions can be chosen from the set of constants µ, ζ, T, C h x , C h xt , A h t , A x t , ψ h , s.We will see that µ and ζ act as sources for ρ and n s , so it is natural to choose their values as the last two conditions.This is in line with their interpretation as intensive parameters in the thermodynamics.Finally, r h is determined by a choice of radial coordinate and can be set r h = 1 via a scale transformation r → r/r h .Ultimately, then, after fixing boundary conditions, the equation of state of a holographic superfluid is specified by two numbers, p(T /µ, ζ/µ).
The constitutive relations of the dual stress tensor and the U(1) current are obtained via variation of the renormalized on-shell action [28,43,44,50,51].Using equations (4.1) and (4.4), the variation can be written (4.12) Here n M is the vector normal to the surface ∂Σ, γ M N is the induced metric on this surface, K M N is its extrinsic curvature with K the trace, and R M N and R (γ) are the Ricci tensor and scalar associated with γ M N .
From this expression, we find that matching the standard Landau-Tisza form [15].Of use are the identities To arrive at these identities, one must use the equations of motion to rewrite the bulk action.Setting (S − S bdy ) I on-shell = (S − S bdy ) II on-shell and evaluating for Euclidean time (t = −iτ with τ = τ + β), these enforce the thermodynamic relations The latter of these two relations imposes the conformality condition ⟨T µ µ ⟩ = 0 and leads to p being a function of T /µ and ζ/µ instead of T, µ, ζ independently.
Variation of the Euclidean on-shell action yields the first law.To see this note Using the asymptotics of the bulk fields, this variation can be written for fixed β as [15] In Figures 3 and 4, we present results in the entropy ensemble.To compare to the susceptibilities used in Section 3, one must perform a Legendre transform, leading to the identities where We have now specified how, given a µ and ζ, we fully specify a bulk spacetime solution and how the subleading terms in (4.11) can be matched onto the parameters that appear in the equilibrium currents.By varying µ and ζ, we can obtain an entire family of solutions and obtain the static susceptibility matrix for the holographic superfluid.In Figures 3 and  4, we plot the thermodynamic parameters for a wide range of T /µ and ζ/µ.

Linearized fluctuations
To identify the out of equilibrium constitutive relations, we consider linearized perturbations of the holographic superfluid.We work in radial gauge and consider plane wave perturbations with wavevector pointed along the superfluid velocity,  The equations of motion for these fluctuations are easily obtained but messy, so we omit their explicit expressions here.Nevertheless, we find that 4 of the equations are second order and 4 are first order, necessitating that we specify 12 boundary conditions.We now explain how to do so.
To describe the relaxation back to thermal equilibrium, we care about fluctuations which decay in time.Holographically, this is equivalent to imposing ingoing boundary conditions at the black hole horizon.This amounts to the following near horizon expansions, Not all of these constants are independent, Fixing the radial dependence of the fluctuations to take an ingoing form amounts to 8 boundary conditions.In addition, we see that we are allowed to choose four independent constants, e.g.σ h i , σ h r , g h xx and a h x .This fully specifies the solutions of the linearized equations of motion.
In the UV, the fields have the expansion These expansions imply that the phase field behaves as There are potential terms at O(r 0 ) which can appear if ψ s ̸ = 0, but we will never consider this case.Using this expansion for δφ and (4.12), we identify Expanding the fields near the boundary, we have Variations of the other currents are obtained similarly i ) . (4.31) The variation of the currents contain pieces that are proportional to the leading term in the asymptotic expansion for the fluctuations (i.e.terms with a superscript "(0)").These contribute contact terms in the response functions and are responsible for enforcing the Ward identities.To see this, we note that as r → ∞, the equations of motion impose that 0 = 3µ h t + ka (1)   x From these relations, one can see In the second line, we have identified δψ s = σ (0) r as the fluctuation of the source for the condensate ⟨O ψ ⟩ which explicitly breaks conformal invariance.Furthermore, using the definition of δJ µ and δK gives the variation of the Ward identity δ⟨∇ µ J µ − K⟩ = 0 . (4.34) Finally, using the definitions of δ⟨T µν ⟩, we can check that As we did not consider operators which explicitly break conformal invariance in our hydrodynamic theory, the last term does not appear in (2.9).By choosing δψ s = 0, our results will overlap with the hydrodynamics presented earlier.
While they are important for ensuring that the linearized fluctuations of the holographic superconductor respect the perturbations of the Ward identities, contact terms can differ between various hydrodynamic approaches to obtaining the retarded Green's functions, potentially introducing ambiguities into Kubo formulae. 5In particular, when computing Kubo formulae, we should always consider which will give rise to the same transport coefficient independent of approach.We find that for this quantity, we can neglect the contact terms entirely and instead consider variations of the fields δv I defined by tx , 3h (3)  xx , 3h (3)  yy , −a t , a (1)  x , 2(σ (1)  r + iσ with respect to sources δs I defined by t , a (0) x , σ (0) r + iσ In other words, When working in radial gauge, we have shown that there are only four linearly independent solutions (specified by a h x , g h xx , σ h r , σ h i ) with ingoing boundary conditions which is less than our eight sources, δs I .In order to find the sufficient number of independent solutions, we supplement our ingoing solutions with the pure gauge solutions, which are defined in terms of four functions β M = {β t , β x , 0, β r } M and Λ.We take these to have the following behavior The four pure gauge solutions and the four solutions specified by σ h i , σ h r , a h x , g h xx form eight linearly independent solutions which we label ⃗ X I (r).Here, I = 1, ..., 8 labels one of the solutions and ⃗ X(r) = {a µ (r), h µν (r), σ i (r), σ r (r)} is a vector of the linearized fluctuations corresponding to that solution (in other words, ⃗ X I (r) is an 8 × 8 dimensional object).Since we are only considering linearized fluctuations, a solution, Y I (r), corresponding to a particular configuration of sources δs I in the UV can be obtained by the linear transformation It is clear from this expression that the C IJ are functions of the sources.By considering solutions with all but one source non-vanishing δs I = δ I J , and identifying δv I by expanding the solutions near r → ∞ we can then obtain The transport coefficients are then simply obtained by taking the appropriate limits of this expression, for instance  In Appendix A we give more details on the construction of the correlators and the numerical simulations leading to the results shown in this section.
In Figure 5, we plot the dissipative transport coefficients for a wide range of T /µ and ζ/µ.These results could also be obtained working directly at zero frequency, along the lines of [54,55].
Note that in addition to producing contact terms which enforce the Ward identities, the holographic retarded Green's functions also capture the non-linear physics of the bulk Einstein equations.It is only in the low energy and small wavelength limit, ω/T, k/T ≲ 1, where linear response is appropriate, that we expect the results to match the hydrodynamic expressions.

Quasinormal modes
We have shown how to obtain the thermodynamic and dissipative transport coefficients via holography for holographic superfluids.Furthermore, by varying the quadratic onshell action, we have seen how linearized fluctuations of the currents as obtained from holography manifestly satisfy the linearized Ward identities.It is therefore not surprising that, in the low energy and long wavelength limit where linear response is an appropriate approximation, we can obtain the dispersion relations of the hydrodynamic modes from a holographic calculation as well.
The hydrodynamic modes are poles of the retarded Green's functions whose dispersion relation ω(k) vanishes as k → 0. As we have illustrated that, up to contact terms, the retarded Green's functions are given by (4.43), the poles of the retarded Green's functions are obtained when there exists a solution when δv I ̸ = 0 but δs I = 0 [49].We have similarly argued that, for a given ω and k, there exist solutions δv I with arbitrary sources δs I .However, setting δs I = 0 will generically result in δv I = 0 except for certain sets of ω and k which are called quasinormal modes.These modes generically have complex frequencies, and, in the hydrodynamic limit where k/T ≲ 1, the modes which have |ω(k)|/T ≲ 1 are expected to have dispersion relations which agree with hydrodynamics. 6As expected we have found four hydrodynamic sound modes whose frequencies we denote as ω i , (i = 1, . . ., 4).In Figures 6a, 6b, 7a, and 7b, we show that, indeed, the results agree very well.In particular, in Figure 6b we show that the mode ω 2 becomes unstable at the critical value of the superfluid velocity (3.12) derived in the hydrodynamic theory.For comparison to previous results, e.g.[28], we note that in the limit of zero superfluid velocity ω 1 and ω 4 coincide with superfluid first sound and ω 2 and ω 3 coincide with superfluid second sound.
The recipe for obtaining the hydrodynamic description of holographic superfluids is straightforward in principle, but due to the non-linear nature of the Einstein equations, solutions most often must be found numerically.In Appendix A, we outline the specific way in which we formulated this numerical recipe.0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 -0.15 -0.14 -0.13 -0.12 -0.11 Figure 6: Plots of the real and imaginary components of ω/µ at k/µ = 3.2 × 10 −7 obtained from hydrodynamics (lines) and from quasinormal modes (points) for one of the hydrodynamic sound modes.The color of the lines correspond to the curves in Figure 1.We compare the hydrodynamic results to the quasinormal mode results only on one curve for ease of reading.The mode which signals the instability is ω 2 .We also plot the critical value of ζ obtained from (3.12) in open circles and show that it coincides the sign change of the real and imaginary components of ω 2 .

Discussion
In this work, we have studied the instabilities of relativistic superfluids at finite superflow, reconciling the familiar formulation of such instabilities established long ago by Landau un--27 - obtained from hydrodynamics (lines) and from quasinormal modes (points) for one of the hydrodynamic sound modes.The color of the lines correspond to the curves in Figure 1.We compare the hydrodynamic results to the quasinormal mode results only on one curve for ease of reading.
der the assumption of a weakly-coupled bosonic quasiparticle spectrum, [2]; and more recent results on strongly-coupled superfluids without quasiparticles obtained using gauge/gravity duality, [20], which reported a dynamical instability triggered by the crossing to the upper complex frequency plane of a gapless pole of the retarded Green's functions.More precisely, we have showed that the instability is thermodynamic in nature, and caused by the divergence of the second derivative of the pressure with respect to the (relative) superfluid velocity, (1.2).Assuming a bosonic quasiparticle spectrum, evaluating this quantity and checking when it diverges leads to the usual Landau criterion (1.1), [6].On the other hand, working in the hydrodynamic regime, we also showed that precisely when the thermodynamic derivative diverges, one of the hydrodynamic mode crosses to the upper half plane.This is exactly what happens in holographic superfluids, as we verified in this work.Thus, we have given a unified description of superfluid instabilities, which apply equally well at weak and strong coupling, and as shown in [6], irrespective of any boost invariance.In passing, we have also demonstrated that the dispersion relation of the gapless poles of holographic retarded Green's functions matches the corresponding hydrodynamic predictions at sufficiently low wavenumbers, as expected.
We note that the criterion (1.2) as formulated makes no prediction for the endpoint of the instability.On the other hand, it applies irrespective of its specific mechanism, that is, whether it is caused by exciting rotons in Helium 4 or the nucleation of vortices in holographic superfluids as in [21].In fact, in [21], relaxation of the superfluid velocity through the nucleation of vortices was demonstrated when the initial state was prepared with a superfluid velocity located in the darker grey region of the phase diagram in Figure 1, freezing normal velocity and temperature fluctuations.It would be interesting to include backreaction of these fluctuations, and to probe the black-colored region of the phase diagram.There, sound modes become complex (the 'two-stream' instability discussed in [45][46][47]), the second derivative of the pressure with respect to the superfluid velocity becomes positive again, while instead the heat capacity changes sign.
Another interesting direction for future work is to study the fate of such instabilities at very low temperatures, when the superfluid is coupled to a quantum critical sector, [58,59].In this regime, the dynamical instability may no longer be within the regime of validity of hydrodynamics, if it occurs for values of the superfluid velocity large compared to the temperature.If the critical superfluid velocity is not large compared to the chemical potential, it may still be captured by a zero temperature effective field theory [60], where the derivative (1.2) can still meaningfully be computed.

A Details of the numerics
In this Section we describe the numerical methods used to construct our solutions.We work with a Fefferman-Graham type coordinate u = 1/r with a horizon located at u = 1.We discretize the u coordinate with a Chebyshev-Gauss-Lobatto collocation grid, and solve the resulting equations using a Newton-Raphson relaxation method.All fields are functions of the radial coordinate u only.To allow for temperature variations at fixed µ and ζ, we allow for a variable horizon location u h in our numerics.We work in d = 2 with spacetime coordinates x µ = {t, x, y}.Choosing a superfluid velocity along the x direction, without loss of generality, we can write the metric ansatz as In writing the ansatz in this manner, we have chosen a frame in which the normal velocity, u µ , points only along the time direction.In other words, the boundary conditions on the metric are The Maxwell field in this case has components only along the t and x directions.Boundary conditions on this field are ζ /μ = .109,T/μ = .0023 and regularity of all other fields.The temperature and entropy associated with the black hole are where we have set G N = 1/16π.Note that u h can be set equal to 1 via a scale transformation and as a consequence, the phase space of the background is fully specified by µ/T and ζ/T .Nevertheless, we will keep u h around since it simplifies the derivatives with respect to the temperature, T .
The static susceptibilities are obtained via variations of the equilibrium solutions with respect to µ, T, ζ.We do not need to consider variations with respect to the normal velocity u µ as these susceptibilities are determined via the other variations in terms of Onsager relations.To accurately obtain the variations, we discretize µ, u h , ζ over independent Chebyshev-Gauss-Lobatto grids and solve for the background fields, g M N , A M , ψ as We take these to have plane wave dependence and an asymptotic behavior as u → 0, The retarded Green's functions, GR AB , are obtained by finding four independent sets of solutions, X I for I = 1...4, corresponding to the boundary conditions as follows.Take the largest power ω n appearing in (A.15) and define the extended eigenvector  and matrix In terms of the new matrices, M (0) and M (1) , we are back to a standard eigenvalue problem.
Often, inversion is not necessary to find ω, but the algorithms are still time consuming especially with a large number of grid points.Furthermore, when evaluating the M (i) in terms of numerical solutions, floating point errors can accumulate and fixed high-precision numerics are necessary, causing further slowing.For our purposes, we worked with fixed 60-point precision and monitored numerical accuracy by varying grid sizes.To mitigate this slowdown, another approach is to promote ω to a function of u by introducing the additional equation of motion ω ′ (u) = 0.The eigenvalue problem is then solved instead as a boundary value problem with a Newton-Raphson algorithm.This is a much faster approach when one has a good idea of the form of the quasinormal modes and spectrum and wants to increase accuracy.For instance, we used this method to monitor the convergence of our quasinormal modes as we increased the grid size from N = 40 to N = 400 points at a very low temperature, shown in Figures 8 and 9.The majority of numerical results reported in the main text were obtained at a grid size of N = 120 points, though for low temperatures, we went up to N = 400 points.

B Probe superfluid
In the main text, we considered the case of a superfluid with momentum and energy fluctuations.We can also consider the case where the dynamics of these conserved quantities are decoupled from the charge sector.This simple case is instructive so we have chosen to include it in this Appendix.
The hydrodynamic modes are easily obtained by turning off fluctuations of the temperature and normal fluid velocities.We can then fix u µ = δ µ t .In this case, only three transport coefficients are left, namely σ 0 , ζ 2 and ζ 3 .We can write the dissipative currents in the collinear limit as The spectrum of linearized fluctuations at small k is then given by two sound modes with dispersion relations For larger ζ, v + < 0 and Γ + < 0. This signals a dynamical instability.At the same time we see that the instability is also thermodynamic, for the diagonal susceptibility χ ξξ diverges and changes sign.In [6], we demonstrate that this same criterion coincides with the critical current in superconductors [62].Furthermore, the coincidence of dynamical and thermodynamic instabilities has also been observed in [46].
We now match this to a holographic superfluid.To do so, we take the so-called probe limit, rescaling while taking Q → ∞.In this limit, the complex scalar and gauge fields do not backreact on the geometry, so that we consider their propagation on a fixed background.In most cases, the equations of motion for holographic superfluids must be solved numerically.However, in d = 4, there exists an analytic solution [14,17] when the mass sits at the Breitenlohner-Freedman bound m 2 L 2 AdS = −4 [63,64] and the condensate is very small ⟨O⟩ ≪ 1. Defining u = 1/r, we fix the background metric to be where we have chosen the wavevector to point along the superfluid velocity.Importantly, while we can use a gauge transformation to set Ψ = Ψ in the background, the fluctuations will be independent.As u → 0, the fluctuations can be expanded a t → δA The equations of motion for the fluctuations are three second-order equations for a x , σ r , and σ i and one first-order equation for a t .As such, in addition to fixing the four ingoing boundary conditions, we are also able to fix a h x , σ h r , and σ h i at the horizon, though not a h t which can be found as a function of the other three constants [65].We choose three linearly independent sets of boundary conditions {a h x , σ h r , σ h i } 1 = {1, 0, 0}, {a h x , σ h r , σ h i } 2 = {0, 1, 0}, {a h x , σ h r , σ h i } 3 = {0, 0, 1}.(B.15) which lead to three linearly independent solutions X i = {a t , a x , σ r , σ i } I for I = {1, 2, 3}.Green's functions are obtained in the standard way from linear response [52]-the Green's function for two currents J a and J b is given by the response of the current J a to the source s b of the current J b , As opposed to the main text, we will see there are no contact terms.Here, the sources are the UV boundary values of the fluctuations, δA t , δA x , Σ (0) , and Σ(0) .Clearly, to obtain all Green's functions, we want to independently vary each source.This is not possible with only three independent solutions, however, we are able to add a fourth solution via gauge invariance [65], X 4 = {−iω, ik, 0, ψ}. (B.17) With this solution, any solution for the fluctuating fields can be written and in particular, we can find the solution corresponding to a particular choice of δA x , σ r ,

Figure 1 :
Figure1: The phase diagram for the holographic superfluid we study (also reported in[6]) as obtained from hydrodynamics and quasinormal modes.The curves on the upper left are slices through the phase diagram for which we report detailed hydrodynamic results.Here, we choose Q = 1 where it is known that there is only one value of the condensate at fixed ζ/µ and T /µ[15].Finally, the star marks the two-stream instability shown in Figure2.

Figure 3 :
Figure3: Thermodynamic quantities along the curves 1-10 in Figure1that are used to calculate the hydrodynamic dispersion relations.For convenience, in the last plot, we have include the values of µ that we use to construct numerical solutions when r h = 1.

Figure 4 :
Figure 4: Quantities appearing in the susceptibility matrix along the curves 1-10 in Figure 1 that are used to calculate the hydrodynamic dispersion relations.

Figure 5 :
Figure 5: The dissipative quantities along the curves 1-10 in Figure 1 that are used to calculate the hydrodynamic dispersion relations.Note that η = s/4π.

Figure 7 :
Figure7: Plots of the real and imaginary components of ω/µ at k/µ = 3.2 × 10 −7 obtained from hydrodynamics (lines) and from quasinormal modes (points) for one of the hydrodynamic sound modes.The color of the lines correspond to the curves in Figure1.We compare the hydrodynamic results to the quasinormal mode results only on one curve for ease of reading.

ζFigure 8 :
Figure 8: (Left) Comparison of the components of the susceptibility matrix as a function of grid size N to the values at N = 400 showing exponential convergence.(Right) Comparison of the dissipative coefficients as a function of grid size N to the values at N = 400 showing exponential convergence.The value of ζ + η converges slowest and is the largest source of error in the hydrodynamic prediction.

Figure 9 :
Figure 9: (Top left) A comparison of the real and imaginary components of one hydrodynamic mode as a function of grid size N to the value at N = 400.At N = 400, ω/µ = −(.0721477...)(k/µ) − i(.0580989...)(k/µ) 2 + O(k/µ) 3 (Top right) A comparison of the hydrodynamic prediction for the real part of this mode to the real part of the QNM as a function of grid size N .The points represent Re(ω hydro ) including the value for the dissipative coefficients obtained using Kubo formulae at the same value of N whereas the crosses represent the value for Re(ω hydro ) when the dissipative coefficients are set to zero.The strong convergence in the absence of dissipative coefficients reflects the strong convergence of the susceptibility matrix in Figure 8. (Bottom) Comparison of the imaginary component of the hydrodynamic prediction to the imaginary component of the QNM as a function of grid size N (the dissipative coefficients are non-zero).

Figure 10 : 1 .
Figure 10: (Top left) The dissipative transport coefficients from the Kubo formulae are numerically extracted by extrapolation to ω → 0 (thick line) where the function is constant.At low temperatures and small grid sizes, this may not be possible and larger grid sizes are necessary.Here, we see that N = 120 has good agreement with N = 400 by sight.From Figure 8, they differ by a relative magnitude of ∼ 10 −1 .