Incoherent hydrodynamics of density waves in magnetic fields

We use holography to derive effective theories of fluctuations in spontaneously broken phases of systems with finite temperature, chemical potential, magnetic field and momentum relaxation in which the order parameters break translations. We analytically construct the hydrodynamic modes corresponding to the coupled thermoelectric and density wave fluctuations and all of them turn out to be purely diffusive for our system. Upon introducing pinning for the density waves, some of these modes acquire not only a gap, but also a finite resonance due to the magnetic field. Finally, we study the optical properties and perform numerical checks of our analytical results. A crucial byproduct of our analysis is the identification of the correct current which describes the transport of heat in our system.


Introduction
The holographic conjecture predicts that in a certain large-N limit, large classes of conformal field theories possess a dual classical gravitational description. Apart from its fundamental implications about quantum gravity, it provides a powerful tool to study strongly interacting regimes of quantum field theories which are inaccessible by standard perturbative techniques.
Over the last decade, the duality has been used to study aspects of strongly coupled systems. One of its exciting applications concerns condensed matter systems at finite temperature, chemical potential and magnetic field [1][2][3][4][5][6][7][8]. In that context, the discussion was sparkled by the discovery of electrically charged black hole instabilities which lead to superfluids/superconductors [9][10][11] from the field theory point of view. In this case, the order parameter is given by the expectation value of a complex operator which breaks an internal U (1) symmetry.
Soon after the discovery of holographic superfluid phases, black hole instabilities which spontaneously break translations were found in [12]. These phases are expected to play a crucial role in understanding particular physical aspects of various condensed matter systems which exhibit instabilities such as charge and spin density waves, including the cuprate superconductors. In this paper we wish to construct the effective theory of long wavelength excitations in holographic phases with spontaneously broken translations.
In order to make contact with realistic condensed matter systems, one needs to tackle the extra complication of the ionic lattice which relaxes the momentum of charge and energy carriers in the system. Momentum relaxation is an essential ingredient in discussing the low frequency transport properties of real materials. In order to accomplish this holographically, we need to deform our UV conformal field theory by relevant operators with source parameters which break translations. In other words, apart from the spontaneous, we also need to implement explicit breaking of translations.
The construction of these inhomogeneous black hole backgrounds and the study of the corresponding thermodynamics is technically challenging mainly due to the fact that unstable modes naturally lead to inhomogeneous backgrounds where the only expected symmetry left is time translations. However, for certain classes of holographic theories with a bulk action which is invariant under global U (1) symmetries one can follow a Q-lattice construction [13] in order to implement both the holographic lattice as well as the order parameter that spontaneously breaks translations by simply solving ODEs. This system was introduced in [14,15] where the focus was on the transport properties and the derivation of analytic formulae for the DC transport coefficients.
In [16] it was shown that the spontaneous breaking of the global U (1) in the bulk introduced additional diffusive hydrodynamic degrees of freedom to the system, which are separate from the universal ones associated to the conservation of heat and electric charge in the system. From that point of view, the system we are studying is different from the modulated phases of holography where apart from translations, no additional symmetry breaking occurs. In this work, our aim is to generalise the results of [17] in order to include an arbitrary number of internal broken symmetries as well as a finite magnetic field. 1 Similar systems with spontaneous breaking of translations via an internal symmetry have been studied before in [21]. In the absence of an explicit lattice and at zero magnetic field, the longitudinal hydrodynamic modes included one pair of sound and two diffusive modes. One of the diffusive modes can be accounted to the incoherent thermoelectric mode while the second one was associated to the diffusive mode of the internal symmetry breaking described in [16]. Similarly, the transverse sector of the system contains a single pair of sound modes. It is known that a finite magnetic field has the effect of combining the transverse and longitudinal sound modes to produce a gapped mode and a mode whose frequency was growing quadratically with the wavenumber [22,23]. Based on the hydrodynamic models of [24,25], reference [26] argued that the corresponding constant of proportionality is complex, using numerical techniques in a holographic massive gravity model; see also [27] for related work in effective theories for weak explicit background lattices.
In contrast, in our work, we consider phases with strong explicit background lattices, and we analytically show that all hydrodynamic modes remain diffusive even in the presence of a magnetic field. Another aspect of our effective theory is the inclusion of explicit deformation parameters which perturbatively pin the density waves in the system. As one might expect, such deformations introduce a collection of gaps for some of the diffusive modes in our theory. Interestingly, we find that at finite magnetic field pinning also introduces resonance frequencies. As we show, from the retarded Green's functions point of view these show up as poles in the lower half plane.
In section 2 we discuss the class of holographic models we are considering along with some important aspects of their thermodynamics. In section 3 we start by introducing the model of hydrodynamics that provides an effective description of the long wavelength excitations, meanwhile identifying the correct current that describes the transport of heat in our system. We then move on to include the effects of pinning in order to compute the resulting gap and resonance of the density waves, as well as compute the retarded Green's functions and extract their optical properties. We conclude the section by discussing how to decouple the Goldstone modes from the U (1) and heat currents, and by deriving the dispersion relations of our hydrodynamic modes. Section 4 contains a number of non-trivial numerical validity checks of the effective theory of section 3. We summarise our most important observations and conclude in section 5. Finally, the appendix contains technical details of the analytical calculations of section 3.

