Dissipation in Holographic Superfluids

We study dissipation in holographic superfluids at finite temperature and zero chemical potential. The zero overlap with the heat current allows us to isolate the physics of the conserved current corresponding to the broken global $U(1)$. By using analytic techniques we write constitutive relations including the first non-trivial dissipative terms. The corresponding transport coefficients are determined in terms of thermodynamic quantities and the black hole horizon data. By analysing their behaviour close to the phase transition we show explicitly the breakdown of the hydrodynamic expansion. Finally, we study the pseudo-Goldstone mode that emerges upon introducing a perturbative symmetry breaking source and we determine its resonant frequency and decay rate.


Introduction
Holography provides a powerful framework to study large classes of strongly coupled systems at finite temperature. One of its most useful applications concerns the study of phase transitions and broken phases of thermodynamic systems. Moreover, the duality allows us to study large-N systems at all energy scales. This includes the hydrodynamic limit of finite temperature systems which is expected to be universal.
In this paper we will be interested in the long wavelength excitations of the Goldstone mode that emerges in superfluid phases of matter which spontaneously break a global U (1) symmetry. One of our aims is to include the first non-trivial dissipative effects which will result to a finite decay rate for the expected "second sound". The second aim of this paper is to introduce perturbative sources which explicitly break the global symmetry and write down the effective theory of the resulting pseudo-Goldstone mode.
One important aspect of the phases we will consider is that they will be at zero chemical potential and electric charge. This will allow us to isolate the long wavelength dynamics of the condensate avoiding its mixing with the other hydrodynamic degrees of freedom of the system.
Holography has been extensively used to construct and study superfluid phases after their original discovery [1][2][3]. From the bulk point of view, the boundary theory global U (1) is mapped to a gauged U (1) symmetry acting on a complex scalar field. The vector field that gauges the U (1) in the bulk is the gravitational dual of the conserved field theory Noether current operator. Below a critical temperature T c the charged scalar develops a perturbative instability yielding a new branch of broken phase black holes.
The effective theory of the Goldstone mode in the limit of long wavelengths is a topic which has been explored before. The starting point is the non-dissipative fluids which were first described in the framework of the Landau-Tisza theory [4,5]. Dissipative effects at finite temperature were subsequently studied in the literature starting from [6,7] and more recently in [8], in the context of a systematic classification 1 .
Due to the absence of chemical potential and the presence of rotational and time reversal symmetry, the only information we will need to extract from our holographic theories is two transport coefficients. We will be able to express these in terms of horizon data related to finite temperature static configurations. This, is in addition to the charge and current susceptibilities which are fixed at the level of thermodynamics and don't capture dissipative dynamics. A particularly useful tool in extracting these 1 See also [9] for more recent developments in the framework of superfluids. coefficients will be the symplectic current for gravity [10]. At a practical level, this can be seen as the generalisation of the Wronskian for ordinary differential equations. However, it is powerful enough to be applicable in situations where the spacetime is inhomogeneous. Given the explicit expressions for our transport coefficients, we will be able to study our system near the phase transition. As we will see, the speed of second sound will go to zero as we approach the transition and the mode will become diffusive. However, its diffusion constant is not continuously connected to that of the incoherent current pole which exists in the normal phase.
Apart from the purely spontaneous case, we will also consider a scenario in which the global U (1) is an almost exact symmetry 2 . To implement the explicit breaking of the U (1) in a controlled manner, we will introduce a perturbative static source for our complex operator that condenses at low temperatures. As a result, the corresponding current becomes only partially conserved and the propagating pseudo-Goldstone mode acquires a mass as well as a finite gap. Our analysis is based on a modified Ward identity for the now almost conserved electric current and is therefore applicable beyond the details of our specific model.
In section 2 we discuss the class of holographic models we will employ in order to realise our scenario along with some relevant aspects of their thermodynamics. In section 3 we construct the theory of hydrodynamics which captures the dynamics of our Goldstone mode along with the technical aspects of our computation. We conclude the section with the computation of the Green's functions and the Kubo formulae for the two transport coefficients that enter our description. In section 4 we introduce perturbative static sources which break the global symmetry in a controlled manner as well as spacetime dependent ones for the complex scalar operator. We use the sourced Ward identity for the currents to extract the retarded Green's functions of the system. From the poles we extract the dispersion relations of second sound pseudo-Goldstone mode. Section 5 is devoted to numerical checks of some of the analytic results we obtain in the previous sections.

