Boost-invariant superfluid flows

We present some exact solutions to the ideal hydrodynamics of a relativistic superfluid with an almost-conformal equation of state. The solutions have stress tensors which are invariant under Lorentz boosts in one direction, and represent superfluid generalisations of the Bjorken and Gubser flows. We also study corrections to the flows in first-order hydrodynamics, arguing that dissipation is dominated by the shear viscosity. We present some simple numerical solutions for these viscous corrections. Finally, we estimate the size of corrections to the flows arising when the spontaneously broken $\mathrm{U}(1)$ symmetry responsible for superfluidity is only approximate, giving the corresponding Goldstone boson a small non-zero mass. We find that the massless solutions can still provide good approximations at sufficiently small spatial rapidities.


Introduction
Hydrodynamics is an effective field theory describing long wavelength, low frequency excitations of fluids, and is thus applicable to a diverse range of physical systems. In high energy physics, viscous hydrodynamics has proven to be an excellent framework for describing the experiments on heavy-ion collisions at both the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) [1,2]. Furthermore, it has been claimed that the hydrodynamic approximation may even be helpful in collisions involving protons [3,4].
The hydrodynamic equations are difficult to solve, typically requiring numerical simulation. However, analytical solutions can sometimes be found for flows that preserve a large degree of symmetry. One such solution is the so-called "Bjorken flow" [5]. Bjorken flow arises when the colliding nuclei are modelled as infinite planes, leading to a flow that is invariant under Lorentz boosts along the beam axis, and translationally-invariant transverse to the beam axis.
The latter approximation is of course unrealistic: the radius of a nucleus is on the order of 10 fm, comparable to other scales of typical heavy-ion collisions (for instance, the hydrodynamisation time is of the order of 1 fm/c [6]). A generalisation of Bjorken flow with no translational invariance in the transverse plane was found by Gubser in ref. [7]. However, this solution comes with its own unrealistic simplifications, such as that the initial state and dynamics respect relativistic conformal invariance, and that the collision is perfectly central.
Despite the highly simplifying assumptions of the Bjorken and Gubser flows, they have been extremely useful as a tool for building intuition and as a starting point for understanding more complicated regimes. In addition to providing analytical insight into the hydrodynamic equations, they also serve as a check for the numerical performance of hydrodynamic codes (for example, see ref. [8]).
The aforementioned flows involve arise in the hydrodynamics of a normal fluid. However, if the up and down quarks were massless, quantum chromodynamics (QCD) would have a U(1) × SU(2) L × SU(2) R global symmetry -chiral symmetry -which at low temperatures would be spontaneously broken to U(1) × SU(2) V . In the presence of spontaneously broken continuous symmetries, new hydrodynamic modes appear (the Goldstone modes, which we will refer to as pions), and these should be included in the hydrodynamic equations [9]. The resulting fluid is called a superfluid.
In reality, the up and down quarks are massive, explicitly breaking chiral symmetry. However, since the scale of the explicit symmetry breaking is small, the corresponding Noether current is approximately conserved, and so a suitably modified form of superfluid hydrodynamics should still apply, with possible phenomenological implications. For example, refs. [10][11][12][13] have advanced the critical dynamics near the pseudo-chiral phase transition as a possible explanation of an observed enhancement of soft pion production in Pb-Pb collisions [14][15][16][17][18].
In this paper we will exploit the symmetry arguments used in refs. [5,7] to find analytic solutions to the equations of ideal superfluid hydrodynamics. For simplicity we will consider only the spontaneous breaking of a U(1) symmetry, and thus there will be only one pion. The solutions we find will also apply to systems with spontaneously broken non-abelian symmetries, with the U(1) being the diagonal component of the full broken symmetry group. We leave generalisations of the solutions involving non-abelian effects to future work. We will also consider simple dissipative corrections to the analytic flows, which turn out to require some mild numerics.
Superfluid generalisations of Bjorken flow have been considered previously in the literature [19,20], and we should comment on the differences between our flows and those found in these earlier works. Ref. [19] found an analytic solution for superfluid Bjorken flow in low-temperature QCD by approximating the pion decay constant f by its zero temperature value f π . Accounting for a realistic temperature dependence of f , ref. [19] also found some numerical solutions to the hydrodynamic equations. We extend these results by finding an analytic flow including the leading order low-temperature dependence of f . Ref. [20] included fluctuations of the chiral condensate (the sigma meson) in the list of hydrodynamic degrees of freedom, allowing their theory to apply also above and near the critical temperature. Unfortunately, the resulting equations are sufficiently complicated that only numerical solutions may be found.
In practice, the solutions we find are unlikely to be directly physically applicable to heavy-ion collisions. Both because we neglect fluctuations of the order parameter, and due to the simple equation of state we choose, the solutions will describe superfluid flows with local temperatures well below the approximate chiral symmetry breaking temperature T c . In heavy-ion collisions, for T T c the hadrons chemically decouple [21,22], and at yet lower temperatures they also kinetically decouple [23], leading to a breakdown of the hydrodynamic description. However, we hope that the solutions we present are at least of theoretical interest, and that they may serve as a starting point for a more realistic treatment of the physics at temperatures closer to T c . Perhaps they may also be used to test future numerical hydrodynamic codes that account for superfluid contributions.
The rest of the paper is organised as follows. In section 2 we will write down our model superfluid equation of state, and review the equations of superfluid hydrodynamics. In section 3 we will find both boost-and translationally-invariant solutions to these equations, mimicking Bjorken's treatment for a normal fluid. Flows with no translational invariance will be found in section 4, following Gubser's method. In section 5 we will estimate the size of corrections to our flows arising from small, non-zero pion mass. Finally, in section 6 we will summarise our results and comment on possible directions for future work.