Setup
In this section we introduce the holographic model that captures all the necessary ingredients that we would like to include in our theory. For this reason, we consider holographic theories which in addition to the metric, they contain a gauge field A µ and N Y + N Z complex scalar fields Y J and Z I with a global U (1) N Y +N Z symmetry. Essentially, this is a generalisation of the model that we considered in [17] to include an arbitrary number of complex scalars in the bulk.
The gauge field will be used to introduce the chemical potential µ and the magnetic field B in the dual field theory. The first N Y complex scalars are going to implement the explicit lattice and should therefore be relevant operators with respect to the UV theory. The remaining N Z will provide the density wave order parameters in our system. For simplicity, we consider only four bulk spacetime dimensions, corresponding to a 2 + 1 dimensional conformal field theory on the boundary, but all our results can easily be generalized to higher dimensional theories as well.
The class of theories we are considering is described by the two-derivative bulk action 2 with G I , W J , τ and V being only functions of the moduli b I = |Z I | 2 and n J = |Y J | 2 . In this case, the global internal symmetries are represented by the invariance under the global transformations Z I → e iθ I Z I and Y J → e iω J Y J . The equations of motion are In order for our bulk theory to admit an AdS 4 solution of unit radius, we demand the small Z I and Y J expansions In this case, the conformal dimensions ∆ I and∆ J of the dual complex operators For the AdS 4 vacuum we use a coordinate system in which the metric reads In these coordinates, the near conformal boundary expansion for the scalars takes the form Provided that the operators O Y J are relevant with∆ J < 3, the Q-lattice construction [13] picks the deformation parameters y J s (t, . At this point we see that we have a set of N Y wavevectors k J si which are determined externally as part of the sources related to the explicit breaking of the lattice. In our theory, the operators O Z I take a VEV spontaneously suggesting that for the bulk fields Z I we should have z I s = 0. In this case, the bulk field Z I are going to be zero above a certain critical temperature. As we lower the temperature of the system, they start developing instabilities which will yield new branches of black holes breaking the global bulk N Z U (1)'s spontaneously. However, in section 3.3 we will turn on z I s perturbatively in order to study the pinning of the density waves.
Specifically, within the Q-lattice ansatz for the backgrounds, we have that up to potential contact terms. Thus, we have a set N Z wavevectors k I i associated to the order parameters O Z I . These are dynamically chosen by the system in a way that the free energy is minimised, in contrast to k J si which are fixed by us as part of the boundary conditions. For the particular holographic systems we are studying, the free energy is minimised when k I i = 0. However, following the logic of [17,28], we will still consider background solutions with k I i = 0. In this way the order parameters of the spontaneous breaking also break translations apart from the internal U (1)'s.
It is useful to define the real operators which have a constant expectation value The above suggests that from a microscopic point of view, the operator is the right object to focus on in order to study the gapless fluctuations of the system. In order to make this point clearer, we parametrise the spacetime fluctuations of the VEVs O Z I according to s ]. In order to solve the bulk equations of motion more efficiently, we find convenient to make the field transformations bringing the bulk action to the form A consistent ansatz that captures all the necessary ingredients we discussed so far for the thermal state is (2.14) According to (2.12), the constants c I in our ansatz (2.14) for the background black holes translate to an overall phase for the complex scalars Z I . The absence of explicit sources in our asymptotic expansion for the corresponding field does not fix it and we have to leave it arbitrary. These are essentially the gapless modes associated with the symmetry breaking in the bulk that we wish to promote to hydrodynamic ones in section 3. We choose our coordinate system so that the conformal boundary of AdS 4 is approached as we take r → ∞. In this case, the asymptotic expansions of the functions in our ansatz (2.14) take the form where we chose to only show the terms where constants of integration of the relevant ODEs appear. The constants of integration g (3) ij that appear in the expansion of the metric have to satisfy δ ij g which is the gravitational constraint and yields the conformal anomaly for the stress tensor.
The constant of integration R in (2.15) represents a global shift in the radial coordinate r. This is fixed by demanding that the finite temperature horizon is at r = 0. Demanding our background solutions to be regular imposes the near horizon expansion The equations of motion (2.2) lead to the following equations for the phases of the complex scalars associated to spontaneous breaking At this point, it is useful to note that the fields σ J and χ I are not well defined when either ψ J or φ I are equal to zero. This certainly happens close to the conformal boundary and in order to avoid misinterpretations with the holographic dictionary, we discuss asymptotic expansions in terms of the complex fields Y J and Z I through (2.12). This is well defined in the regime of perturbation theory that we are interested in. The asymptotic expansions for the perturbations of φ I and χ I are (2.18) From these expansions and using (2.12) we obtain where we have used thatz I s = 0, or equivalently φ I s = 0, in the phases we are interested in. This shows that, up to contact terms, δ S I = ∆ I − 3 2 φ I v δc I , and that 2 δφ I s is a source for Ω I while 2 ζ S I is a source for S I , consistent with equation (2.11) and the discussion below it.
In the next subsection we discuss aspects of the thermodynamics of our broken phase black holes. This will give us the opportunity to define certain quantities that will appear later in the context of hydrodynamics.