Setup
In this section we will introduce our holographic superfluids. Our main results on the hydrodynamics of the superfluidity Goldstone mode are independent of the specific scenario in which we will realise the neutral superfluid phase transition. For concreteness, we will consider a holographic CFT which is deformed by introducing a source φ (s) for the neutral relevant operator O φ with bulk dual φ. Below a critical temperature T c , which is set by the scale φ (s) , a boundary operator O ψ which is charged under a global U (1) symmetry condenses. From the bulk point of view, this condensation leads to black holes with non-zero profile for the dual bulk scalar ψ.
To realise this setup in holography, we consider a gravitational theory in the bulk consisting of a metric, a neutral scalar φ and a complex scalar ψ which is charged under a local U (1) symmetry and a Maxwell field A µ gauging the symmetry. The local symmetry in the bulk corresponds to a global U (1) symmetry on the boundary theory and the gauge field A µ is dual to the corresponding Noether current operator J µ .
The system is described by the bulk action with the covariant derivative D µ ψ = ∇ µ ψ+iqA µ ψ and the field strength F = dA. We will be mostly interested in configurations which are dual to superfluids and therefore spontaneously break the U (1) symmetry. These backgrounds will have complex scalar with non-trivial modulus and for such configurations it will be beneficial to perform the field redefinition ψ = ρ e iθ . This brings the action to the equivalent form (2. 2) The corresponding equations of motion are In order to implement our holographic scenario, we will assume that for small φ and ψ the potential and the gauge coupling admit the expansions (2.4) Under these conditions, our theory admits the geometry of AdS 4 as solution with metric and trivial scalars and gauge field. Given the expansion (2.4), the bulk scalars φ and ψ correspond to operators O φ and O ψ whose dimensions ∆ φ and ∆ ψ are fixed by the mass terms according to ∆ φ (∆ φ − 3) = m 2 φ and ∆ ψ (∆ ψ − 3) = m 2 ψ . Without loss of generality, a suitable ansatz which captures all the necessary ingredients is This ansatz leads to a non-linear system of ODEs for the radial functions that appear in it. Notice that the above choice of coordinates does not fully fix the radial coordinate r which we are still free to shift by an arbitrary constant. We will choose to fix this freedom by requiring that the horizon of Hawking temperature T is located at r = 0. This allows us to write the near horizon expansion Ultimately, the backgrounds we wish to consider correspond to thermal states of the theory deformed by a source of the relevant operator O φ and a spontaneous VEV for the charged operator O ψ . Given these considerations, the appropriate near conformal boundary expansions at r → ∞ are where we have the leading behaviour as well as the various constants of integration for our system. Given that we are looking for phases where O ψ takes a VEV spontaneously, we will either set the sources ρ (s) equal to zero or, as we will do in section 4, keep them as perturbative deformations. The VEVs of the scalar operators then are fixed by φ (v) and ρ (v) 3 . The constant shift R is an artefact of the way we chose to fix our radial coordinate by placing the horizon at r = 0. At this point, it is worth understanding the asymptotics of the phase field θ. In the absence of (or for perturbatively small) source ρ (s) a perturbation of δθ admits the asymptotic expansion (2.9) In the complete absence of a source for ψ, the relevant bit for us is the constant term in r, which parametrises the phase of the order parameter. Indeed, for the VEV of the condensed operator we will have that O ψ ∝ ρ (v) e i θ (v) . From the boundary point of view, we can simply change the global phase without any cost in energy. The fact that we can freely rotate the VEV O ψ without changing the energy of our system reflects the existence of a Goldstone mode in the boundary theory. This suggests that the dynamics of the Goldstone mode is captured by the VEV of the operator where O ψ is the thermal state VEV. In the absence of a source for ψ, for the fluctuations we can write In section 4 of this paper will be interested in hydrodynamic perturbations in which the source δθ (s) will be much smaller than the scale of δθ (v) . Given this information, in this language the perturbative source for the scalar will be δs ψ = i(ρ (v) δθ (s) +δρ (s) δθ (v) ) while the VEV will be given by O ψ ∝ (2∆ ψ −3) ρ (v) (1+iδθ (v) ) after a perturbation. Therefore the source of the operator O Y will be simply s Y = (ρ (v) δθ (s) + δρ (s) δθ (v) ) while for its VEV we can write equation (2.11).
A central point of interest to our paper is the fluctuations of the bulk gauge field A around the background black holes of equation (2.6). From the equations of motion (2.3), we see that fluctuations of the gauge field A µ and phase field θ are captured by the last two equations. In fact, after defining the one-form field B µ = ∂ µ θ + q A µ , we can simply consider the equation of motion of a massive vector field in the bulk where we have defined the two form W = dB. The penultimate equation of (2.3) now reads and is simply a consequence of taking the divergence of equation (2.12). For this reason we will only need to consider equation (2.12). Close to the conformal boundary, in the absence of background sources for complex scalar, the 1-form field B ν admits the expansion where v α = ∂ α θ (v) + q µ α is a gauge invariant combination of the superfluid velocity ∂ α θ (v) and the source µ α for the U (1) current. Moreover, as we will explain later the constants of integration j α satisfy the Ward identity for the currents which are given by equation (2.20). Note that in the thermal states we are interested in we have µ α = 0; we will only consider non-zero µ α when computing the thermodynamic susceptibilities.