Thermodynamics and ideal hydrodynamics
In a normal fluid the hydrodynamic degrees of freedom are the densities of conserved charges, since they relax only slowly to equilibrium; a conserved charge relaxes over macroscopic timescales since it must spread throughout the system, while other quantities typically relax over microscopic timescales. In a superfluid, the gradient ∂ λ ϕ of the massless Goldstone boson also relaxes over macroscopic timescales, so must be included among the microscopic degrees of freedom, see e.g. ref. [9].
In the superfluids of interest in this work, the conserved currents are the stress tensor T µν and a U(1) current J µ . Their equations of motion are the continuity equations In ideal superfluid hydrodynamics, the currents take the form [9] T µν = εu µ u ν + p∆ µν + µf 2 (u µ ξ ν + ξ µ u ν ) + f 2 ξ µ ξ ν , where u µ is the normal fluid velocity, ∆ µν ≡ g µν +u µ u ν (with g µν being an arbitrary metric) is a projector orthogonal to u µ , ξ µ ≡ ∆ µν ∂ ν ϕ is proportional to the relative superfluid velocity, and ε, p, µ, ρ and f 2 are the energy density, pressure, chemical potential, charge density, and pion susceptibility, respectively. In addition, we must impose the Josephson condition on the pion, which in ideal hydrodynamics reads 3) The Josephson condition may be motivated in a number of different ways: it arises from the canonical conjugacy of the pion and the charge density [9]; it may be derived by demanding that the entropy current is conserved at the ideal order [24,25]; finally, it arises as a natural consequence of a higher form symmetry associated to winding of the superfluid phase [26,27].
To complete the ideal hydrodynamic equations, we must specify an equation of state p(T, µ, ξ), where ξ ≡ ξ λ ξ λ is real since ξ λ is spacelike (as it is by definition orthogonal to the timelike vector u µ ). The entropy density s, charge density ρ, and pion susceptibility f 2 are determined through the first law dp = s dT + ρ dµ − f 2 ξ dξ, while the energy density is given by ε = T s + µρ − p. For clear discussions of the thermodynamics of superfluids, see refs. [9,28].
Throughout this work, in order to obtain concrete solutions to the equations of superfluid hydrodynamics we will use a model equation of state of the form 1 The first term describes a normal fluid with conformal equation of state. This is one of the simplifying assumptions mentioned in the introduction that limits applications to QCD, in which the equation of state at µ = ξ = 0 takes a rather more complicated form [30,31]. 2 The simplicity of the T 4 term will be crucial in allowing us to derive exact solutions to the hydrodynamic equations. The second and third terms in equation (2.4) represent the leading corrections to the pressure at small non-vanishing chemical potential and superfluid velocity. The restriction to small chemical potential is merely a simplifying assumption, but it will at least turn out to be self-consistent; when we come to construct boost-invariant superfluid flows with this equation of state, we will find that they require us to set µ = 0. On the other hand, the restriction to small ξ is physically motivated. At large ξ superfluidity breaks down due to the creation of vortex filaments [33,34]. This has been observed in superfluid helium, for example [35].
We will take the susceptibilities in equation (2.4) to have the form 3 for some constants (f π , k, k ). These are the leading order behaviours of the susceptibilities at low temperature in O(N ) linear sigma models, for example [36]. Since we are now also making a low temperature expansion, we expect the model equation of state (2.4) with susceptibilities (2.5) to be appropriate for superfluid flows with the hierarchy of scales Our equation of state is a special case of the more general equation of state for nuclear matter discussed in ref. [29], obtained by treating the normal fluid component as being conformal. 2 A conformal normal fluid may be more appropriate for models with chiral symmetry-broken phases that are deconfined, such as can occur in the holographic Sakai-Sugimoto model [32]. 3 Our motivation for the notation fπ is that in QCD, f (T = 0) = χ(T = 0) is the pion decay constant.
In the QCD literature the susceptibility f is usually referred to as the spatial pion decay constant, and is often denoted fs, while χ is the square of the temporal pion decay constant ft. See for example ref. [36].
The equation of state that we have chosen contains only one dimensionful parameter, f π , and so for f π = 0 the superfluid is conformal. Indeed, the energy density following from the equation (2.4) is and we find the trace of the stress tensor to be T µ µ = −ε + 3p + f 2 ξ 2 = −f 2 π (∂ϕ) 2 , which vanishes when f π = 0. Notice that this expression for T µ µ is precisely what one would obtain for a free, massless scalar field ϕ with action S = −f 2 π d 4 x √ −g (∂ϕ) 2 . We will see that f 2 π plays very little role in the hydrodynamic equations, and so in section 4 we will be able to leverage the powerful conformal-symmetry based method of ref. [37] to construct boost-invariant superfluid flows with non-trivial energy density profiles, even when f π = 0.