Thermodynamics
In this subsection we would like to consider the thermodynamics of the background black holes we are interested in. In order to do this we need to regularise the bulk action (2.1) by adding suitable boundary terms which act as counter-terms [29,30]. The purpose of these terms is dual, the first is to render the total on-shell action finite. The second is to make the variational problem well defined, provided we have a unique way to fix the boundary conditions on our bulk fields.
It is often the case that such terms are not unique and for the purposes of our paper it is enough to list the following terms [30] (2.20) Further counter-terms can be added [31] but these will introduce extra contact terms in the retarded Green's functions that we wish to compute from the bulk theory. In order to compute the free energy of the system we need to consider the Euclidean version of the total action I E = −iS tot . We then need to evaluate the value of I E on the solution with the analytically continued time t = −it E and the periodic identification t E ∼ t E + T −1 . Since our system extends infinitely in the spatial field theory directions, the total free energy W F E = T I E is not meaningful and we instead consider the free energy density w F E where , s and ρ denote the energy density, entropy density and electric charge density respectively. Apart from the thermodynamic data T , µ, B and the explicit lattice data ψ J s , k J si our solutions also depend on the wavenumbers k I i which are related to the spontaneous breaking. Even though different values of c I in (2.14) yield different solutions, the free energy is independent of those in the spontaneous case, when φ I s = 0. 4 The first variation of the free energy yields the first law where the electric charge and entropy densities can be computed from the black hole horizon data and M is the magnetisation and by τ (0) we denote the value of τ when φ I and ψ J are evaluated on the horizon. In order to compute the variation w i I of the free energy with respect to the wavenumber k I i , we simply vary the total action S tot and use the background equations of motion to find 5 Notice that in the spontaneous case this is coming entirely from the variation of the bulk action (2.1) and it is a finite number. Potential contributions from finite counter-terms other than those listed in (2.20) would be possible but these would constitute contact terms which we can ignore for our purposes.
In order to obtain the electric magnetisation M ij of the dual field theory, we need to perform a straightforward variation of the bulk action (2.1) with respect to the magnetic field B. Apart from the electric magnetisation, our backgrounds are also going to have a non-trivial thermal magnetisation M ij T . Similarly to the electric magnetisation, this is not immediately obvious from the backgrounds in (2.14) since the homogeneity of our solutions prevents the appearance of explicit heat magnetisation currents. In order to define it, we would need to introduce a larger background ansatz than (2.14) in order to include NUT charges in our metric [35]. Instead of doing that, we simply give the expression for its value in terms of the background Apart from the quantities that appear in the first variation of the free energy, we also find it useful to introduce a set of susceptibilities through the second variation of the free energy An expression for the susceptibilities ν i I , β i I and w ij IL in terms of the background can be obtained by simply varying e.g. equation (2.24). Moreover, the susceptibilities ν i I and β i I can also be found by varying the densities in equation (2.23) with respect to the spontaneous wavenumbers k I i . Even though it is not obvious from the bulk expressions that these two approaches lead to the same result, this is guaranteed by the thermodynamic Maxwell relations. This observation will become important in section 3, when we derive the constitutive relations for the currents and the Josephson equation in a derivative expansion.
In order to discuss the constraints which we need to impose, we choose to work in a radial foliation and we define the normal one form n = dr normal to constant r hypersurfaces. We now use our equations of motion in (2.2) to define L µ = E µ ρ n ρ and C = C ρ n ρ with E µν = L µν − 1 2 g µν L ρ ρ . The four gravitational constraints are simply L µ = 0 and the Gauss constraint is C = 0.
The constraints that we need to satisfy can be imposed on any hypersurface with e.g. constant radial coordinate r. Close to the conformal boundary, they are equivalent to the Ward identities of charge conservation, diffemorphism and Weyl invariance of the dual conformal field theory. Similar to [17], we will derive an effective hydrodynamic theory in section 3 by utilizing a subset of these constraints on a surface close to the background black hole horizon at r = 0, in terms of the constants that appear in (2.30) and (2.31).
We now define the horizon currents Building on [35,36], we derive the horizon constraints that the above currents should satisfy 6 Close to the conformal boundary at r = ∞, the expansion of our functions is where we have chosen to only show the undetermined terms involving the constants of integration of the second order ODE's we need to solve. We have also included the constants ζ S I which represent the sources for the operators S I we discussed in section 2 and which we need to fix. For the components which we are free to choose by using diffeomorphism and gauge invariance, we only show their desired behaviour close to the boundary as they offer no additional information as far as the constants of integration are concerned. The remaining 9 + 2N Z + 2N Y constants of integration δg , δc I and δσ J(v) together with the 13+2N Z +2N Y coming from the 6 It is relatively straightforward to check that the r and t components of the gravitational constraints are equivalent in the r → 0 limit. horizon expansion, will fix a unique solution of the 9+2N Z +2N Y second order ODE's and the 4 constraints. In the next section we construct solutions which correspond to the late time, long wavelength hydrodynamic modes of the boundary theory.

Linearised hydrodynamics
In this section we study the hydrodynamic limit of the fluctuations that we introduced in subsection 2.2. In subsection 3.1 we discuss the construction of these modes up to second order in the derivative expansion. This allows us to write an effective hydrodynamic theory for the conserved currents of the system and its gapless modes related to spontaneous breaking in the bulk in subsection 3.2.
Having a complete effective theory for our fluctuations, in subsection 3.3 we examine the gap induced for the spontaneous density waves sliding modes, by perturbatively small deformations for the operators Ω I . To illustrate, we specialise to the isotropic case where we find that apart from a gap, the theory also develops resonance frequencies. In subsection 3.4 we study the retarded Green's functions of the operators in our theory at finite frequency and we give the precise way that the poles of subsection 3.3 influence the transport properties.
We then move on to give more general, model-independent, Kubo formulas for some of the transport coefficients in subsection 3.5. There, we also define heat and electric current operators which decouple from the Goldstone modes, and we discuss some of their properties. Finally, in subsection 3.6 we give an algebraic equation whose solutions yield the dispersion relations of our hydrodynamic modes. Even though we are not able to find the solutions in closed form, we can show that all our modes are purely diffusive at the order we are working. An interesting outcome of our results is that, after correctly identifying the heat current, the thermodynamic coefficients w i I all drop out from physically interesting quantities.

Hydrodynamic perturbations
In their infinite wavelength q i → 0 limit, the hydrodynamic modes we wish to study will reduce to a uniform distribution of energy, electric charge and phase shift of the complex scalar Z I which break the global symmetries in the bulk. In order to keep track of our expansion we scale q i → εq i with ε a small number and expand the frequencies and radial functions in the bulk according to where δX(r) can be any of the functions that appear in the ansatz for the perturbation (2.28).
A key point of our construction is the leading piece of the ε expansion which according to our earlier discussion has to reduce to The functions X b represent the background fields of the black hole in equation (2.14).
In order to generate the perturbations which satisfy the correct boundary conditions (2.30), (2.31) and (2.34), at the same time with a simple partial derivative with respect to T , µ and c I we also need to perform the perturbative coordinate and gauge transformation [36] As expected, the boundary condition requirements for our perturbations do not uniquely fix g(r) in the bulk. It is enough to choose it such that close to the conformal boundary it vanishes sufficiently fast while close to the horizon at r = 0 it approaches Choosing the perturbations as in (2.28) leads to an inhomogeneous system of differential equations coming from the bulk equations of motion (2.2) at the perturbative level. It is clear from the above construction that the seed solution (3.2) satisfies the corresponding homogeneous system of equations [17]. This suggests that we can add them at each order in the ε expansion (3.1) and therefore consider the split, with δX [n] a solution to the inhomogeneous problem which is sourced by lower order terms of the solution. Following closely the analysis of [17], we can show that the eigenmodes of the system necessarily have ω [1] = δT [0] = δµ [0] = 0. Therefore the variation of the temperature and chemical potential starts at order O(ε). The next to leading part of the bulk solution δX [1] will only be driven by a shift in the background phases of the complex scalars according to δχ . Moreover, since we are examining the equations of motion at order O(ε), it is only the first derivatives of the varying exponential that enter the source of the inhomogeneous part δX [1] . Effectively we can say that up to order ε we have The above implies that δX [1] will simply be the change of the background solution under δk I i = i εq i δc I . The same pattern will appear at all orders in the ε expansion, leading us to the further split of the solution to the inhomogeneous problem according to [17] In the end, the whole solution is determined by the variations δT , δµ and δc I and so does the horizon fluid velocity v i , the local chemical potential , and the vector potential δa (0) j . More specifically, we can identify (3.7) The above identification will prove useful in the next subsection where we give the effective theory of hydrodynamics that governs the fluctuations up to and including δX [1] in our expansion.