Thermodynamics
An important ingredient in studying the thermodynamics and the linear response of our system is the free energy density w. In order to compute it, we first need to regularise our bulk action (2.1) by introducing appropriate couterterms which render the total action finite and the boundary value problem well defined. These form a boundary action which is defined on a hypersurface ∂M of constant radial coordinate r near the conformal boundary. An appropriate choice of boundary action includes the terms with γ αβ the induced metric on ∂M , R bdr the associated Ricci scalar and K the extrinsic curvature scalar of the constant r hypersurface. In fact, for the case with ∆ φ = ∆ ψ = 2 on which we focus in our numerics section 5, these are all the terms needed.
In order to compute the free energy we need to evaluate the regularised action S tot = S bulk + S bdr on the Euclidean version of our backgrounds (2.6) with t = −i τ .
Since we are dealing with an infinite field theory system, this is an infinite quantity and we should be discussing densities instead. After dropping the integration over the boundary spatial non-compact coordinates, this yields the density I tot which is related to the free energy density w F E via w F E = T I tot with T the Hawking temperature of our black brane. Notice that, in the absence of sources, the terms which involve the complex scalar ψ in (2.15) will not contribute to the thermodynamics of our system. However, these are needed in order to extract the VEV of the complex operator O ψ .
In order to compute the on-shell value of the Euclidean bulk action I bulk , it is useful to write the integrand as a total derivative. After exploiting the existence of the Killing vector ∂ t and using a Komar type of argument, we can write Being in the grand canonical ensemble and at zero background chemical potential, the free energy density of our system is where is the energy density and s is the Bekenstein-Hawking entropy. This can be expressed in terms of the horizon data (2.7) as In the case of e.g. ∆ φ = ∆ ψ = 2 the bulk on-shell action (2.16) combines with the boundary action (2.15) to yield Even though we have not included a source µ α for the current, this is certainly an important part of thermodynanics for a superfluid. The response of the system to an external source includes a non-trivial expectation value for the conserved U (1) current For our purposes, we will need this information in the context of linear response. In addition to the free energy w, we are also interested in the thermodynamic susceptibilities of our system. In order to compute them, we consider the static perturbations δB (t) and δB (x) with zero superfluid velocity. In terms of the asymptotics (2.14), we have δv ν = q δ t ν δµ t and δv ν = q δ x ν δµ x . Close to the conformal boundary these perturbations will admit the expansion In terms of the thermodynamic susceptibilities χ QQ and χ JJ we must have δj t = −χ QQ δµ t 4 and δj x = −χ JJ δµ x . Close to the event horizon, regularity demands the expansion δB (t) The quantities a (0) t and a (0) x are certainly part of the thermodynamics of a superfluid. Later, we will see that together with the horizon data of (3.12) they will determine the transport coefficients of our superfluid when we express the current in terms of the superfluid velocity ∂ α θ (v) .
Later in the paper, we will also want to understand our hydrodynamic expansion close to the phase transition. For the case where the transition is second order, close to the critical temperature This will let us see the important fact that the Goldstone mode does not connect continuously to the diffusive mode of the incoherent current which exists above T c . More interestingly, we will see that the source for this discontinuity is that our hydrodynamic expansion breaks down close to the critical temperature making the radius of convergence infinitesimally small.
By integrating the bulk equations of motion (2.12) we obtain the relations As we can see, the right hand side of the above equations still involves the perturbation of the one-form field. Close to the phase transitions these equations yield This is indeed the behaviour we would have expected close to a superfluid phase transition with the charge susceptibility χ QQ remaining finite and the current susceptibility χ JJ approaching zero.

Linear Response
In this section, we will study long wavelength fluctuations of the U (1) current of our system. The thermal state above the critical temperature is simply an electrically neutral hot plasma corresponding to a CFT that has been deformed by a relevant operator while preserving Poincare invariance. In the absence of electric charge, the electric current fluctuations are dominated by a diffusive mode in the hydrodynamic limit [15][16][17]. The existence of this mode in a neutral plasma is due to the existence of an incoherent current in strongly coupled CFTs [18,19]. Below the critical temperature the hydrodynamics of the electric current will be dominated by the Goldstone mode associated to the spontaneous breaking of the global U (1) of the boundary theory. In section 3.1 we will introduce the symplectic current in the bulk which will be the main technical tool we will use for our computations. In section 3.2 we will construct the hydrodynamic perturbations which couple to the Goldstone mode of our boundary theory. This will allow us to write constitutive relations for the expectation value of the current operator in terms of the phase of our order parameter. Using these and current conservation we will extract the dispersion relation of the second sound mode in terms of thermodynamic susceptibilities and two transport coefficients which are expressed in terms of horizon quantities. Finally, in section 3.3 we will include sources for the current in our analysis. We will give expressions for the retarded Green's functions of the system and Kubo formulae for the transport coefficients of our superfluid.