First-order hydrodynamics
The first-order dissipative corrections to superfluid hydrodynamics are rather involved: the most general constitutive relations of first-order superfluid hydrodynamics contain twenty transport coefficients [38]. We will consider only parity invariant superfluids, eliminating six of these coefficients [38]. Further, the assumption of small relative superfluid velocity discussed above eliminates a further nine, leaving only the five transport coefficients discussed in ref. [33]. Adopting the notation of ref. [24] for the dissipative corrections, one then obtains the constitutive relations 4 where T µν ideal and J µ ideal are the ideal currents given in equation (2.2), {η, κ, ζ 1,2,3 } are the five transport coefficients, 2u (µ ξ ν) ≡ u µ ξ ν + ξ µ u ν , and The five first-order transport coefficients are not completely free. Demanding that the entropy current has postive divergence imposes that η, κ, ζ 1 , ζ 3 ≥ 0 and ζ 2 2 ≤ ζ 1 ζ 3 [24,39,40]. Further, taking the trace of the stress tensor in equation (2.8), one finds We thus see that a conformally invariant theory has ζ 1 = ζ 2 = 0. In section 3.2 we will use this fact to argue that the shear viscosity η is the most important first-order transport coefficient for the physics of superfluid Bjorken flow.

Ideal superfluid
The natural coordinate system for studying Bjorken flow may be obtained starting from four-dimensional Minkowski space in cylindrical coordinates, where we will refer to the z direction as the beam axis, while (x ⊥ , θ) are polar coordinates in the transverse plane. One then defines the Milne time τ and spatial rapidity w through t = τ cosh w and z = τ sinh w. We use the symbol w for the spatial rapidity, rather than the more conventional η, in order to avoid confusion with the shear viscosity. The Milne time is invariant under boosts along the beam axis, while the rapidity is shifted to where v is the boost velocity. In the new coordinates, the metric is We now seek solutions of the ideal superfluid hydrodynamic equations that are invariant under boosts in the z direction, rotations about the z axis, and translations perpendicular to z. Together, these conditions imply that observables will depend only on τ . Let us take the normal fluid velocity to be u = ∂ τ . The combination of the Josephson condition (2.3) and symmetry then demand that ξ = ξ w (τ )∂ w . The fact that ∂ λ ϕ = ξ λ + µu λ is a gradient implies that ∂ τ (τ 2 ξ w ) = 0, so we may take ξ w (τ ) =ξ/τ 2 , for some constantξ. If we assume that µ = µ(τ ), then we can straightforwardly solve ξ λ = ∂ λ ϕ − µu λ for the Goldstone boson field, ϕ =ξw − dτ µ(τ ).
We can now evaluate the continuity equations (2.1). Conservation of the U (1) charge and the momentum, ∇ λ J λ = 0 and ∇ λ T λw = 0 respectively, imply that For non-zeroξ and µ, this implies that ∂ τ (χ/f 2 ) = 0. Since χ and f 2 are generically different functions of the local temperature [36,41], we conclude that superfluid Bjorken flow is usually consistent only for µ = 0. Concretely, for χ and f 2 given by equation (2.5), equations (3.3) imply that µ = 0 unless k = k . We will therefore set the local chemical potential to zero for the remainder of this section, satisfying the charge and momentum conservation equations (3.3) trivially. The solution for the Goldstone mode is then ϕ =ξw. A notable example of a theory with k = k is actually chiral perturbation theory, where it has been shown that f 2 π and χ differ by a term O(T 4 ) at low temperatures [42,43]. Thus, at the level of the low-temperature approximation used in equation (2.5), the equation of state from chiral perturbation theory permits superfluid Bjorken flow with non-zero µ. Unfortunately, we have not been able to find an analytic solution for the µ = 0 case. 5 This case is studied numerically in ref. [19].
With µ = 0, the only non-trivial conservation law is the conservation of energy, ∇ λ T λτ = 0, which implies The general solution is where τ 0 is an integration constant. Substituting this solution into equation (2.7) we obtain the corresponding energy density, At late times, this energy density decays as The first term in the expansion is the usual Bjorken flow energy density, with superfluidity providing subleading corrections. Notice that the zero temperature pion decay constant f π appears only in the first subleading correction to the late-time energy density, which decays as τ −2 . We then obtain a series of further corrections decaying as larger powers of τ −1 , controlled by the finite temperature correction to the susceptibility f 2 . The term in equation (3.7) proportional to τ −2 is reminiscent of a similar term found in the superfluid Bjorken flow of ref. [19], however it has a different origin. Their flow haŝ ξ = 0, but a non-zero chemical potential µ ∝ τ −1 (allowed because they set f 2 = χ = f 2 π ). This leads to an energy density containing a term ∝ µ 2 ∝ τ −2 .
Expanding equation (3.5) at small proper time τ → 0, one finds that the local temperature grows linearly as T (τ ) ≈ 4nτ /(kξ 2 τ 2 0 ) for non-zeroξ. Thus, the system appears to heat up with τ at early times, and then cool down at late times. However, this may just be an artifact of our definition of the local temperature. In particular, no such heating is apparent in the energy density, which also decreases at early times as Finally, we make a brief comment on the behaviour of the U(1) current J µ in our solutions. In the fluid frame, the charge density u µ J µ = 0, while there is a non-zero current in the w-direction that decays with increasing τ . The situation is rather more intricate in the "laboratory frame" (3.1). In cylindrical coordinates (t, z) the current is with T given by equation (3.5) (and τ 2 = t 2 −z 2 ). At late Minkowski time t the contribution from T 2 is negligible, and the current (3.9) takes a rather simple form J = f 2 π (z ∂ t + t ∂ z )/(t 2 − z 2 ). For fixed t, we see that the current density diverges as J t → ±∞ as the lightcone at z = ±t is approached. The solution thus describes two shockwaves of equal and opposite infinite charge density travelling outwards at the speed of light. The total charge in the z ≥ 0 region increases with t, due to the non-zero J z . At early and intermediate t, the behaviour of the current is rather more intricate due to the temperature dependence in equation (3.9). In particular, the charges of the shockwaves can flip sign at intermediate times.