The effective theory
An important point of our construction is the way that we choose to impose the gravitational and Gauss constraints. As we discussed in section 2, these constraints should be imposed at once on any hypersurface at e.g. constant radial coordinate r. From the dual field theory point of view, the natural choice for this hypersurface would be near the conformal boundary as they become equivalent to the Ward identities of diffeomorphism, Weyl and global U (1) invariance. More specifically, the constraints L b = 0 and C = 0 give with F = dA the field strength of the external source one-form A a andȳ J s ,z I s are the sources for the complex scalar operators.
Contracting the stress tensor Ward identity with a vector Λ b gives The thermal gradient ζ and electric field E perturbations enter the boundary metric g ab and external field A a according to along with the source δz I s for the scalar field We are now going to make the choice Λ = ∂ t and perturbatively expand the contracted Ward identity to give the electric charge and heat conservation with δQ a = −δT a t − µ δJ a . Equations (3.12) define two conserved currents at the level of first order perturbation theory. From the point of view of the effective theory, this is a good starting point in order to give a closed system of equations, provided that we can express these currents in terms of the hydrodynamic variables δμ, δT and δĉ I [17]. However, we will see soon that, in the phases we are interested in, the current δJ a H which describes the transport of heat is different from δQ a . As we have argued in the previous subsection, the time derivatives scale according to ∂ t ∝ O(ε 2 ) while for the spatial derivatives we have ∂ i ∝ O(ε). This suggests that we need to consider the charge densities up to order O(ε) and the transport currents up to order O(ε 2 ) in (3.12).
We will write our theory in position space where all our functions will be denoted by hats. Moreover, from now on, we find it useful to define hatted thermodynamic quantities which are local functions of the dynamical temperature, chemical potential and phasons. For example, we defineŵ ≡ w(T + δT , µ + δμ, k I i + ∂ i δĉ I ). In appendix A.2 we relate the boundary currents δJ a , δQ a to the horizon currents we have defined in equation (2.32) in the ε-expansion. More specifically, we can write where we have defined the divergence free magnetisation currents At this point it is important to identify the correct current δJ i H that describes the transport of heat. For this reason, using the first law (2.22), we write the conservation equations as 7 The conservation equation (3.15) implies that up to magnetisation current contributions, the currents that correctly describe the transport of heat and electric charge are δĴ i H = δQ i + I w i I ∂ t δĉ I = δQ i (0) and δĴ i = δĴ i (0) . In terms of operators, we can writeĴ At first sight, it might seem surprising that the horizon heat current correctly describes the transport of heat. However, one might have expected this to happen since at the level of thermodynamics the entropy density of the system is determined by the horizon. This also ties well with the common lore in holography that dissipation is captured by horizon physics. The above discussion suggests that the physically relevant current to discuss iŝ J i H rather thanQ i . In the end, we would like to build our theory of hydrodynamics around the conserved currentsĴ i H andĴ i and the light modes associated to the operators S I . The corresponding triplet of sources is ζ i , E i , ξ S I for our choice of 7 It would be interesting to go beyond linear hydrodynamics and examine whether demanding positivity of entropy production (according to an appropriately defined entropy current) leads to constraints on the transport coefficients defined below, as happens generally [37,38], as well as in related contexts [21,26,[39][40][41].
operators and by examining the source terms in the deformed action we can easily show that Note that, as discussed in section 2.1, thermodynamically preferred phases satisfy w i I = 0, in which caseĴ i H =Q i . Given these results, we are able to express the boundary transport currents in terms of δT [1] , δµ [1] , δc I [0] and v i [2] at order O(ε 2 ). However, in [17] we have shown that within perturbation theory we can choose to impose the constraint L i = 0 on a hypersurface close to the horizon. This allows us to integrate out the horizon fluid velocity v i [2] and therefore obtain local expressions for the currents in terms of our hydrodynamic variables For the specific holographic model we are considering, the transport coefficients can be expressed in terms of horizon data and thermodynamic susceptibilities according to where we have used the definitions from appendix A.3, and indices in N are raised and lowered with the horizon metric g (0)ij . We now turn our attention to the pseudo gapless degrees of freedom δc I related to the density waves in our system. In order to introduce pinning to our system we turn on a perturbative deformation δφ I s ∝ O(ε 2 ) in the background asymptotics (2.15). The constitutive relations (3.18) remain unchanged [17], while in appendix A.1 we integrate the equation of motion (2.17) to obtain the effective Josephson relation with the transport coefficient  [21,24,25,42]. It would be interesting to find a holographic model which realizes that.
In the following subsections we study the gap of hydrodynamic modes after turning on the pinning parameters δφ I s as well as the dispersion relations without pinning. Finally we compute the finite frequency retarded Green's functions which will help us extract the transport properties of the system.

Pseudo-gapless modes
In order to identify the gapped modes in our effective theory we need to study perturbations with the sources switched off and with wavevector q i = 0. They belong to the Goldstone mode sector, which leads us to consider the ansatz δT = 0, δμ = 0, δĉ I = δc I 0 e −iδω g t .
where we have defined the matrix 8 (3.25) 8 We remind the reader that capital indices are not being summed over.
In order for equation (3.24) to have non-trivial solutions, the matrix multiplying the vector of amplitudes should not be invertible. Equating the determinant of this matrix to zero gives an algebraic equation which determines the N Z different values for the gaps δω g in terms of the eigenvalues of the matrix M −1 .
In order to illustrate the effect of the magnetic field on the gaps of the theory we consider a simple case with N Z = N Y = 2 and which apart from the U (1) 4 the model has a Z 2 × Z 2 symmetry which exchanges Z 1 ↔ Z 2 and Y 1 ↔ Y 2 . This allows us to consider an isotropic background which can be achieved by choosing The symmetries of the model along with the choice of boundary parameters leads to data of integration in which the internal indices can be suppressed, allowing us to write This simplifies the quantities where we see that the magnetic field introduces an antisymmetric piece in the matrix B ij . In the above ω c = ρB/s , (3.29) can be identified with the cyclotron mode frequency [1,2]. As a consequence, the matrix M will now have complex eigenvalues leading to the gaps for the pseudo-massless modes. It is easy to see that in zero magnetic field, these modes lie on the lower imaginary semi-axis and they agree with the expressions that were obtained in [17]. Moreover, due to the isotropy of the model and background we are considering, they also lie at the same point. However, the characteristic frequency of our nearly gapless modes has a resonant frequency apart from a gap at finite magnetic field.
It is interesting to consider the extreme limits for the behaviour of the poles (3.30) of our simple example. In the limit where B is the smallest parameter 9 we have a perturbatively small correction to the results of [17] (3.31) We therefore see that for small magnetic fields the two modes split and they move horizontally in opposite directions in the complex plane. In the opposite limit, where B is the largest scale in the system we have 10 32) and the two frequencies become degenerate once again, at the value which is given by the k → 0 limit of (3.30). The same result is obtained if we keep the filling fraction ρ/B finite while taking B → ∞. It is interesting to note that the expression (3.32) is the same with [16] but in a different thermal state. This situation is relevant to weak lattices and close to the phase transition at T ∼ T c . A third, distinct possibility which is relevant to weak lattices and magnetic fields is when σ 0 B ρ and B k 2 In this case we see that the resonant part dominates the poles related to the phase relaxation. The expression (3.33) agrees with the prediction from hydrodynamics in the pseudo-spontaneous regime once we identify the pinning frequency ω 2 0 ∼ k 2 Ω δφ s [23,24,43].
Note that the full expression (3.30) holds for finite magnetic field; in particular, we nowhere assumed that B is perturbatively small. All the background quantities however, including the horizon values Ψ (0) , Φ (0) , depend implicitly on B as well as k s . Thus, the scaling of δω g with B in the above expressions can be determined analytically (or, when not possible, numerically) from the properties of the ground state. 11 In particular, it is not straightforward to compare our results to [26], which found a B 1/2 scaling for large B using a holographic massive gravity model. 9 To be precise, the regime where (3.31) is valid is σ 0 B ρ and ω c k 2 Φ (0) . 10 In the regime σ 0 B ρ and B k 2 s Ψ (0) , k 2 Φ (0) . 11 The only exception is (3.31), since we expect the various background quantities to be continuously connected to the corresponding ones of the B = 0 state.