Symplectic Current
A powerful tool we will use is the symplectic current for the bulk theory. This will allow us to deal with the one-form field fluctuations in an elegant way. Its benefit lies in the fact that it is a conserved current in the bulk and will allow us to relate boundary quantities to horizon ones. As we will see, this logic will allow us to express the boundary electric current in terms of derivatives of the phase of our order parameter, the superfluid velocity, and the sources.
In order to introduce, it we imagine that we have two perturbations δB <1> µ and δB <2> µ which solve our bulk equation (2.12). We can then construct the vector density which can be shown to be divergence free, This current is not a consequence of a conservation law as it is based on two independent solutions. It merely comes from the fact that the equations of motion of our classical theory in the bulk are derived from a local classical action. In the hydrodynamic limit, the dominant term in (3.2) will be the one involving the radial derivative which will let us relate the symplectic current on the black brane horizon to the symplectic current on the conformal boundary boundary. This is especially helpful when one of the two solutions we are using in the construction in the construction of (3.1) is considered to be known. For our purposes, the role of the known solution will be played by the static perturbations δB (t) µ and δB (x) µ we discussed in section 2.1. At a philosophical level, we will need to consider as known all the static black holes which can be used to describe the thermodynamics of our system. That would include the black holes corresponding to finite chemical potential as well as persistent supercurrents. The reason is that, as we will see, the susceptibilities χ QQ and χ JJ will play an important role even though we will be studying transport in thermal states of zero chemical potential and external vector field. For the purpose of extracting those, the knowledge of the static perturbations δB (t) µ and δB (x) µ is sufficient.