Dissipative corrections
We would now like to study dissipative corrections to the flow in equation (3.6). As for the ideal case, the assumed symmetries force observables to be functions only of τ . We continue to take the normal fluid velocity to be u = ∂ τ , and consequently the relative superfluid velocity is still of the form ξ = ξ w (τ )∂ w . With this ansatz, two of the five transport coefficients discussed in section 2.2 have no effect on the flow: it is straightforward to show that when T and µ are functions only of τ , Comparing to the first-order constitutive relations in equation (2.8), we see that ζ 3 and κ drop out of the hydrodynamic equations completely.
We are then left with three transport coefficients: the shear viscosity η, the bulk viscosity ζ 1 , and ζ 2 . In principle, we would need to do a microscopic calculation to compute how these transport coefficients depend on T , µ and ξ. However, as noted in section 2.2, for a conformal superfluid ζ 1 = ζ 2 = 0. Although our superfluid is not conformal, the only dimensionful parameter is f π , which drops out of all equations when µ = ξ = 0. We therefore expect ζ 1,2 = 0 in this limit. 6 In the case of interest, namely µ = 0 and small non-zero ξ, we thus expect the shear viscosity η to be much larger than ζ 1,2 . For the remainder of this section, we will therefore set ζ 1,2 = 0 and concentrate on the effects of the shear viscosity. We will further approximate the shear viscosity by its conformal form η ≈ηT 3 whereη is a dimensionless constant, neglecting any corrections at non-zero ξ.
The net effect of these approximations is that we are keeping terms up to O(ξ 2 ) in ideal hydrodynamics, while throwing away everything except the normal fluid contribution to first-order hydrodynamics. In principle, given formulas for η, ζ 1 , and ζ 2 at non-zero ξ, there is no barrier to extending the treatment in this section to account for dissipative superfluid corrections to the flow. However, we expect these corrections to be doubly small in the physically interesting regime, being suppressed in both the small-ξ and derivative expansions. With µ = 0 and the only non-trivial transport coefficient being η =ηT 3 , the Josephson condition implies ϕ =ξw for some dimensionless constantξ, as in the ideal case, while the only non-trivial conservation equation is the conservation of energy ∇ λ T λτ = 0, which implies 1 This equation is solved implicitly by the solution to the algebraic equation 12) where τ 0 is an integration constant, chosen such that the leading order behaviour of T (τ ) at late times matches that of the inviscid solution (3.5). We cannot solve this transcendental equation (3.12) for T (τ ) analytically, but it is straightforward to obtain numerical solutions, for example with the Newton-Raphson method. Substituting the resulting temperature into equation (2.7), we then obtain numerical results for the energy density as a function of τ . For example, in figure 1 we show results for the energy density for various different values ofη, for the following sample values of the parameters: n = π 2 /30, k = 1/6,ξ = 1/10, and τ 0 = 1/f π . 7 At small τ the curves appear to coincide in figure 1a at ε ≈ f 2 πξ 2 /2τ 2 , due to the dominance of the final term in the energy density (3.6). At large τ the curves coincide again, as can be seen from the following argument. We can solve equation (3.11) perturbatively in 1/τ to determine the late-time behaviour of the local temperature. Substituting the result into equation (2.7), we find that the energy density decays at late times as ε(τ ) = 3n We see that the leading order correction to ε(τ ) due to non-zero shear viscosity occurs at the same order in the late-time expansion as the leading-order correction from non-zero ξ.
Finally, we note that the superfluid component regularises a singularity present in the Bjorken flow of a viscous normal fluid. The normal fluid local temperature, obtained by solving equation (3.11) withξ = 0, is (3.14) As pointed out in ref. [7], for example, this local temperature becomes negative at early times, when the second term dominates over the first. This is not a huge problem, it is just a symptom of the expected breakdown of hydrodynamics at early times when gradients become large. However, it is perhaps pleasing that no such singularity exists in the 7 The choices of n and k were motivated by O(N ) sigma models, for which n = (N − 1)π 2 /90 and k = (N − 2)/12 [44]. We obtain the chosen values of n and k by taking N = 4, as relevant for QCD. The value k = 1/6, which is also found in chiral perturbation theory [42,43]. The qualitative behaviour of ε(τ ) does not depend sensitively on the choices of n and k.  superfluid case, as visible from the lack of a zero in the subtracted energy densities plotted in figure 1b. Concretely, solving equation (3.11) perturbatively at small τ , we find that at leading order the local temperature scales linearly in τ as The proportionality coefficient is (perhaps contrary to first appearance) real and positive for all realη. However, notice that at intermediate τ there is a very large deviation between thē η = 0 curve in figure 1 and those corresponding to all but the smallestη = 0. The zeroand first-order terms in the hydrodynamic expansion are thus comparable in this regime, so we cannot trust hydrodynamics along the whole flow except for small values ofη.