Finite frequency response
In this subsection we study linear response at zero wave number q i = 0. In order to achieve this we turn on all our sources with a finite frequency ω and look for a solution to the conservation law equations (3.15) and Josephson relation (3.21). A suitable ansatz for this purpose is δT = δT [1] e −iωt , δμ = δµ [1] After combining the transport heat current (3.18) and the Josephson relation (3.21), we obtain Plugging this solution back in the constitutive relations (3.15) we obtain the VEVs for the scalar fields S I and the transport currents in terms of the sources ζ i , E i and ξ S I according to In order to extract the retarded Green's functions, we need to consider the derivative of the VEVs with respect to the time dependent sources. Since we are only considering linear response, we can write 12 for any operator C in our theory. After a little algebra we obtain the expressions The non-trivial frequency dependence in the above thermoelectric conductivities includes only the horizon quantitiesᾱ ij H ,κ ij H , due to the fact that the Goldstone couples only to the heat current, (3.21). Moreover, from the above expressions we see that the gaps of the previous subsection determine the poles of the retarded Green's functions (3.38). From these expressions it is clear that the pseudo-gapless modes we discussed in section 3.3 couple to the conserved currents of our system. It should be emphasised that the Green's functions in equation (3.38) attain the most general form possible. In particular, G S I S K is simply controlled by the existence of a pole related to the pseudo-spontaneous symmetry breaking (subject to the presence of the gap), while S I couple to the currents only through their time derivatives, giving an additional factor of ω in the numerator. Thus, we expect to see the same structure in more general theories of this type.
Given the time reversal symmetry of the theory, as a non-trivial check, we can see that the expressions above satisfy the Onsager relations G CD (ω) , with ε C,D = ±1 depending on how the operators C and D transform under time reversal. In particular, we find that which hold by construction. The explicit expressions (3.38) show that the limits ω → 0 and M → 0 do not commute, as also discussed in [17,44]. Specifically, taking the DC limit ω → 0 while the gap remains finite, reduces the conductivities to the corresponding horizon conductivities, first obtained in [35], as well as the expected Goldstone mode susceptibility χ IK with the rest of the correlators vanishing. However, when we first take δφ I s → 0, transport at zero-frequency gets modified by the Goldstone modes while G S I S K diverges as (3.43)

Decoupling the Goldstone modes
Given the results of the previous subsection we can proceed to obtain Kubo formulas, which can be taken as the fundamental definition of the corresponding transport coefficients in a generic theory with the symmetry breaking pattern we are considering. 13 In purely spontaneous phases the susceptibilities diverge [45], but here this divergence is regulated by the perturbative explicit source δφ I s .
We first extract the transport coefficient Λ I K as which is a finite quantity, given by the combination of transport coefficients shown in (3.35) in our specific holographic model. We can then express γ i I and λ i I as The order of the limits is important, as was also noticed in [24]. We first need take the gap to zero in order to include the effects of the Goldstone modes in the low frequency regime, and then take ω → 0. Similarly we can definẽ since, as explained in section 3.2, only the heat current enters the Josephson relation. However, more generally, equations (3.45),(3.46) will hold in a generic theory in which we do not have explicit expressions for the low energy Green's functions, as long as the expressions (3.45), (3.46) remain finite as δφ I s → 0. It would be interesting to also define modified electric and heat current operators J i ,J i and J i H ,J i H which decouple from the Goldstone modes. This is satisfied as long as we demand that the Green's functions irrespective of the order of limits. Within the hydrodynamic regime, (3.45), (3.46) imply that indeed satisfy Note that these currents are related by Thus, from a holographic perspective, the combinations (3.49) isolate the horizon contribution to the electric and heat currents. As expected, the finite-frequency poles related to the pseudo-gapless modes cancel out in the above Green's functions, which turn out to be frequency-independent for low frequencies up to ω g . The Green's functions for the time-reversed currentsJ i ,J i H are simply the timereversed versions of (3.52) and so they also satisfy Onsager relations similar to (3.39). This can be seen by combining (3.51) and (3.40).
We can proceed further by recalling the relation (3.48) between the transport coefficients entering in (3.49). We then see that the combinations do not include contributions from the Goldstone modes, and can be solely expressed in terms of the original currents J i , J i H . As before, J i dec +J i dec is a well-defined vector operator.
For the retarded Green's function we find which turns out to be given by the horizon quantity σ ij 0 defined in (3.20). Similarly We thus observe that the horizon quantity σ ij 0 defined in (3.20) corresponds to the conductivity of the part of the U (1) current which decouples from the heat current J H and the Goldstone modes. In the special case of time-reversal invariant backgrounds with B = 0, we have that Then both decoupled combinations (3.53) reduce to the current considered in [46,47]. Furthermore, in the absense of a background lattice, the latter can be identified with the incoherent current which decouples from the conserved momentum operator [48].
Finally, note that all of the above results hold in the strong holographic lattice limit that we are considering in our paper, where the low frequency transport properties are determined by the Goldstone modes and the momentum non-conservation poles are outside our hydrodynamic regime. However, the explicit sources δφ I s also relax momentum apart from the massless modes with a relaxation rate ∼ (δφ I s ) 2 . So, had we not included such a background lattice, the momentum poles would dominate over the Goldstone mode poles in the hydrodynamic regime.