Second sound
In this section we would like to construct the bulk perturbation corresponding to the Goldstone mode of the boundary theory: the second sound. As we discussed in section 2.1, the Goldstone mode involves the phase of our order parameter and therefore the associated conserved current due to the global symmetry on the boundary. In the absence of electric charge, the current decouples from the stress tensor of the theory within linear response. This will allow us to clearly isolate the dynamics of the Goldstone mode from the rest of the hydrodynamics modes of our system. From the bulk point of view we will only need to examine perturbations of the one-form field B µ . The spacetime translational symmetries allow us to study Fourier modes. Thus, in order to study perturbations in the hydrodynamic limit, we will consider long wavelength excitations of the form for some small number ε. The function S(r) is chosen so that it drops faster that O(1/r 3 ) close to the conformal boundary and it therefore doesn't interfere with holographic renormalisation. However, close to the horizon, it is chosen so that it approaches S(r) → 1 4πT ln r and the combination t + S(r) is regular and ingoing. Note that we picked the momentum k to point in the direction x without loss of generality given that the background is isotropic. An important ingredient in extracting the second sound is the absence of boundary sources. In general, the asymptotic expansion close to the conformal boundary is The leading terms are a gauge invariant combination of the external one-form source and the superfluid velocity. Absence of sources is equivalent to demanding so that v α = ∂ α θ (v) in the notation of equation (2.14). Moreover, charge conservation implies that At this point it will be illuminating to count constants of integration. This is crucial in order to see that we are doing the right thing in terms of boundary conditions since we are after specific quasi-normal modes. The equation of motion (2.12) yields two second order equations for δb t and δb x while δb r can be solved algebraically. In order to find a unique solution we will therefore need to fix four constants.
Close to the conformal boundary we have the constant δĉ, while δĵ t and δĵ t give one more constant since they have to satisfy the constraint (3.6). Close to the horizon, in-falling boundary conditions fix the expansion We therefore have an extra two constants of integration coming from the near horizon expansion, giving an overall total of 4 constants. However, we are only solving a linear system of homogeneous equations. Since we are not introducing any sources, we can use the scaling symmetry of the problem to set any of the above constants of integration to one leaving with only three. The fourth constant is precisely the frequency ω which will ultimately become a function of the wavenumber k that we are free to choose. This procedure will fix the dispersion relation for ω( k) for the quasi-normal modes. For our purposes however, we are after a particular mode which is hydrodynamic and has ω = 0 for ε = 0. This simply corresponds to the trivial solution for the one form field. This allows us to consider the hydrodynamic expansion (3.8) For the radial function in the ansatz (3.3) we will write the expansion where δB (t) t (r) and δB (x) x (r) are precisely the thermodynamic perturbations we discussed in equation (2.14) with For the subleading terms we impose the near conformal boundary expansion which is simply a reorganisation of the perturbation with boundary conditions as in equation (3.5). Close to the horizon, we impose in-falling boundary conditions. This gives us the expansion, At leading order in the ε expansion the components of the conserved current are which follows from the definition of the charge and current susceptibilities. At that order, we see that equation (3.6) gives with c 2 s = χ JJ /χ QQ being the universal speed of second sound. To specify ω [2] , we will need to fix the currents δĵ α [2] . In order to do this we will employ the conservation (3.2) of the symplectic current (3.1). As we explained in section 3.1, the symplectic current is particularly useful when we compare the solution we are after to a solution we already know. For this reason we will consider it for two different choices of pairs of perturbations. For both cases, we will choose δB <1> to be the perturbation in (3.3).
Our first choice for δB <2> is the static perturbation δB (x) of section 2.1. This will allow us to specify the component δĵ x . Expanding the symplectic current in ε we have One would expect that the leading term would be O(ε) since our perturbation δB <1> starts at order ε and δB <2> starts at zeroth order. However, the leading term in the ε-expansion of δB <1> in equation (3.9) is proportional to δB <2> and the would-be leading term in the symplectic current turns out to vanish. More explicitly we have We now consider the divergence (3.2) and keep only terms up to O(ε 2 ) giving Integrating this equation from the horizon up to the conformal boundary, we obtain the relation where we have defined 5 Choosing now δB <2> to be δB (t) , we obtain 5 For this transport coefficient, a similar result was obtained in [20]. Here we clarify the significance of the quantities that appear in this equation as data related to the geometries dual to the thermal state after the inclusion of the superfluid velocity in the ensemble variables. (3.20) Keeping once again only terms up to order O(ε 2 ) in equation (3.2) and integrating from the horizon up to the conformal boundary, we obtain q δĵ t [2] = iω [2] χ QQ δĉ + e 2g (0) τ (0) δB We have used the near horizon expansion (3.12) and we need to specify the constant of integration δB . We can do this by considering the near horizon limit of the radial component of the equation of motion (2.12) along with the expansion (3.12). This fixes Moving on to next to leading order in equation (3.6) along with (3.23) and (3.18), we have where we have defined (3.26) In addition to the constitutive relations for the currents, the VEV of our charged scalar operator is at leading order in the ε expansion. Here, we have introduced the VEV of the scalar in the thermal state O ψ b = (2∆ − 3) ρ (v) . Using this notation, the dispersion relation for our sound mode is (3.28) Note that these constitutive relations are based on the two thermodynamic susceptibilities χ QQ and χ JJ as well as the two transport coefficients σ d and Ξ. In the next section we will introduce sources for the external vector field which will allow us to compute the retarded Green's functions of the current. As a result, we will manage to write Kubo formulae for our transport coefficients σ d and Ξ. Before moving on, we would like to examine the behaviour of our hydrodynamics close to the critical temperature T c . As one can see from the constitutive relations (3.33), the coefficients that enter in our hydrodynamic expansion include the susceptibilities χ JJ and χ QQ . The behaviour of these coefficients close to T c was given in equation (2.24). The first derivative with respect to the temperature of both of them exhibits the expected discontinuity close to the transition with the former approaching zero. From that, we can easily see that the speed of second sound fixing the linear part of the dispersion relation(3.14) behaves like c s ≈ (T c − T ) 1/2 close to the phase transition.
The dissipative part of the constitutive relations (3.25) is determined by the transport coefficients σ d and Ξ. The holographic expressions that we obtained for these are given in equations (3.19) and (3.26) respectively. The behaviour of σ d is fairly simple to understand as τ (0) goes to a fixed value as we increase temperature and a (0) x is simply unity at T c . On the other hand, equation (3.26) is telling us that Ξ blows up close to the transition since ρ (0) ∝ (T c − T ) 1/2 and everything else remains finite in our holographic expression. We interpret this fact as the breakdown of the hydrodynamic expansion making the radius of convergence smaller and smaller as we approach phase transition. One would have expected such a breakdown since the phase δc becomes a not well defined field close to the transition. This is easy to see directly from the definition of O Y itself in equation (2.10). Note that there is another mode in the theory which becomes gapless close to the transition which is associated with the modulus of the complex VEV O ψ [21]. However, this mode completely decouples from the system we are examining as it is captured by the bulk field ρ in our parametrisation.
Despite this fact, the speed of sound and the attenuation in the dispersion relation (3.28) remain finite since, according to equation (2.24), the current-current scalars behave like χ JJ ≈ (T c − T ). We can easily see that close to the transition the dispersion relation (3.14) becomes which remains finite as we take the T → T c limit. Note that this dispersion relation holds for k T c − T . Above the phase transition where our system is in its normal phase, the U (1) current hydrodynamic excitations are given by the constitutive relations, with δµ the local chemical potential and σ d the incoherent conductivity. We therefore see that above the critical temperature we will have a single diffusive mode determined by the diffusion constant D inc = χ −1 QQ σ d . So, even though our sound modes tend to become diffusive close to the transition, the attenuation constant does not continuously connect to the diffusion constant of the normal phase incoherent mode at the critical temperature.