Review of the normal fluid case
In the previous section we described a boost-invariant superfluid flow that is also invariant under translations transverse to the beam axis. For real applications, for example to heavyion collisions, this translational invariance is unrealistic, since the fluid will be expanding from some central collision region.
In ref. [7], Gubser describes an elegant method for deriving simple, transversely expanding, boost-invariant flows of normal, uncharged fluids with conformal equation of state p = ε/3. Gubser's method works by replacing the transverse translational invariance with a generalised conformal invariance. In this subsection we will review Gubser's method. In section 4.2 we will show how the method may be applied to superfluids.
The first step is to define a vector field [7] where q is a constant with dimensions of inverse length, a is the vector generating a unit translation in the θ = 0 direction, and b is the vector generating a unit special conformal transformation in the same direction. Thus, ζ is a conformal Killing vector of the Minkowski metric (3.2), where L ζ is the Lie derivative along ζ. Given a tensor field T ν 1 ν 2 ... µ 1 µ 2 ... , we say that it transforms covariantly under ζ if it satisfies where the weight α[T ν 1 ν 2 ... µ 1 µ 2 ... ] is a real constant. For instance, equation (4.2) implies that the metric has weight α[g µν ] = −2. The inverse metric g µν correspondingly has α[g µν ] = 2, ensuring that the Kronecker delta δ µ ν = g µλ g νλ has vanishing weight. Gubser's flows are constructed by demanding that the fluid velocity and energy density transform covariantly under ζ. Applying L ζ to the normalization condition g µν u µ u ν = −1 and using α[g µν ] = 2, it is straightforward to show that the fluid velocity must have weight α[u µ ] = −1, so that Assuming rotational invariance around the beam axis -implying that u θ = 0 -and boosting to a frame where u w = 0, the solution to the differential equation (4.4) is [7] (4.5) Having fixed the fluid velocity, we now turn to the conservation of the stress tensor. If we demand that the energy density transforms covariantly under ζ with weight α[ε] ≡ α ε , this implies [7] ε =ε (g) whereε(g) is a so far arbitrary function. The variable g is the only independent scalar combination of τ and x ⊥ satisfying L ζ g = 0, up to trivial redefinitions, so any scalar with vanishing weight must be a function of g only. Assuming a conformal equation of state p = ε/3, at the level of ideal hydrodynamics the stress tensor is T µν = ε (4u µ u ν + g µν ) /3. With the fluid velocity given by equation (4.5) and the energy density taking the form (4.6), conservation of the stress tensor implies that the energy density satisfies 3(1 + g 2 )ε (g) + 4g(α ε − 2)ε(g) + 4qτ (α ε − 4)ε(g) = 0. (4.7) Plainly this equation is only consistent with the assumption thatε(g) is a function only of g and not τ if α ε = 4. This value may also be obtained from dimensional analysis: ifε(g) is dimensionless, then α ε = 4 correctly assigns the energy density dimensions of (length) −4 . We may straightforwardly solve equation (4.7) with α ε = 4, and substitute into equation (4.6) to find the energy density [7] ε(τ, x ⊥ ) =ε whereε 0 is an integration constant.
Ref. [7] also computed the first-order dissipative correction to equation (4.8). Since the fluid is assumed to have conformal equation of state, the bulk viscosity vanishes and so the only non-zero first-order transport coefficient is the shear viscosity η, which contributes to the stress tensor as in equation (2.8). Further, by dimensional analysis, the shear viscosity in the conformal fluid must take the form η =ηε 3/4 for some dimensionless coefficientη. The stress tensor therefore reads with the fluid velocity given by equation (4.5). Conservation of this stress tensor leads to a modified version of equation (4.7), with solution [7] ε(τ, , (4.10) where we leave implicit that g is the function of τ and x ⊥ appearing in equation (4.6).

Superfluid case, ideal hydrodynamics
We will now show how to generalise the conformal symmetry based-approach reviewed in section 4.1 to a superfluid with equation of state (2.4). We will once again seek flows in which all physical quantities transform covariantly under ζ as in equation (4.3). The procedure for determining the normal fluid velocity is unchanged, so that u is still given by equation (4.5). We then fix the superfluid velocity ξ by demanding that it transforms as in equation (4.3) with weight α[ξ µ ] ≡ α ξ . Combining this with the condition that u µ ξ µ = 0 that arises from the definition ξ µ ≡ ∆ µν ∂ ν ϕ, one finds We now impose the condition that ∂ λ ϕ = ξ λ + µu λ is a gradient, which implies that ∂ µ (ξ ν + µu ν ) − ∂ ν (ξ µ + µu µ ) = 0. If we assume that µ transforms with weight α[µ] ≡ α µ , so that µ =μ(g)/τ αµ , then the various components of this gradient condition imply (α ξ − 2)ξ(g) = 0, (α µ − 1)μ(g) = 0,ξ (g) = 0. (4.12) Thus, non-zero solutions forξ andμ are only possible if α ξ = 2 and α µ = 1, respectively, and the only possible solution forξ is that it is constant. Now we can attempt to solve the conservation equations. We will assume that the local temperature transforms covariantly under ζ with weight α T , so that we may write T =T (g)/τ α T . As for the Bjorken flow described in section 3, we find that simultaneous conservation of J λ and T λw is impossible unless the local chemical potential vanishes. Concretely, taking a linear combination of the equations ∇ λ J λ = 0 and ∇ λ T λw = 0, we find the condition (k − k )ξT (g)μ(g) (1 + g 2 )T (g) + α T (g + qτ )T (g) = 0. (4.13) This equation makes clear that for the generic case k = k , flows with non-trivial superfluid componentξ = 0 and non-trivial local temperature are only possible ifμ(g) = 0. We then find that the conservation equations for the remaining components of the stress tensor are only consistent with the assumption thatT (g) is a function only of g if α T = 1, which is also what one would expect from dimensional analysis. In this case, conservation of the stress tensor implies (1 + g 2 ) 12nT (g) 2 + kξ 2 T (g) + 2g 4nT (g) 2 + kξ 2 T (g) = 0. (4.14) Equation (4.14) may be solved explicitly forT (g). The resulting expression for the local temperature T =T (g)/τ is whereΣ is related to the function appearing in equation (3.5) bỹ where τ 0 is an integration constant, chosen such that T reduces to the Bjorken flow result (3.5) in the q → 0 limit. Substituting the solution for the local temperature into equation (2.7) with µ = 0, we then obtain the energy density Since this is a rather complicated function of τ and x ⊥ , we plot the typical dependence of ε on x ⊥ of the energy density in figure 2, for several values of τ . The plots were made by setting n = π 2 /30, k = 1/6, q = f π ,ξ = 1/10, and τ 0 = 1/f π , the same values as in figure 1. In both the superfluid and normal fluid cases, the energy density exhibits a single peak as a function of x ⊥ , which at early times qτ < 1 is located at x ⊥ = 0. For qτ > 1 the peak is located at qx ⊥ = q 2 τ 2 − 1. 8 Note that since x ⊥ is the radial coordinate orthogonal to the beam axis, the peak at non-zero x ⊥ is really a ring of energy density surrounding the beam axis, that expands outward with a speed asymptotically approaching the speed of light as τ → ∞. Away from the peak, the superfluid energy density is dominated by the x ⊥ -independent contribution to (4.17) proportional to f 2 π . The form of the late-time decay of the energy density depends on the value of x ⊥ . At the peak at qx ⊥ = q 2 τ 2 − 1, the energy density takes the relatively simple form where x ⊥ q ε/q 4 Figure 2: Sample logarithmic plots of the energy density of ideal superfluid Gubser flow as a function of transverse distance x ⊥ , at different Milne times τ . To make these plots, we have fixed the parameters of the model and the solution to n = π 2 /30, k = 1/6, q = f π , ξ = 1/10, and τ 0 = 1/f π . The solid blue curves show the energy density (4.17). The dashed orange curves show the same energy density (4.17), but with the x ⊥ -independent contribution f 2 πξ 2 /2τ 2 subtracted, in order to make the x ⊥ -dependence clearer. The dotted grey curves show the normal fluid result (4.8), withε 0 chosen such that the two results agree at qτ = 0.1 and x ⊥ = 0, for comparison. Away from this peak, for large τ at fixed x ⊥ the energy density decays as 9 The dominant term in both cases is proportional to f 2 π /τ 2 and is independent of x ⊥ , consistent with the behaviour plotted in figure 2. We see that first subleading correction at late times decays more rapidly away from the peak.

Dissipative corrections
We now consider the effect of dissipative corrections on the superfluid Gubser flow of the previous subsection. The symmetry arguments fixing the velocities are unchanged from the ideal case, such that u is given by equation (4.5), and ξ = τ −2ξ ∂ w with constantξ. With this form of ξ µ one finds ∇ µ f (T ) 2 ξ µ = 0, so that the transport ceofficient ζ 3 makes no contribution to the hydrodynamic equations for the flows that we are considering. We will consider only flows with µ = 0, so the coefficient κ also makes no contribution. 10 As argued in section 3.2, of the three remaining transport coefficients {η, ζ 1 , ζ 2 } we expect the shear viscosity η to be most important for dissipation at small relative superfluid velocity. Moreover we expect that the shear viscosity will be dominated by the normal fluid contribution, such that it may be approximated by η ≈ηT 3 for some dimensionless constantη. Making this approximation, and neglecting the contributions of ζ 1,2 to the hydrodynamic equations, we find that the only independent hydrodynamic equation is We are not able to solve equation (4.21) analytically, however it is straightforward to solve numerically. To do so, we first solve the equation perturbatively for large values of g, findingT (g) = n q 2 τ 2 0ξ 2 k where τ 0 is an integration constant, chosen such that the leading term in this expansion matches the leading term in the large-g expansion of the ideal solution (4.15). We then use this expansion to set boundary conditions at some large value of g, and then numerically integrate to smaller values of g. The resulting solution for T =T (g)/τ can then be substituted into equation (2.7) to obtain the energy density. For example, in figure 3 we show the result of this numerical procedure for a sample value of the dimensionless shear viscosityη = 1, together with the inviscid solutions shown 9 Note that the expansion (4.20) diverges in the limitξ → 0, despite the fact that the full expression for the energy density in equation (4.17) is finite in this limit. The reason for this noncommutativity of limits is that the functionΣ(g) appearing in equation (4.17) tends to zero in the τ → ∞ andξ → 0 limit, at a speed that depends on the order of limits. 10 If one were to consider a flow with µ = 0, transforming with weight αµ = 1 as required by dimensional analysis, one would find ∆ µν ∇ν µ T = 0 for the fluid velocity in equation (4.5). As a result, κ again drops out of the hydrodynamic equations. in figure 2. The rest of the parameters in the model were set to the same values as in the latter figure, namely n = π 2 /30, k = 1/6, q = f π ,ξ = 1/10. For theη = 0 curves we used τ 0 = 1/f π , while for theη = 1 curves we tuned τ 0 ≈ 74/f π , such that the value of the energy density at the local maximum (as a function of x ⊥ ) matches theη = 0 result at late times. From the figure, we see that non-zero shear viscosity tends to lead to reduced energy density at early times, and much more rapid decay of the energy density at large x ⊥ . These two effects are related by the conformal symmetry: from the form of g in equation (4.6) we see that both small τ and large x ⊥ correspond to large positive g.
As for Bjorken flow, we find that the superfluid component regularises the pathology pointed out in ref. [7]. Whenξ = 0, equation (4.21) may be solved to obtain [7] T (g) = qε with integration constant ε 0 , leading to the energy density (4.10). This solution forT (g) becomes negative at large positive g, corresponding to small τ or large x ⊥ , signalling a breakdown of hydrodynamics in these regimes since the viscous contribution outweighs the ideal fluid result. No such behaviour is visible in the superfluid energy density plotted in figure 3, where a change in the sign ofT would lead to a zero in the subtracted energy density (the dot-dot-dashed purple curve in the figure).

Comments on introducing a pion mass
In QCD, chiral symmetry is only approximate, giving a small mass m to the pion ϕ. To what extent do the massless flows found in sections 3 and 4 provide good approximations to the true flows for small non-zero m? A necessary condition for the validity of the approximation is that the kinetic term in the pion Lagrangian dominates the mass term, (∂ϕ) 2 m 2 ϕ 2 , when evaluated on the m = 0 solution. All of our massless solutions have ϕ =ξw, and so this condition is equivalent to m 2 w 2 τ 2 1, where the factor of τ 2 arises from g ww = τ −2 in the metric (3.2). For given values of m and w, this implies an upper bound on τ , beyond which we cannot trust the massless-pion approximation. This is dangerous, since we would typically only trust hydrodynamics at late times where we expect gradients to become small. For instance, when the pion is the real pion of QCD with m ≈ 1.4 fm −1 in natural units, we expect the massless-pion approximation to hold for where we have restored a factor of the speed of light c. Since 1 fm/c is of the order of the time scale for hydrodynamisation in heavy-ion collisions [6], one cannot hope to accurately describe the regions with |w| ∼ O(1) by setting m = 0. Instead, we expect that the massless pion approximation will only be accurate for sufficiently small values of the rapidity w.
We will now estimate the size of massive corrections to superfluid Bjorken flow in more detail. Throughout this section we consider only ideal hydroydnamics, for simplicity. A small pion mass modifies the conservation equation of the U (1) current to ∇ µ J µ = f 2 m 2 sin ϕ, (5.2) while the Josephson condition and conservation of T µν are unchanged. Using the constitutive relations written in equation (2.2) with ρ = χµ and ξ λ = ∂ λ ϕ − µu λ the conservation equation becomes where we have used the Josephson condition to replace µ with −u α ∂ α ϕ.
The dependence of f 2 and χ on the local temperature couples this equation to the conservation equation for the stress tensor. We have not been able to find explicit solutions to the resulting set of coupled equations for the hydrodynamic variables, so for simplicity we will approximate the susceptibilities by their zero temperature value f 2 ≈ χ ≈ f 2 π . The conservation law (5.3) then becomes simply the sine-Gordon equation This approximation is not very interesting from the point of view of hydrodynamics, since we have essentially decoupled the pions from the hydrodynamic flow. However, it is hopefully good enough to estimate the size of massive corrections to our solutions. We will solve equation (5.4) perturbatively at small m by expanding the pion field as whereξw is the m = 0 solution found in sections 3 and 4, and ϕ 1 is the leading-order correction at non-zero m. Substituting this expansion into equation (5.4) and discarding terms O(m 4 ) and higher, one finds The general solution to this sourced wave equation is where W ± (τ ±w ) are left-and right-moving waves along the beam axis. 11 Note that the non-linear dependence of ϕ 1 on w indicates that boost-invariance is broken by non-zero m 2 . The relative size of the O(m 0 ) and O(m 2 ) contributions to ϕ depends on the explicit form of the functions W ± . We will motivate a particular form of these functions by demanding that the contribution of the pion to the stress tensor is approximately boost-invariant 11 Moving to the cylindrical coordinate system (3.1) by setting τ = √ t 2 − z 2 and w = 1 2 log t+z t−z , one finds τ e ±w = t ± z. near w = 0. Concretely, substituting the expansion (5.5) into the sine-Gordon stress tensor T µν = f 2 π ∂ µ ϕ∂ ν ϕ − 1 2 g µν (∂ϕ) 2 + m 2 g µν cos ϕ , we find the non-zero components and T x ⊥ x ⊥ = x 2 ⊥ T θθ = −T τ τ . If we impose approximate boost invariance near w = 0 by setting ∂ w T µν | w=0 = 0, this yields the conditions ∂ τ ∂ w ϕ 1 | w=0 = ∂ 2 w ϕ 1 | w=0 = 0. Substituting the solution (5.9) into these conditions yields two second-order ordinary differential equations for W ± , which may be solved to obtain where c is an integration constant. 12 If we further impose the physically reasonable condition that ϕ → 0 as t → ∞ for fixed z, then we find c = 0 in equation (5.9). Then, expanding ϕ 1 (τ, w) for small w one finds m 2 ϕ 1 (τ, w) ≈ −ξm 2 τ 2 w 3 /6. The massless-pion approximation will be self-consistent only if this is much smaller in magnitude than the m = 0 solution ϕ 0 =ξw, and so we recover the estimate m 2 τ 2 w 2 1 made at the beginning of this section. We emphasise that this conclusion depends strongly on the functions W ± , however. For instance, if one were to choose boundary conditions that set W ± = 0, and also assume thatξ 1, then to leading order at smallξ one finds instead |m 2 ϕ 1 /ϕ 0 | ∼ m 2 τ 2 .

Discussion
In this work we have presented a number of boost-invariant solutions to the equations of motion of superfluid hydrodynamics, for a superfluid with a simple, yet physically motivated equation of state (2.4). In the superfluid generalisation of Bjorken flow presented in section 3, we found that the superfluid component of the flow leads to subleading terms in the energy density at late times beginning at O(τ −2 ), the same order in the late-time expansion as viscous effects appear.
In section 4, by applying Gubser's conformal symmetry based approach [7], we found flows expanding transverse to the boost direction. The superfluid component of the flow provides a contribution to the energy density that is independent of the transverse direction and dominates at late times, plus additional subleading corrections to the flow at large 12 There is only one integration constant in equation (5.9) because we have set boundary conditions to eliminate a term which is independent of τ and linear in w -since such a term can be trivially removed by a redefinition of the O(m 0 ) integration constantξ -and a constant term independent of both τ and w that may be removed by the shift symmetry present when m = 0. transverse direction. The former contribution is rather physically unsatisfying, and in future work it would be pleasing to find a modification of the flow in which the superfluid contribution to the energy density is also expanding.
To find solutions to the hydrodynamic equations, we have treated the spontaneously broken symmetry giving rise to the superfluid as exact, so that the resulting Goldstone bosons are massless. In QCD chiral symmetry is only approximate, and thus the pions have non-zero masses, with important consequences for transport. In section 5 we estimated the size of massive corrections to the above flows, concluding that the massless-pion approximation should hold for the region of the flow at small spatial rapidities w 1. A more complete treatment of our flows would include the pion masses throughout. However, one would need to take into account the complication that a non-zero pion mass is incompatible with boost invariance.
There are many other possible directions for further work. For instance, one could study second-order dissipative effects, or perturbatively study anisotropy in the transverse plane (as would be introduced by non-zero impact parameter in a heavy-ion collision) using the methods of ref. [45].
Throughout this work we have considered only superfluids arising from the spontaneous breaking of a U(1) symmetry. As discussed in section 1, our solutions apply also to nonabelian superfluids, with our U(1) being a subgroup of the full non-abelian symmetry. 13 A natural question is whether a non-abelian superfluid permits more general boost-invariant flows with multiple pions excited.
One might attempt to generalise superfluid Bjorken flow along the lines of ref. [37], in which a complex time translation was applied to the Bjorken stress tensor, resulting in a complexified stress tensor T µν C . Since the hydrodynamic equations are invariant under time translations and linear in the conserved currents, T µν ≡ Re T µν C remains an exact solution. For an appropriately chosen translation, the new stress tensor is approximately that of Bjorken at small rapidities, while providing a better fit to data at large rapidities. Applying the aforementioned procedure to the stress tensor of superfluid Bjorken flow would lead to a conserved T µν , however, it seems unlikely that this would be the stress tensor of a superfluid: since the superfluid stress tensor in equation (2.2) is non-linear in ξ µ , the relative superfluid velocity that would be read off from this new T µν would probably not be expressible as ξ µ = ∆ µν ∂ ν ϕ for some scalar ϕ.
Finally, as mentioned in section 1, possibly the largest barrier to the physical interpretation of our results is the restriction to temperatures much smaller than f π , clashing with the chemical and kinetic freezeouts that occur at low temperatures in heavy-ion collisions. It is thus desirable to extend our treatment to larger temperatures by adding new terms to the equation of state (2.4). Such terms would further break conformal symmetry, so we do not expect it to be possible to repeat the analysis of section 4. It should be straightforward to extend the translationally-invariant case discussed in section 3, although solving the resulting equations may require numerics. An additional complication is that sufficiently close to T c one must also take into account fluctuations of the order parameter, taking one to the regime studied in ref. [20].
We hope that the solutions presented in this work, and perhaps some of the generalisations that we have suggested, provide helpful intuition for superfluid effects in highly relativistic hydrodynamic flows, and that they may be a useful tool for benchmarking future numerical codes.