Hydrodynamic modes
In this subsection we wish to extract the dispersion relations of the hydrodynamic modes in our system at zero pinning. Similarly to the previous section, we switch off all the sources and we look for solutions of the form which solve the conservation law equations (3.15) and Josephson relation (3.21).
Similarly to the previous subsection, the resulting system of equations reduces a linear system of equations for the vector of amplitudeŝ where we have defined the matrix which further implies that the dispersion relations are independent of the sign of B, since the determinant is invariant under transposition. However, the Kernel ofŜ which solves equation (3.57) will depend on it and thus the actual modes will change as we change the sign of B. By directly exploiting these properties, we show below that the frequencies ω(q i ) of our hydrodynamic modes are pure imaginary. Finally, we notice that the transformation q i → λ q i , ω → λ 2 ω and δc I → λ −1 δc I is a symmetry of equation (3.57). This shows that all our modes are diffusion-like with ω(λ q i ) = λ 2 ω(q i ).
In contrast, the effective field theory of [24], as well as the holographic model of [26], both include a real quadratic part in the dispersion relation for the magnetophonon mode in the purely spontaneous or pseudo-spontaneous symmetry breaking regime. In our case, the modes are purely diffusive as a result of being in the strong translational symmetry breaking regime.
Let us now prove that the modes are pure imaginary. We find convenient to split the mode vector which solves (3.57) according to (3.61) Then for the background with B → −B there is a different mode with |ṽ 1 and |ṽ 2 but with the same dispersion relation ω(q i ). This can be justified by using the transformation property (3.60) and the comment below it. We can write From the above systems and after using (3.59) we obtain the relation showing that iω has to be a real number. For B = 0, the matrices X, W, Σ and Θ are symmetric. We also know that the vectors |v 1 and |v 2 coincide with the vectors |ṽ 1 and |ṽ 2 . This observation shows that if the matrices X and W are positive definite, then iω > 0. In other words, at zero magnetic field, thermodynamic stability implies dynamical stability in the hydrodynamic sector that we focussed on.

Numerical checks
In this section we numerical confirm the results presented in section 3, and in particular the formula for the dispersion relations of the hydrodynamic modes coming from equation (3.58), the gap (3.30) and the optical conductivities (3.38). This is achieved by focusing on the model of [44], which is a truncation of the general bulk action (2.1) down to the four-dimensional Einstein-Maxwell theory coupled to six real scalars, φ, ψ, χ i and σ i with i = 1, 2 The variation of the above action gives rise to the following field equations of motion The simplest solution to the above equations is the unit radius vacuum AdS 4 , which is dual to a d = 3 CFT with a conserved U (1) charge. In this work we choose to place the CFT at finite temperature and deform it by a chemical potential, an external magnetic field and a background lattice. Within this theory, we are interested in thermal states that correspond to density waves. Putting all the ingredients together, the solutions we are after are captured by the ansatz (2.14), which we rewrite here for convenience where I = 1, 2, i = 1, 2. For simplicity we choose k I i = k i δ Ii , k I si = k si δ Ii . Let us now move on to discuss the boundary conditions. In the IR, we demand the presence of a regular Killing horizon at r = 0 by imposing the following expansion U (r) = 4π T r + . . . , which is specified in terms of 6 constants. In the UV, we demand the conformal boundary expansion U → (r + R) 2 + · · · + W (r + R) −1 + . . . , V 1 → log(r + R) + · · · + W p (r + R) −3 + . . . , Just like in [44], the scalar fields (ψ, σ) are taken to constitute the anisotropic Q-lattice in which both translational invariance and the two U (1) ψ symmetries are explicitly broken, while the density wave phase is supported by (φ, χ) and breaks the the two U (1) φ symmetries spontaneously. As such, the thermal states of interest correspond to taking ψ s = 0 and φ s = 0. Thus, this expansion is parametrised by 8 constants.
Overall we have 14 constants appearing in the expansions, in comparison to the 11 integration constants of the problem. Thus, for fixed γ, δ, B, k i , k si and temperatures below a critical one T < T c , we expect to find a 3 parameter family of solutions, labelled by ψ s , µ, T .
In figure 1 we plot the critical temperature, T c , as a function of k = k 1 = k 2 for a particular choice of parameters. This is obtained by considering linearised fluctuations around the normal phase of the system (φ = 0, χ = 0) and exhibits the usual "Bell Curve" shape.  Figure 1: Plot of the critical temperature at which the background Q-lattice becomes unstable as a function of k for (k s1 , k s2 , ψ s , γ, δ, µ, B) = ( 3 10 , 3 10 , 4, 3, 1, 1, 1 10 ). We see that the most unstable mode corresponds to k = 0.