Green's functions
In this section we will introduces sources for the conserved U (1) current in our system which should be of order δs α ≈ O(ε). The construction of the hydrodynamic perturbation in the bulk is identical to that of equation (3.9). However, this time the frequency is fixed by the external sources and our ultimate goal is to express the boundary current of the system in terms of the sources and the phase of the condensed operator.
The above suggests that we will simply need to replace the boundary conditions (3.5) by with δŝ α ∝ O(ε). Effectively, the whole analysis of section 3.2 goes through after replacing −iω δĉ → δŝ t − iω δĉ and iεk δĉ → δŝ x + iεk δĉ. In terms of spacetime coordinates on the boundary, this is equivalent to replacing the partial derivatives of the phase by the gauge invariant combination ∂ α δc → δs α + ∂ α δc . (3.32) One can of course check this explicitly but here we will just state the final result for the expressions of the components of the boundary current q δj t = −χ QQ (δs t + ∂ t δc) + Ξ ∂ t (δs t + ∂ t δc) , (3.33) After having obtained the constitutive relations (3.33), we can impose the continuity equation (3.6) to obtain an equation for the Goldstone mode. Doing this in momentum space, gives us an expression for δĉ in terms of the sources δŝ α which we can plug back in (3.33). From that linear relation between the VEVs for the current and the sources we read off the retarded Green's functions where we have defined Notice that the last equation in (3.34) is compatible with the Onsager relation given that G is even in k. Moreover, the positions of poles of the Green's functions (3.34) are set by the roots of the denominator of the function (3.35). As expected, these are located precisely on the curves ω = ω ± (εk) set by the dispersion relations we found in section 3.2 for the superfluid sound mode. The final result we would like to present in this section is the Kubo formulae for our transport coefficients σ d and Ξ. By using the Green's functions (3.34) we can write (3.36)

Scalar Sources and Pinning
In this section we would like to explicitly deform our theory by a perturbatively small pinning parameter δρ (s) and study the resulting pseudo-Goldstone mode. In other words, we would like to explicitly break the global U (1) in a controlled fashion. In order to do this, we will have to deform the backgrounds we have considered so far by adding a small perturbative source ρ (s) = δρ (s) ∝ O(ε 2 ) in the near conformal boundary expansion of equation (2.8).
Before applying our logic to the class of holographic theories we are considering, it is worth revisiting our expectations in field theory terms. Suppose that we couple our theory to an external gauge field A α and that we also introduce a source λ for our charged scalar operators O * ψ . If the resulting theory is invariant under the infinitesimal gauge transformations then we can show that the corresponding Ward identity for the current gets modified to This equation makes clear that the order parameter and the source need to align in the complex plane in equilibrium. In other words we can no longer freely rotate the order parameter in the complex plane after fixing a source without generating a current.
In holography, the small background perturbation we want to consider changes the interpretation of the angle θ (v) that we introduced in equation (2.9). To see this, we need to consider the asymptotic expansion of the charged scalar close to the conformal boundary. To do this, we first need to better understand the asymptotics of the vector field B µ to include sources θ (s) for the Goldstone mode. In the absence of background sources for the charged scalar field, we have B α = ∂ α θ (s) (r + R) 2∆−3 + · · · + v α + · · · + q j α (r + R) −1 + · · · , (4.3) where θ (s) , v α and j α are constants of integration, with the equations of motion implying the constraint where we dropped terms of higher order in derivatives of θ (s) . The reason for doing this is that we will be interested in satisfying this equation up to order O(ε 3 ) and such terms would be higher order provided that θ (s) ∝ O(ε 2 ). We will show this later when we consider the proper boundary conditions, compatible with the absence of time dependent sources for the charged scalar operator O ψ . After including the perturbative source terms, the angle θ and the background source term δρ (s) enter the asymptotic expansion of the charged scalar perturbations according to, (4.5) This shows that equation (4.4) is nothing but the Ward identity (4.2) after identifying From the above we see that in order to correctly identify the time dependent perturbative source δs ψ for the scalar field, we need to impose Note that in this notation, the time dependent source of the complex scalar introduces only a source δs Y = δs ψ for the operator O Y that we introduced in section 2. The above equation shows that the source term θ (s) for the perturbation of B µ and the background perturbation source δρ (s) need to be of the same order in the ε expansion. After this observation, the constraint equation (4.4) becomes From the point of view of hydrodynamics, the goal is to have the correct constitutive relations for the currents j α up to second order in ε. That would allow us to satisfy equation (4.8) up to third order in ε provided that we take the source for the charged scalar to be of order O(ε 2 ). This is justified by the fact that the terms we dropped in equation (4.4) are of order O(ε 4 ). This argument shows that the constitutive relations of equations (3.33) and (3.27) are sufficient for this task. Finally, one might worry that we should take in account the fact that the background quantities which enter equation (2.12) will bring their own ε corrections to the perturbation of B µ . However, as we argued above, the correction of the bulk scalar due to the perturbative source δρ (s) will be of order O(ε 2 ). That would induce corrections to our perturbation δB µ at order O(ε 3 ) which is certainly beyond our scope.
In order to compute the retarded Green's functions of the system of our operators, we need to solve the Ward identity (4.8) after identifying θ (v) = δc. By doing so we obtain the explicit expressions where we have defined The poles of the Green's functions in equation (4.9) reveal that the second sound mode has acquired a resonance frequency as well as a gap. More specifically, we find the dispersion relation from which we can read off the resonance frequency ω r = w/χ QQ and the gap 6 ω gap = w Ξ/(2 χ 2 QQ ). It is interesting to notice that, when w < 0, the theory develops an instability. This can be easily seen from the square root in (4.11) which will become imaginary in this case. This instability is simply a mode of the system which wants to align the VEV and explicit deformation of the system in the complex plane.
Finally, we would like to flesh out some of the global effects of the explicit breaking on our observables. In order to do this, we will consider the Green's functions in equation (4.9) at k = 0. In that limit, from the Green's function G O Y J x we see that the transport current J x and the operator O Y completely decouple. Similarly, we will see that J x also decouples from the charge density J t . It is therefore natural to expect that the small explicit breaking will have no effect on the low frequency electric conductivity, However, we see that at k = 0 the charge density and the scalar operator O Y remain coupled turning the electric charge to an almost conserved quantity in the deformed theory.

Numerical checks
In this section we carry out a series of numerical checks of the results derived in the previous sections and in particular, we numerically confirm the sound mode dispersion relation (3.28), the Green's functions (3.34)-(3.35) and the formula for the gap (4.11). We consider the action (2.2) with the following potential and gauge coupling where we note that, for ρ = 0, ∂ φ V = 0 both at φ = 0 and at φ = 1/c. Given the choices above, the corresponding equations of motion admit a unitradius AdS 4 vacuum solution with ρ = φ = A = θ = 0, which is dual to a d = 3 CFT. Placing the CFT at finite temperature corresponds to considering the Schwarzschild black hole, which typically serves as the configuration dual to the normal phase of holographic systems. However, here we choose to deform our boundary theory by a relevant operator, O φ , with scaling dimension ∆ φ = 2. Then, the corresponding backreacted solution dual to the normal phase of our system will be given by black brane with a non-trivial profile for the scalar field φ. As the temperature goes to zero, T → 0, these configurations will approach a flow between the unit-radius AdS 4 in the UV, with φ = 0, and an IR AdS 4 with radius L 2 IR = 12c 2 /(1 + 12c 2 ) supported by φ = 1/c.
To construct these solutions explicitly we consider the ansatz (2.6) with ρ(r) = 0 and the IR and UV boundary conditions, (2.7) and (2.8) respectively, with ρ (0) =ρ (s) = ρ (v) = 0. This boundary condition problem is then solved using a double-sided shooting technique. In figure 1, we plot the logarithmic derivative of the entropy of the system with respect to the temperature, T S (T )/S(T ) for φ (s) = 1 and c = 1. We clearly see that both at very high and very low temperatures the entropy scales like T 2 which is compatible with having AdS 4 on both sides of the RG flow.
On top of these thermal states, we will consider instabilities associated to the scalar field ρ. To ensure that such instabilities exist in our model we need to make sure that the scalar field ρ violates the BF bound associated with the AdS 4 in the IR, i.e.
but is nevertheless stable in the UV: m 2 ρ ≥ − 9 4 . For example, this is the case for c = 1, λ = −3/2 and so we expect an instability to occur for this choice of parameters. Thus, for temperatures below a critical one, we expect a new branch of black holes to emerge characterised by a non-trivial condensate for ρ. To determine the critical temperature at which these instabilities set in we need to study the associated zero mode. In particular, we consider a linearised perturbation around the background constructed above of the form Plugging the above perturbation in the equations of motion, we obtain one second order linear ODE, which we solve by imposing the following boundary conditions at the black hole horizon and asymptotically where we have already set the source for δρ to 0, so that the emergence of the new phase is spontaneous. Overall, the boundary conditions are determined by 2 constants δρ h , δρ v , one of which can be set to 1 because of the linearity of the equation. Consequently, performing the numerical integration will fix the highest value of the temperature, T c = 0.022, for which one can find non-trivial solutions for δρ. The next step is to construct the backreacted solutions corresponding to the broken phase. To achieve this we use the ansatz (2.6) and the IR and UV boundary conditions, (2.7), (2.8). When plugging this ansatz in the equations of motion we obtain a set of three second order ODEs and one first order. Thus, a solution is specified in terms of 7 constants of integration. Looking at the expansion (2.7), we see that it is specified in terms of 3 constants, in addition to the temperature T . Similarly, the asymptotic expansion (2.8) is parametrised by 4 constants in addition to φ (s) -note that we keep ρ (s) = 0 so that the condensation is spontaneous. Overall, in the IR and UV expansions we have a total of 7 constants as well as T and φ (s) , which matches the 7 constants of integration. We proceed to solve this boundary condition problem numerically using double-sided shooting. In figure 2 we show the condensate as well as the free energy as functions of the temperature, T , in support of the phase transition being second order. Here φ (s) = 1 and c = 1, λ = −3/2.

Static perturbations
Having constructed the backreacted black holes corresponding to the broken phase, we now turn our attention to studying perturbations around them. In particular, to check numerically the validity of the analytic expressions for the dispersion relation of the sound mode and the spatially resolved two-point functions, we need to construct the static perturbations as in section 2.1 and extract from them χ JJ , χ QQ . Specifically, considering the perturbations δB (t) t (r) and δB (x) x (r) gives rise to two decoupled linear second order equations, which we solve subject to boundary conditions (2.21), (2.22). Then, given the numerical solutions, the susceptibilities are simply obtained by dividing the corresponding expectation values by the associated sources as extracted by the asymptotic expansion