Quasinormal modes
We now move on to compute quasinormal modes for the backgrounds constructed above. For simplicity, we focus only on isotropic backgrounds characterised by k 1 = k 2 ≡ k, k s1 = k s2 ≡ k s , V 1 = V 2 . We consider perturbations of the form together with (δa t , δa 1 , δa 2 , δφ, δψ, δχ 1 , δχ 2 , δσ 1 , δσ 2 ), where the variations are taken to have the form δf (t, r, x 1 ) = e −iωv(t,r)+iqx 1 δf (r) , (4.8) with v EF the Eddington-Finkelstein coordinate defined as v EF (t, r, Compared to the analytic setup of the problem in section 3, we have chosen S in (2.29) such that S = U −1 , as well as a radial gauge in which all perturbations with an r index vanish. Such a gauge is not compatible with the way we constructed the modes in section 3, but the physical information of the quasinormal modes in the end should of course be the same. Note also that our choice for the momentum q i to point in the direction x 1 is without loss of generality, because the background is isotropic. Plugging this ansatz in the equations of motion, we obtain 5 first order ODEs and 10 second order giving rise to 25 integration constants. Let us now discuss the boundary conditions that we need to impose on these fields. In the IR, we impose infalling boundary conditions at the horizon, which without oss of generality is set at r = 0 δh tt = c 1 r + . . . , δh t x 2 = c 3 + . . . , where the constants c 1 , c 2 , c 3 and c 6 are not free but are fixed in terms of the others. Thus, for fixed value of q, the expansion is fixed in terms of 11 constants, ω, c 4 , c 5 , c 7 , c 8 , c 9 , c 10 , c 11 , c 12 , c 13 , c 14 .
On the other hand, in the UV, the most general expansion with φ s = 0 is given by For the computation of quasinormal modes, we need to ensure that we remove all the sources from the UV expansion up to a combination of coordinate reparametrisations and gauge transformations [δg µν + Lζg µν ] → 0 , where the gauge transformations are of the form for ζ µ , λ constants. This requirement demands that the sources appearing in (4.11) take the form δh (s) where λ 2 = Bζ 3 . Therefore, the UV expansion is fixed in terms of 15 constants: ζ 1 , ζ 2 , ζ 3 , ζ 4 , λ and δh (v) 2 . Overall, for fixed q, B, we have 26 undetermined constants, of which one can be set to unity because of the linearity of the equations. This matches precisely the 25 integration constants of the problem and thus we expect our solutions to be labelled by q and B. We proceed to solve numerically this system of equations subject to the above boundary conditions using a double-sided shooting method. Figure 2 (left) shows the dispersion relations for the four hydrodynamic quasinormal modes in our system for a particular choice of the background configuration. We also illustrate with dashed lines the dispersion relations 14 fixed by the linear system 3.57. Figure 2 (right) shows the analytically predicted diffusion constants from 3.57, and the ones computed numerically from the q → 0 limit of the function iω(q) /2. 15 Let us now make some remarks on figure 2. First of all, we note that all the modes we find are diffusive and purely imaginary, as expected from the analysis in section 3.6. We also see a good quantitative agreement of the numerical solution and the analytical expressions in the regime of validity of hydrodynamics, where q is parametrically smaller than all the dimensionful scales in the system. Actually, we expect the radius of convergence of hydrodynamics to be set by the collision points of the hydrodynamic modes with the first non-hydrodynamic mode [49][50][51]. In figure 2 we chose the parameters such that the lattice is weak; in this case, the lowest lying non-hydrodynamic mode is the momentum relaxation/cyclotron mode. One of the thermoelectric modes, the steepest curve in figure 2, interacts with this non-hydrodynamic mode as q is increased, leading to the quickest deviation from the analytic quadratic dispersion relations. The top curve describes the incoherent thermoelectric mode which decouples from momentum, and agrees very well with the analytic expression even for very large q. 16 This is similar to the case without magnetic field, see [17] for further details.
14 Note that, in order to evaluate the quasinormal modes using the analytic formula (3.58), we need to compute the derivatives w i I and w ij IJ . In order to compute these correctly one needs to consider backgrounds with general k I i , i.e. k 1 1 = k 1 2 = k 2 1 = k 2 2 . 15 Note that, as explained below (4.9), we have chosen the momentum q to point in x 1 , and thus, in this setting, the diffusion matrices D ij become diffusion constants D for each mode. 16 The above characterization of the modes as thermoelectric versus Goldstone is done by examining the system as k → 0; in general all the modes are coupled.

Pseudo-gapless modes and two-point fucntions
In this subsection we outline the numerical computation of the pseudo-gapless modes as well as certain two-point functions involving the currents J, Q in the presence of pinning, φ s = 0. We perform a calculation similar to the one for quasinormal modes, but we now consider fluctuations with q = 0 around a background configuration that has a small but finite source φ s = 0. Looking at the ansatz (4.7), it is consistent to set δh tt , δh x 1 x 1 , δh x 1 x 2 , δh x 2 x 2 , δa t , δφ, δψ = 0. We are thus left we 6 second order and 2 first order equations for the remaining fluctuations, giving rise to 14 integration constants. The IR expansion close to the horizon (r = 0) takes a similar form as above, namely δh t x 2 = c 3 + . . . , δa x 1 = c 7 + . . . , δa x 2 = c 8 + . . . , δχ 1 = c 11 + . . . , δσ 1 = c 12 + . . . , δχ 2 = c 13 + . . . , δσ 2 = c 14 + . . . , (4.15) where the constants c 2 , c 3 are fixed in terms of the others. We see that the expansion is fixed in terms of 7 constants, ω, c 7 , c 8 , c 11 , c 12 , c 13 , c 14 .
On the other hand, the UV expansion changes slightly in comparison to (4.11) because φ s = 0. In particular, it is given by Once again, we remove all the sources from the UV expansion apart from an external electric field E and a temperature gradient ζ in the x 1 direction, up to a combination of coordinate reparametrisations and gauge transformations. This is done by imposing the following constraints on the sources in (4.16) δh (s) Let us first consider the case of the pseudo-gapless modes by setting (E, ζ) = (0, 0). We see that the UV expansion is fixed in terms of 8 constants: 2 . Overall, we have 15 undetermined constants, one of which can be set to unity because of the linearity of the equations. This matches precisely the 14 integration constants of the problem and thus we expect to find a discrete set of solutions, labelled by B. We proceed to solve numerically this system of equations subject to the above boundary conditions using a double-sided shooting method aiming to identify the two pseudo-gapless modes of equations (3.30). Note that the two modes have equal imaginary parts and opposite real parts. In figure 3 we plot the real and imaginary part of these modes as a function of the pinning parameter, φ s , and the external magnetic field, B, and we compare with the analytic formulas which are depicted with dashed lines. We see that the numerical and analytic calculations are in good agreement. The reader is reminded that the analytic computation is perturbative in φ s , but exact in B.
We finally consider the computation of the conductivities. From (4.17) we see that, for fixed (E, ζ), we have 7 constants in the IR and 8 in the UV. Comparing with the 14 integration constants in the problem, we expect to find a 1-parameter family of solutions labelled by ω. Using the linearity of the equations we set (E, ζ) = (1, 0) or (E, ζ) = (0, 1) depending on which source we want to keep. The diffusion currents   are then given by Carrying out the numerical shooting computation, we calculate the (1, 1) and (1, 2) components of the two-point functions For fixed B and φ s , these quantities are plotted in figure 4 with solid lines. In order to compare our numerics with the analytic results of section 3, we use the definition (3.16) to write and thus obtain analytic expressions using (3.38), which are depicted in figure 4 with dashed lines. We see that the two are in good quantitative agreement at small frequencies.
The reader is reminded that in this calculation we have set ζ I S = 0 and we only included sources in the x 1 direction; thus we can not compute the (2, 1) and (2, 2) components of the the two-point functions.

Discussion
In this paper we constructed the effective theory of hydrodynamics which captures holographic phases in which translations are broken explicitly and spontaneously. We have significantly extended the construction of [17] to include an arbitrary number N Z of gapless degrees of freedom emerging from spontaneous density waves and we also included a background magnetic field.
A holographic model which incorporates the two Goldstone modes arising from spontaneous breaking of translations in magnetic fields, along with the coupling to the heat current was studied in [26] and a complex quadratic dispersion relation was found. In our setup, the strength of the explicit breaking is large compared to the wavelength of the hydrodynamic fluctuations. In section 3.6 we analytically derived an equation whose roots yield the dispersion relations of the hydrodynamic modes governing our system. Despite not being able to write down the dispersion relations of all of our 2 + N Z hydrodynamic modes in closed form, we prove that they are purely imaginary and diffusive, unlike [26].
In our construction we have also included the corresponding N Z perturbative deformation parameters which pin down the density waves and introduce N Z gaps in our theory. Interestingly, we have shown that apart from the gap, the magnetic field field causes the corresponding poles to move off the imaginary axis due to resonance effects. In section 3.4 we computed the retarded Green's functions of the operators relevant to the hydrodynamic description of the system. As one might expect, the poles due to pinning have a direct effect on the transport properties of our system as can be seen from the explicit form of the Green's functions in equation (3.38).
Finally, an important byproduct in our work is the identification of the correct current in (3.16) which describes the transfer of entropy as can be seen by the conservation equation in the last line of (3.15). Given this definition, the variation of the free energy w i I with respect to the wavenumbers k I i drops out of the corresponding Green's functions (3.38). Moreover, the gaps and the resonance frequencies which can be found by solving the eigenvalue problem equation (3.24) are also independent of w i I .
There are various open questions which one could further explore. It would be interesting to consider second order hydrodynamic perturbation theory and examine what the second law implies for the transport coefficients in phases with spontaneous and explicit symmetry breaking. Additionally, it is important to examine how transport in such phases is constrained from purely field theoric considerations, such as the Ward identities, and also investigate the possible experimental significance of the decoupled/incoherent currents we defined in this paper. Finally, it would be enlightning to move away from homogeneity and explore what kind of novel effects inhomogeneous models with similar symmetry breaking patterns might exhibit. in our ε expansion (3.1), which we carry out in Appendix A.1. Then, in Appendix A.2 we derive the constitutive relations for the currents by relating the field theory and horizon currents densities of our holographic model.

A.1 Perturbations for χ I
After perturbatively expanding the equation of motion (2.17), we have From the form of the solution close to the conformal boundary at r → ∞, we can infer a relation between the sources ζ S I of the operators S I and their vevs δ S I . This will essentially give a Josephson type of equation for the variable δĉ I through equation (3.13). In the next subsections we will solve equation (A.1) in an ε expansion.

A.1.1 Field theory interpretation at order O(ε)
At order O(ε) we obtain the equation We integrate this equation for δχ I [1] while insisting on the near horizon behaviour (2.30). After doing so, we obtain the asymptotic behavior [1] δc I Demanding that the operator S I is not sourced at this order, we must have

A.2 Constitutive relations for the thermoelectric currents
In this subsection we will relate the horizon current densities (2.32) to the boundary quantities δJ i and δQ i that appear in the current conservation equation (3.12).
The bulk electric current is defined as The equations of motion (2.2) imply Following [52], for any vector Λ µ in the bulk we can define the bulk two-form where Λ µ F µν = ∂ ν f + β ν , with β a 1-form and f a globally defined function. After using the equations of motion (2.2), its divergence can be brought to the form We now consider Λ µ = ∂ t , and a general perturbation around the background ansatz (2.14) (not necessarily of the form (2.28)). The bulk heat current is defined as δQ i bulk = √ −g G ir = U 2 √ g g ij ∂ r δg jt U − ∂ j δg rt U − a t δJ i bulk = U 1/2 √ g 2K i t + U 1/2 g ij ∂ t δg rj − a t δJ i bulk , (A. 12) where we have used the result of Appendix B of [36] for the extrinsic curvature component Writingt µ ν = −2K µ ν + Xδ µ ν + Y µ ν , where X = 2K + · · · and Y are additonal terms that come from the counterterms, we recognizet µ ν as the field theory stress tensor, when evaluated on the boundary. Evaluating (A.12) at the boundary, this gives where t µ ν = r 5tµ ν . Note that the contribution from Y i t , as coming from (2.20), and contribution from the term involving a time derivative are subleading even in the precense of sources. This result matches the expression for the boundary heat current obtained from the variation of the action in the presence of the sources as in (2.28) where h µν = g µν −n µ n ν and n is the unit norm normal vector. Furthermore, equation (A.11) implies the radial dependence
We now solve the horizon vector constraint at order O(ε 2 ) to write