Second sound and spatially resolved two-point functions
In order to compute the second sound and the two-point functions we need to go beyond static perturbations. In particular, we consider the linearised perturbation (3.3), with Plugging this ansatz in the equations of motion, we obtain two second order ODEs, along with an algebraic equation for δB r . We now turn to the boundary conditions for these functions. In the IR we impose in-falling boundary conditions at the horizon which is located at r = 0 Thus, we see that the expansion is fixed in terms of 2 constant b t , b x 1 , in addition to ω, k. The UV expansion takes the form where δj t is fixed in terms of the other parameters. Overall, this expansion is determined in terms of the sources δs t , δs x and the 2 parameters δc, δj x 1 , in addition to ω, k. For constructing the second sound, we solve the above equations around the numerical background of the previous section, imposing that the sources vanish δs t = δs x 1 = 0. In addition, due to the linearity of the equations we also impose b x 1 = 1. Thus, for fixed k, the shooting method determines b t , δj x 1 , δc, ω. In figure  3, we compare our numerical results with the analytics of the previous chapter by plotting certain derivatives of the dispersion relation. For small values of k, where our analytic arguments are valid, we see a good quantitative agreement.
For constructing the two point functions we again employ a double-sided shooting technique, but now we set either δs t or δs x 1 to zero and scale the remaining source to one using the linearity of the equations. Thus, for fixed ω, k and say δs t = 0, δs x = 1, the boundary condition system determines b t , b x 1 , δc, δj x , allowing one to compute G xx = ω 2 G(k, ω), G tx = G xt = kω G(k, ω), G tt = k 2 G(k, ω) where In figure 4 we plot G (scaled appropriately) as function of the frequency, either for zero and finite k. In the small frequency limit, these quantities approach the constants χ QQ (top, left), −χ JJ (top, right), σ d (bottom, left) and Ξ (bottom, right), in agreement with the analytic expressions (3.34)-(3.35).

Pseudo-gapless modes
In this subsection we outline the numerical computation of the pseudo-gapless modes in the presence of pinning. We perform a calculation similar to the one for second sound, but we now consider linearised fluctuations (3.3), (5.7) around a background configuration that has a small but finite source, ρ (s) , for the scalar ρ. The only difference to the previous subsection is the expansion of the perturbations in the UV part of the geometry. In this case we find that δB t = 0 + δj t r + R + . . . δB x 1 = 0 + δj x 1 r + R + . . . . (5.10) Note that this expansion differs from (5.8), not only because there are no sources for the perturbations in this case, but also because having ρ (s) = 0 pushes the VEV of the goldstone mode to appear at order 1/(r + R) leading to δj t , δj x 1 being independent constants. The parameter counting follows just like above, suggesting that for fixed k we expect to find a discrete set of solutions. In figure 5, we plot our numerical results for the real and imaginary part of the gap and we overlay them with the analytics of the previous chapter. We see a good quantitative agreements for small values of ρ (s) , as expected.

Discussion
In this paper we studied first order dissipative effects in the hydrodynamic regime of holographic superfluid phases of matter. At zero chemical potential and charge density, the normal and the superfluid collective degrees of freedom remain decoupled. This fact simplified our analysis in extracting the transport coefficients relevant to the superfluid dissipation. For a relativistic superfluid, in principle we would only have to determine four invariant quantities that would fully fix the constitutive relation of the conserved current in terms of the phase [6]. In our case, only two of them were non-trivial and they were both fixed in terms of susceptibilities and black hole horizon data. For a specific class of models the coefficient σ d had been computed earlier in the literature [20,22]. Here, we presented a technique which is applicable also in inhomogeneous black hole backgrounds relevant to phases of matter in which translations are broken either spontaneously or explicitly. Interesting extensions of our work include the reduced hydrodynamics of superfluid phases without momentum conservation [23].
One interesting aspect of our work concerns the coefficient Ξ which appears in the time component of constitutive relations (3.33) for the current. We have shown that this grows like (T c − T ) −1 close to the transition signalling the breakdown of the derivative expansion. More specifically, judging from the next to leading term in the hydrodynamic expansion, we have made explicit that this will converge for wavenumbers with k T c − T . It would be interesting to explore the precise way that the hydrodynanic expansion breaks down by following the logic of e.g. [24][25][26]. Finally, we have studied the hydrodynamics of the system upon introducing sources which break the global symmetry in a controlled manner.
An interesting extension we will report on in the future is the inclusion of a background finite magnetic field [23]. That would imply the existence of background vortices which have been previously studied in the framework of holography in [27]. One step further would include the presence of disorder and the resulting flux pinning which relaxes the supercurrent leading to finite DC electric conductivity at finite temperature.
Finally, holography provides access to a plethora of superfluid phase ground states [28,29]. An interesting direction would be to study the low temperature behaviour of the transport coefficients we computed in this paper. This is possible by extracting the behaviour of low temperature black hole horizons by following the techniques of e.g. [30].