Intertwined orders in holography: pair and charge density waves

Building on [1], we examine a holographic model in which a U(1) symmetry and translational invariance are broken spontaneously at the same time. The symmetry breaking is realized through the Stückelberg mechanism, and leads to a scalar condensate and a charge density which are spatially modulated and exhibit unidirectional stripe order. Depending on the choice of parameters, the oscillations of the scalar condensate can average out to zero, with a frequency which is half of that of the charge density. In this case the system realizes some of the key features of pair density wave order. The model also admits a phase with co-existing superconducting and charge density wave orders, in which the scalar condensate has a uniform component. In our construction the various orders are intertwined with each other and have a common origin. The fully backreacted geometry is computed numerically, including for the case in which the theory contains axions. The latter can be added to explicitly break translational symmetry and mimic lattice-type effects.


Introduction
Behind the unconventional behavior of many strongly interacting quantum systems is an intrinsically complex phase diagram exhibiting a variety of orders, some of which exist at comparable temperatures.A prime example is that of the high T c superconductors, whose peculiar properties seem to be linked to the interplay between different phases, which may not only compete but also cooperate with each other (see e.g.[2,3,4]).Indeed, in certain regions of the phase diagram many of these orders are believed to have a common origin and to be intertwined.
Our interest here is in a very particular quantum phase of matter [5,6], the so-called pair density wave (PDW), which describes a state in which charge density wave (CDW) and superconducting (SC) orders -and possibly spin density wave (SDW) order -are intertwined in a very specific manner.Thus, these are peculiar types of striped superconductors.The most compelling experimental evidence of a PDW is found in the cuprate superconductor La 2−x Ba x CuO 4 (LBCO), but there is also computational evidence suggesting that it might be a robust feature of strongly correlated electron systems [4].The defining properties of PDW order that we will focus on are the following: • The superconducting order itself is spatially modulated in such a way that its uniform component is zero.
• The charge density oscillates at twice the frequency of the condensate.
Thus, a PDW differs from a state with co-existing SC and CDW orders, in which the superconducting condensate does have a homogeneous component.As argued in [4], one way to understand highly complex phase diagrams such as those of the cuprates may be to view the many observed phases as having originated from a "parent phase" which spontaneously breaks a number of symmetries.This is also tied to the idea that certain orders may not compete, but rather cooperate with each other.In this paper we explore these ideas within the framework of holography, by providing an example of intertwined orders which realizes the two PDW features described above.In particular, we extend the work we started in [1] and examine a four dimensional holographic model in which translational invariance and a U(1) symmetry are broken spontaneously and at the same time.As a result, in the dual theory one finds a scalar condensate that is spatially modulated, and similarly an oscillating charge density, exhibiting unidirectional stripe order.Thus, this is a model of a striped superconductor 1 .
Depending on the choice of parameters, our model exhibits either PDW order or co-existing SC+CDW orders.In the former case the oscillations of the condensate average out to zero, O χ ∝ cos(k x) , (1.1) and those of the charge density are induced (in a manner which will be described in Section 4) and oscillate at twice the frequency, ρ(x) = ρ 0 + ρ 1 cos(2k x) . (1.2) On the other hand, in the phase with co-existing SC+CDW orders the scalar condensate has a uniform component, and now oscillates at the same frequency as the CDW.As we will discuss shortly, the model also contains a second "charge" density (which could in principle mimic a spin density) which also oscillates at the frequency of the condensate.
The setup we consider is not of the form of the standard holographic superconductor [7,8], but instead falls within the generalized class of theories examined in [9], where a complex scalar is replaced by two real scalar degrees of freedom, which have the advantage of admitting less restrictive scalar couplings.Indeed, in our construction the spontaneous symmetry breaking of the global U(1) symmetry of the dual field theory occurs via the Stückelberg mechanism2 , and the two real scalar fields χ and θ are not necessarily associated with the magnitude and phase of a complex scalar.While it should be possible to construct models in which the PDW phase is associated with the condensation of a complex scalar -more appropriate for applications to e.g.cuprate high T c superconductors -this would likely entail using a more involved matter sector, and in this paper as a first step we restrict ourselves to this particular framework.We will come back to this issue in Section 7.
Our theory involves gravity coupled to two real scalars, χ and θ, and two vector fields A µ and B µ , with A µ playing the role of the gauge field which provides the charge density of the dual field theory.The high temperature phase is described on the gravity side by the standard electrically charged AdS-Reissner Nordström black brane supported by A µ , with the scalars and the second vector B µ in the vacuum.At sufficiently low temperatures, on the other hand, χ and B µ spontaneously develop non-trivial profiles and become dynamical 3 .When this happens the U(1) symmetry of the high temperature phase is broken spontaneously and a scalar condensate develops in the dual field theory.By allowing for perturbations which are spatially modulated, the instability can be seen to lead to a condensate and charge density which are striped.In particular, depending on whether the scalar χ is charged under the second vector field B µ or not, we find PDW vs. SC+CDW order.The physical motivation of the second vector field in this setup is not obvious.It can be thought of as playing a somewhat passive role, imprinting its own oscillations onto the charge density associated with A µ , i.e. giving rise to the properties of the CDW.On the other hand, the oscillations of its density can in principle be associated with those of a spin density, i.e. generating SDW order in the system.This interpretation, however, requires an explanation for the specific couplings we have chosen in the setup, which we discuss in the next section.At this stage we simply point out that it is possible to obtain some of the key oscillatory features of intertwined SC, CDW and SDW orders, but more work is clearly needed to realize a model which would capture the detailed phenomenology of these highly complex systems.
Holographic superconductors with spatially modulated orders have been studied in a number of setups, starting with [18], in which striped modulations were sourced by a periodic chemical potential, with many generalizations in the literature since then, e.g.[19,20,21,22,23].A study of backreaction by adding a lattice in the form of a periodic potential was reported in [24].Note that in these setups the translational invariance was broken explicitly by introducing some kind of spatially modulated source.The competition between superfluid and striped phases has also been examined in holography, see [25,26] for top down models.The spontaneous formation of striped order in a holographic model with a scalar coupled to two U(1) gauge fields was first studied in [27] and more recently in [28,29,30,31] (but the U(1) symmetry was preserved in these models).Here we have extended such constructions by spontaneously breaking both symmetries at the same time, and focusing on the differences between PDW order and the better studied CDW order in holographic superconductors 4 .What has recently become apparent in holographic models of strongly correlated electrons is the advantage of using multiple vector fields, as they typically lead to richer physics.Our construction provides a further example of this idea.Note that while in our analysis the mass of the vector B µ is not expected to affect the phase structure in a qualitative way, it would play a role for applications to transport.
In this paper we have extended [1] in two main ways.First of all, we have performed an analysis of the fully backreacted geometry, done numerically using the DeTurck method [34].We have also included the perturbation equations for the system, at next-to-leading order in fluctuations, which provide intuition for the numerical results.Second, we have examined the effect of including axions into the model, introduced to mimic the effects of an underlying lattice structure and to provide a mechanism for momentum relaxation.
The outline of the paper is as follows.We introduce our model in Section 2 and examine the critical temperature for the onset of spatially modulated instabilities in Section 3. The difference between PDW and SC+CDW orders is explained in some detail in Section 4, and the numerical analysis of the backreacted striped solutions is presented in Section 5 (no axions) and Section 6 (axions).We conclude with open questions and remarks in Section 7. The full perturbation analysis up to nextto-leading order is included in Appendix A, while Appendices B and C contain, respectively, a discussion of the free energy and a check on the quality of our numerics.

Holographic Setup
We work with a four-dimensional gravitational model which describes two real scalar fields χ and θ coupled to two abelian vector fields A µ and B µ , where This model falls within the generalized class of holographic Stückelberg superconductors5 studied in [9], in which the spontaneous breaking of the global U (1) symmetry happens via the Stückelberg mechanism and the local gauge invariance is encoded in the transformations with an analogous transformation for the second vector field B µ when it is massless.Within the framework of the Stückelberg mechanism χ and θ are not necessarily the magnitude and phase of a complex scalar field, and in particular χ does not have to be positive.As we will see, this will be important in order to realize the main features of PDW order in this particular model.It has implications for the nature of the U(1) symmetry, which in our case turns out to be non-compact, as there are no restrictions on the range of θ.We will come back to this point in Section 7.For our purposes the advantage of the Stückelberg model (2.1) is that it is less constraining than the standard holographic superconductor [7,8], and allows for the type of couplings that are needed to realize the striped orders we are after.One should keep in mind that there may be other ways to engineer similar phases, using a different matter content.
For completeness here we take the scalar χ to be generically charged under both vectors.However, as we will see in detail the case with q B = 0 and q A = 0 will enable us to describe PDW order.Our focus will be on breaking the global U(1) symmetry associated with A µ , the gauge field responsible for the finite charge density in the boundary theory.Thus, here we will always keep q A = 0.Moreover, note that while the current dual to A µ is always conserved, this is not the case for the second vector field B µ when m 2 v = 0. We will consider separately the cases in which B µ is massless or massive.We allow for a mass term because it is expected to play a key role in the transport properties of the system, which we plan to examine in future work, even though it does not affect in any qualitative way the analysis done in this paper.
We have taken all the gauge fields' couplings to depend on the real scalar χ, and have allowed for an interaction ∼ Z AB (χ) between the two fluxes, which, as we will see, will be responsible for the breaking of translational invariance.Finally, the cosmological constant is given by Λ = − 3 L 2 , with L denoting the AdS radius.Einstein's equations resulting from the action above read (2. 3) The equations of motion for the matter fields on the other hand are where we used = ∂ χ , F 2 = F µν F µν , F 2 = Fµν F µν and F F = F µν F µν for short.The scalar couplings and the potential will be chosen so that in the limit χ → 0 they can be expanded in the following way, with (a, b, c, m, κ) constants.The motivation behind these particular choices will become apparent when we discuss the symmetry breaking mechanism.We want to work with a theory whose dual has a finite charge density with respect to the global U(1) symmetry associated with the gauge field A µ .One consequence one shall see immediately is that the expansion (2.8) allows for the standard electrically charged AdS Reissner-Nordström (AdS-RN) black brane supported by A µ , in which the two scalars and the second gauge field are in the vacuum, (2.9) Here r h denotes the horizon, µ the chemical potential and ρ the charge density.The associated temperature is and in the extremal limit T = 0 the horizon radius can be easily seen to be 3 .It is well-known that the extremal, near horizon geometry is AdS 2 × R 2 , with r = r − r h and the horizon radius given by L (2) = L/ √ 6.The AdS-RN solution (2.9) describes the high temperature phase in which the dual theory has a global U(1) symmetry.However, we are interested in backgrounds in which χ and B µ are dynamical and K(χ) = 0, which will describe phases with a broken U(1) symmetry (θ can be gauge-fixed to zero).
Indeed, in our model at sufficiently low temperatures the scalar χ and the second vector field B µ spontaneously develop non-trivial profiles.When this happens, and in particular when K(χ) no longer vanishes, the U(1) symmetry of the high temperature phase (associated with the A µ vector field) is broken spontaneously, and a scalar condensate forms in the dual field theory.Moreover, by considering perturbations which source spatial modulations, one can trigger instabilities to striped superconducting phases.Indeed, we showed in [1] that at finite charge density our model can exhibit spatially modulated phases which spontaneously break translational invariance at the same time as they break the U(1) symmetry.The detailed features of the model will be particularly sensitive to the value of q B , as we discuss next.

Spatially Modulated Instabilities
Here we come back to the analysis done in [1], where we examined the spatially modulated static mode in the spectrum of fluctuations around the unbroken phase (2.9).We are going to first consider instabilities of the zero temperature domain wall solutions which interpolate between AdS 4 and an electrically charged AdS 2 × R 2 solution in the far IR.In particular, we are going to construct PDW or CDW-type modes which violate the IR AdS 2 BF bound.As is well known, the presence of such zero temperature modes is a strong indication that analogous instabilities should exist at finite temperature, and that there should be a region in which one has spatially modulated superconducting order.After examining analytically the perturbations of the zero temperature IR background, we repeat the analysis at finite temperature and focus on the behavior of the critical temperature for the onset of the phase transition, as a function of wave number k.

Instabilities of the
To investigate striped instabilities of the electrically charged AdS-RN black blane (2.9), we start by considering linearized perturbations in the AdS 2 × R 2 background (2.11) describing its extremal IR limit, where we have relabeled r → r for convenience.We turn on the following spatially modulated perturbations, where we introduce ε as a formal expansion parameter which we will use to make sure that the perturbations are indeed small when compared to the background (3.1).
Working at linear level O(ε) in perturbations, we then obtain the two coupled equations where (a, c, κ, m) have been defined in equation (2.8).Note that they are only coupled to each other when c = 0 (recall that c controls the coupling Z AB between the scalar and the two vectors) and that b does not appear in the linearized equations.We consider the following mode where u 1 , u 2 are constants and λ is related to the scaling dimensions of the operators in the one-dimensional CFT dual to the AdS 2 geometry.The linearized equations can then be written in matrix form where we have defined To have a non-trivial solution for u 1 and u 2 , the determinant of the matrix on the left hand of (3.6) should be vanishing, from which we can solve for λ, with (3.9) Notice that we have fixed the chemical potential and taken it to be µ = 1.
The onset of the instability associated with the violation of the AdS 2 BF bound occurs when λ becomes imaginary, i.e. when Thus, in order to have a striped instability, one must identify a non-zero wave number k at which the value of λ is imaginary, for a fixed choice of Lagrangian parameters.Such BF-bound violating modes are then associated with spatially modulated phases.
It is easy to check that by choosing parameters appropriately, the BF bound can be violated.For an explicit example, note that when 2 and c = 3, we find a mass which is clearly below the AdS 2 BF bound for a certain range of k.

Finite temperature instabilities
The zero temperature instability analysis of the IR AdS 2 region indicates that similar spatially modulated unstable modes should be present at finite temperature, in the black brane background (2.9).Indeed, we are now ready to calculate the critical temperature below which the AdS-RN geometry becomes unstable, as a function of wave number k.When the scalar field instabilities are associated with a non-zero value of k, the condensate is spatially modulated.
In analogy with the AdS 2 analysis, we turn on the fluctuations given in (3.2).To identify the critical temperature for the instability, it is again sufficient to work to linear order in perturbations.After expanding around the AdS-RN background (2.9) we obtain the two coupled ODEs, with the coupling between the two modes once again proportional to c.We require the fluctuations to be regular at the horizon r = r h , with the expansion (3.12) The UV expansion (as r → ∞) is instead given by where ∆ χ and ∆ B are the scaling dimensions of the operators dual to χ and B µ , respectively 6 .The sources for the operators, encoded by the parameters w s and b s , are going to be turned off, since we want to break the U(1) symmetry and translational invariance spontaneously and not explicitly.
After fixing the parameters in the model, for a given wave number k there will be a normalizable zero mode appearing at a particular temperature, which is therefore the critical temperature associated with the phase transition.We are going to take the mass of the scalar field χ to be m 2 L 2 = −2 throughout this section, so that the scaling dimension of the dual operator is ∆ χ = 2.We have considered two different cases for the vector field B µ .First, we have taken it to be massless, corresponding to a conserved current in the dual field theory, with a UV expansion given by Next, we have allowed for a mass term m 2 v L 2 = 0.11, in which case the current dual to B µ is no longer conserved and the expansion is instead given by The associated critical temperatures are shown in Figure 1, and correspond to having set w s = b s = 0.The figure exhibits the typical bell curve behavior as a function of wave number -the curves are peaked at non-zero values of k, indicating that the condensate is indeed driven by the momentum-dependence of the spatially modulated instabilities.We emphasize that the mass of the vector does not play any significant role in the behavior of T c -by examining more explicit examples one finds that the two cases do not exhibit any qualitative differences.A more important role is played by the strength c of the coupling Z AB ∼ cχ between the two vector fields, as expected, since the latter is responsible for couplings the modes to each other.As seen in Figure 2, the finite-momentum peak becomes less and less pronounced (and shifts towards zero) as |c| decreases.What this is indicating is that, while for small |c| one may still have a superconducting instability, it is not going to be striped.As a result, in this model the coupling must be non-zero and sufficiently large, in order for the phase transition to occur at a finite value of k.This singles out the coupling c between the scalar and the two vector fields as the leading parameter controlling the onset of the striped phase transition.However, note that when |c| becomes too large the instability once again disappears, because the BF bound can no longer be violated.

PDW vs. SC+CDW
By performing a perturbative analysis and working to leading order in perturbation theory, we have just found the static, normalizable zero modes in the electrically charged AdS-RN background, and obtained the critical temperature at which the geometry becomes unstable to the formation of a spatially modulated phase.From the bell-type behavior of Figures 1 and 2 we see that there is a particular wave number k c associated with the critical temperature T c at which a new branch of black branes appears spontaneously.The wave number k c thus characterizes the breaking of translational symmetry.
As we have shown numerically and discuss in detail in the next section, for temperatures just below T c the scalar operator O χ dual to χ acquires spontaneously a spatially modulated expectation value, breaking the (non-compact) global U(1) symmetry of the boundary theory.Therefore, this spatially modulated phase is always associated with a non-vanishing condensate, i.e the new phase is a type of striped superconductor.Moreover, the "charge" density ρ B dual to B µ becomes striped, and this, in conjunction with O χ , induces a modulation in the charge density ρ A associated with A µ .The crucial feature in this model is that the frequency of oscillations of ρ A as compared to that of O χ , which is one of the determining factors for whether the state is a PDW or a SC+CDW, depends on whether the charge q B vanishes or not 7 .In particular, when q B = 0 we find that the time component of the vector operator J t A corresponding to A t becomes spatially modulated with which implies that the associated period of oscillations is precisely one half of that of the scalar condensate given in (4.1).Meanwhile, the time component J t B of the operator corresponding to B t oscillates as in with the same frequency as the scalar condensate.In addition, when q B = 0 the scalar condensate does not have a homogeneous component, so that its oscillations average precisely to zero.The behavior of ρ A and O χ match precisely with those of a PDW phase, as discussed in the introduction.The situation is different when q B = 0.In this case the charge density ρ A associated with the A µ vector field oscillates at the same frequency as the scalar condensate, and the latter develops a homogeneous component.Thus, the q B = 0 case describes a phase with coexisting SC and CDW orders.Notably, while the second vector field B µ does not determine the type of order (PDW or SC+CDW) developed in the system, it can in principle be associated with spin degrees of freedom, with its modulated density ρ B describing SDW order.Recall that if we set m 2 v = 0, the dual theory has a second conserved current corresponding to B µ .In the numerical analysis of the next section, for numerical convenience we will focus precisely on this massless case, and leave the massive case to future work (we anticipate that it may be useful when discussing the transport behavior).
Our numerical results can be explained in part by inspecting the structure of the perturbations at next-to-leading order.In particular, to quadratic order in ε, the perturbative parameter which is proportional to 1 − T /T c , we have where we are singling out the scalar and vector fields for the sake of space.We refer the reader to Appendix A for the full perturbation analysis.We find that the order O(ε 2 ) components of δχ and δB t are sourced by O(ε) terms with a prefactor proportional to q A q B , and therefore vanish when q B = 0.In particular, note that this implies that the homogenous perturbations χ (1) (r) and b (1) t (r), which fall under Set (I) of Appendix A, can both consistently be chosen to vanish when q B = 0.This causes the scalar condensate O χ modulations, which to leading order are ∝ cos(k x), to average out to zero.By the same argument the oscillations of the charge density ρ B also average to zero.Also, their period agrees with that of the scalar condensate, which is consistent with SDW order in a PDW.On the other hand, when q A q B = 0 the two uniform components χ (1) (r) and b (1) t (r) must be turned on, and we thus lose PDW order -the SC condensate now has a homogeneous component.
To understand the behavior of the charge density ρ A we must inspect the perturbations in the two channels (III) and (IV) of Appendix A. First of all, we emphasize that the modes a (1) t (r) and a (2) t (r) are sourced by the O(ε) (leading order) perturbations in δχ and δB t , and are therefore non-vanishing.Indeed, the ∼ cos(2k x) oscillation of ρ A is a next-to-leading order effect, which is sourced by the leading order oscillations ∝ cos(k x) of χ and B t .In this sense, it is induced, as needed in a phase with PDW order.Precisely because when q B = 0 the ∼ cos(2k x) oscillations in the scalar condensate are absent, the frequency of the ρ A oscillations is twice that of the scalar field.On the other hand, when q B = 0 the frequency of the oscillations of ρ A is the same as that of the scalar condensate 8 , which now has a uniform component.Thus, what we have is a co-existing SC+CDW state, and not a PDW.
Note that from the perturbation expressions (4.4) and the fact that ε ∼ 1 − T /T c , when q B = 0 we expect the following scaling behavior near T c , Here O 1 and ρ B1 are the cos(k x) components of O χ and ρ B , O 2 and ρ B2 are their cos(2k x) components, and finally ρ A1 is the cos(k x) component of ρ A .This is precisely observed in our numerics, as we show in Figure 3.In summary, the perturbative analysis suggests that the system behaves like a PDW when q B = 0 (no uniform component to O χ , and doubling of oscillation frequency), while when q B = 0 it describes a SC+CDW state (the uniform contribution ∝ χ (1) (r) is sourced and the frequencies match), as confirmed explicitly by the numerics.Table 1 summarizes the key differences between two cases.

VEVs
q B = 0 Comparison of various VEVs between q B = 0 and q B = 0.The former case corresponds to PDW order, while the latter to a state with co-existing SC+CDW orders.

The Striped Black Brane Solutions
In this section we construct the non-linear solutions corresponding to the striped black branes.To be more explicit, we take the couplings (2.8) to be given by which correspond to The critical temperature for the striped instability as a function of wave number for this choice of parameters is shown by the solid blue curve of Figure 1.For our numerical analysis it is convenient to introduce a new radial coordinate, in terms of which the AdS-RN black brane geometry (2.9) reads ( Note that the horizon is now located at z = 0 and the AdS boundary is at z = 1.The black hole temperature, expressed in terms of r h , is given in (2.10).

The ansatz and boundary conditions
We adopt the following ansatz for the striped black brane geometry, where the eight functions x. Notice that we recover the AdS-RN solution (5.3) by choosing their background values to be We will consider two different cases, corresponding to q B = 0 and q B = 0. Since we are looking for solutions with a regular horizon at z = 0, all functions depend smoothly on z 2 .Thus, at the horizon one can easily impose Neumann boundary conditions of the type ∂ z φ(0, x) = 0, and similarly for the remaining components.Additionally, we impose the Dirichlet boundary condition Q tt (0, x) = Q zz (0, x), so that the temperature of the black brane (5.4) is still given by (2.10).
We employ the DeTurck method [34] to solve the resulting PDEs.The key point behind the method is to introduce a new term into the Ricci tensor, where the DeTurck vector field is given by and Γµ λσ (ḡ) is the Levi-Civita connection associated with a reference metric ḡµν .We choose the reference metric to be the metric of the AdS-RN solution (5.3).To avoid Ricci solitons [34], we should ensure that the quantity ξ µ ξ µ is zero to the required numerical accuracy (see Appendix C).
To know the asymptotic behavior near the AdS boundary, we start from the scaling dimensions of the operators in the dual field theory.Labeling the eight unknown functions by i +δQ i , we are going to assume that the perturbations are of the form where v i denote constants.After substituting into the equations of motion, we obtain a matrix equation of the form where M ij is an 8 × 8 matrix that depends on the exponents δ i .In order to have non-trivial solutions we demand the matrix determinant to vanish, which in turn allows us to solve for the possible values of δ i , where δ 1,2,3,4,5,6 correspond to, respectively, the scalar, the two vector fields and the metric.The last two pairs δ ± 7,8 are due to the dynamical gauge fixing procedure [31].After imposing the source free condition for the scalar χ and vector field B t and fixing the boundary metric to be Minkoswski, the asymptotic expansion in the UV behaves as where we have imposed the following boundary conditions at z = 1, (5.10) The coefficients (φ v , ρ a , ρ b , q tt , q xx , q yy ) appearing in the boundary expansion are all dimensionless.Finally, the equations of motion lead to the following conditions which means that the energy momentum tensor of the dual field theory is traceless and conserved (see the discussion in Appendix B).
In our numerics we first use the standard pseudospectral collocation approximation to change the PDEs into non-linear algebraic equations.More precisely, we adopt Fourier discretization in the x direction and Chebyshev polynomials in the z direction.The resulted system is then solved using a Newton-Raphson method.

Numerical solutions
As expected from the analytical analysis, at sufficiently low temperatures the expectation value of the operator dual to χ is non-zero, and we have a phase with a striped condensate.As a typical example, throughout we are going to focus on the k/µ = 1 branch, with the corresponding critical temperature being T c = 0.016µ.We will work in the grand canonical ensemble with µ = 1 and choose the convention 2κ 2 N = 1.We plot the temperature dependence of the VEV of the scalar condensate in the left panel of Figure 4, for q B = 0.The right panel on the other hand shows the temperature dependence of the averaged free energy, which indicates that the broken phase is indeed thermodynamically preferred over the normal phase.Representative profiles for the bulk fields for q B = 0 are shown in Figures 5  to 8. Note that one can clearly see that the spatial modulations are imprinted on the horizon (at z = 0), and decrease in overall magnitude as the UV is approached.Indeed, the behavior of the strength of the striped oscillations, clearly visible from the figures, is in contrast to what is expected from explicit UV lattices of the type [35,36].In our model the striped feature is clearly strongest in the IR, and is therefore a relevant deformation of the UV field theory.
Figure 5 shows some of the functions appearing in the metric.The left panel of Figure 6 shows the dependence of the scalar field χ ∝ φ on both the holographic coordinate z and the "striped" spatial coordinate x, while in the right panel we plot the VEV of the scalar condensate as a function of x.We see from the latter that when q B = 0 the oscillations in the condensate O χ average out to zero, as required by a PDW state.In Figure 7 we plot the profile for the gauge field A t ∝ α (left panel) and for the corresponding charge density ρ A (right panel).The analogous plots for the second vector field are shown in Figure 8.We show the charge density profiles for the two vector fields as a function of the x coordinate in Figure 9, where they have been superimposed on top of each other.It is clear from the figure that when q B = 0 the  first vector field A oscillates twice as fast as the second one, and also twice as fast as the scalar condensate, confirming the other crucial feature of the PDW.
The results are qualitatively different when the second charge is non-zero, q B = 0. Bulk configurations of the vector fields and their corresponding charge densities for q B = 0.5 are shown in Figures 10 and 11.It is clear that both vector fields share the same period, which is in sharp contrast to the q B = 0 case.In this case one finds a homogeneous component to the scalar field condensate, which in turn implies that its oscillations no longer average out to zero.
The difference between the two cases, q B = 0 vs. q B = 0, is shown in Figure 12, where we plot the corresponding profiles for O χ superimposed on each other.The dotted line describes the profile for q B = 0.5, with the average value of the oscillations given by the horizontal dotted line.Finally, we plot the charge densities ρ A and ρ B when q B = 0 in Figure 13, for two different temperatures.The crucial difference with the case shown in Figure 9 is that the two densities now have the same period,  as already visible from the corresponding bulk fields in Figure 10 and Figure 11.Thus, we have reproduced the main features of a PDW (q B = 0) versus a CDW+SC (q B = 0) state.

Axions as Sources of Momentum Dissipation
The normal phase in our model, for which the gravitational background is the standard AdS-RN black hole, is dual to a translationally invariant system with a finite charge density.Therefore, it cannot relax momentum and would exhibit an infinite electric DC conductivity in the normal phase.Also, in our model translational invariance is broken spontaneously at low temperatures.Thus, even if the system were not superconducting, one would still expect to find generically an infinite conductivity, at least when the behavior of the Goldstone modes in the system is captured at large N by the gravitational description.Thus, to construct realistic models of dissipation in real materials, one needs to incorporate a mechanism for momentum relaxation.A simple way in holography is to introduce massless scalar fields which are linear in one of the spatial coordinates [37], where the strength of momentum relaxation is identified with the associated proportionality constant.Here we examine the role of such axionic scalars, on the details of the phase transition and on the background itself.We should stress that the mechanism for generating striped order and the charge density modulations is the same as described in the previous sections, and is not due to the axions9 .To this end, we would like to add two axions (ψ 1 , ψ 2 ) to our model (2.1), The normal phase is now given by with r h the horizon and µ the chemical potential.Notice that the axions depend on the spatial coordinates linearly, breaking translational invariance and giving rise to structures in two cases, with and without axions.Alternatively, instead of using axions one could introduce a lattice via a spatially modulated chemical potential.This was done recently in [38], where the authors observed the commensurate lock-in between the stripes and the underlying lattice.
momentum relaxation, whose magnitude in this setup is going to be controlled by the parameter τ .The only role of the axions is to affect the blackening function and therefore the black brane temperature, which is now given by The analysis of striped instabilities proceeds similarly to that of the standard AdS-RN background.To take into account the axions one must turn on the following modes up to second order, (1)  xx (r) + h (2)  xx (r) cos(2k x)] , δg yy = ε 2 [h (1)  yy (r) + h (2)  yy (r) cos(2k x)] .
Substituting them into the equations of motion and expanding to linear order in ε around the normal background (6.2), we obtain two coupled linear ordinary differential equations which are the same as those in (3.11), but with the original blackening function f replaced by f .As usual, the critical temperature at the onset of striped order is determined by looking for zero mode solutions.More precisely, we impose regularity conditions at the horizon r = r h and source free conditions at the AdS boundary r → ∞.The critical temperature will now depend on the strength of momentum dissipation τ , and its dependence on wave number is shown in Figure 14.We find bell curves similar to those in Figure 2.However, the presence of momentum dissipation suppresses the striped instabilities as compared to the case in which τ = 0 (the top curve in the figure).Indeed, notice that as τ increases the magnitude of the peak decreases, and at the same time its location shifts to the right, i.e. to larger values of k.We note that a similar feature in the behavior of T c was seen in the axionic model examined in [39].While it would be interesting if this was a somewhat universal feature, it is not clear to us at this stage whether this is the case, and we suspect that it is instead model dependent.
We now construct the fully backreacted geometry that describes the striped phase.We choose the same couplings as (5.1).Again, it is convenient to use the new coordinate (5.2), in terms of which the background in the normal phase is still given by (5.3) but with F replaced by the quantity explicitly impacted by the axions.The horizon and boundary are still located, respectively, at z = 0 and z = 1, and the temperature is given by (6.3).
Our ansatz for the striped black brane geometry is where the nine functions The normal solution (6.5) is found by setting , and once again we look for solutions with a regular horizon at z = 0, so that all the unknowns depend smoothly on z 2 .As we did for the case without axions, we impose one additional Dirichlet condition Q tt (0, x) = Q zz (0, x) such that the temperature of the black brane is still given by (6.3).In our numerics, we impose the following boundary condition at the conformal boundary z = 1, In particular, notice that the last condition ψ(1, x) = 0 means that the leading terms which are dual to the sources are just ψ 1 = τ x, ψ 2 = τ y, the same as in the normal phase.In the spatial direction x, we impose periodic boundary conditions.The coupled PDEs are solved using the DeTurck method, and representative solutions for the bulk fields are given in Figure 15.Their behavior qualitatively matches the perturbation analysis in (6.4).While the phase structure is not impacted in any qualitative way by introducing ψ 1 and ψ 2 , we expect to find interesting changes in the transport behavior in the system.
Indeed, we note that -in the absence of a superconducting condensate-when the translational symmetry breaking is spontaneous, the electrical conductivity should be divergent even along the direction of symmetry breaking [40].This can be explained by the presence of Goldstone modes in the system, which can be excited at no cost.Without any explicit breaking of translational symmetry, the momentum relaxation rate Γ and pinning frequency ω o are both vanishing, so the DC conductivity is infinite.Our model with axions provides a simple way to introduce a finite momentum relaxation rate Γ.If the U(1) symmetry is preserved, so that one doesn't have a superconducting condensate, then the axions would get rid of the divergence in the conductivity due to the Goldstone modes.However, when the U(1) symmetry is broken and the condensate forms, we expect to obtain an infinite DC conductivity, associated solely with the superconducting excitation.It will be interesting to see if one can see the pinning effect due to the CDW excitation from the behavior of the optical conductivity, in the presence of such axionic contributions.

Final Remarks
Our goal in this paper was to examine a gravitational model which would realize holographically certain features of a PDW, a phase in which charge, superconducting and spin orders are intertwined with each other in a spatially modulated manner.Specifically, we wanted to obtain a quantum system in which the scalar condensate and charge density are both striped, and are such that the average value of the superconducting order parameter O χ vanishes (i.e.there is no uniform component), and the frequency of the charge density modulations is twice that of the condensate.Our model indeed realizes both of these features -two defining properties of PDW order -by spontaneously breaking the global U(1) symmetry of the boundary theory, at the same time as translational invariance.Thus, the resulting ordered phases are intertwined with each other and are understood to have a common origin.
We recall that our construction (2.1) is not of the form of the standard holographic superconductor, but instead falls within the generalized class of models advocated for in [9], in which the spontaneous breaking of the global U(1) symmetry happens via the Stückelberg mechanism.The local gauge invariance is now encoded in the transformation (2.2), and in these theories the two scalars χ (the field which condenses) and θ (which can be gauged fixed to zero) are not necessarily associated with, respectively, the magnitude and phase of a complex scalar field.In particular, in order to realize PDW order in our setup, χ must oscillate about zero, and is therefore allowed to be negative -it cannot be directly identified with the magnitude of a complex scalar Ψ = χe iθ .In turn, this means that θ is not a priori restricted to be periodic and therefore the global U(1) symmetry that is being broken is non-compact (without imposing additional ad hoc restrictions on the range of θ).
Even though locally we expect to see the same features as in a theory with a compact U(1) symmetry (infinitesimally the two are the same), the global properties will be different, and clearly one should be working with theories in which the symmetry is compact, as is usually the case.We expect to be able to modify our model to ensure that the U(1) is indeed compact, and still obtain a PDW phase, but we anticipate that this requires adding a more involved matter sector.Thus, we view our construction as only a first step towards realizing more sophisticated and realistic phase diagrams involving PDW order.We have adopted this particular toy model because of its tractability, and take it as a proof of principle that the PDW features outlined above can indeed be realized.Guided by the intuition developed here, the model can then be refined and improved to be more relevant for phenomenology, and to yield more realistic order parameters, along the lines of e.g.[41].We postpone such efforts to future work.
Recall that the lack of a uniform component of O χ in our construction was associated with taking the charge q B to be vanishing.When q B = 0 on the other hand O χ contains a homogeneous, non-oscillating piece, and the CDW oscillations have the same period as those of the scalar condensate.Thus, we have co-existing SC and CDW orders, as opposed to a PDW.While the density ρ A associated with the first vector field A µ has the natural interpretation of a charge density (leading to CDW order), the physical interpretation of the density ρ B of B µ is less clear.Naively, it obeys some of the properties of a SDW, i.e. it is also striped and always oscillates at the same frequency as the scalar condensate.However, to be able to claim that ρ B is really describing a density of spins (leading to SDW order), one must justify physically the particular couplings and interactions chosen in our model, which are needed to realize these symmetry breaking patterns.Perhaps more natural couplings can be introduced, which would allow for a more realistic identification between the degrees of freedom of our model and those in a real high T c superconductor exhibiting these kinds of intertwined orders.However, the second gauge field can also be thought of as a "spectator" field, whose only role is to seed the oscillations of the charge density ρ A .Indeed, it is important to emphasize that in our model CDW order is induced from the oscillations of the scalar condensate and of ρ B , as can be understood at least qualitatively by inspecting the perturbation analysis.
An interesting question is that of the behavior of the conductivity in these systems.As we discussed in Section 6, the high temperature phase in our model is dual to a translationally invariant system, which lacks a mechanism to dissipate momentum.Moreover, the low-temperature phase describes a superconducting condensate and therefore exhibits an infinite conductivity.However, even if we turned off the scalar condensate and focused only on the CDW state of the broken symmetry phase, we would still expect the conductivity to be generically infinite.This can be traced back to the fact that the translational symmetry in our model is broken spontaneously -the resulting Goldstone modes can in principle be excited at no cost, leading to an infinite conductivity.To remedy this problem we have added axions to our model -a simple way to generate momentum dissipation in the system -and have examined the full backreaction.While the axions don't affect in any qualitative way the PDW (or SC + CDW) features explored in this paper, we expect them to play a crucial role in the transport properties of the system.In particular, if we were to restrict our attention to the case in which the U(1) symmetry is preserved, i.e. no superconducting condensate, we would expect the axionic case to lead to a finite conductivity.It would be interesting to explore the conductivity in this framework (as well as the pinning effect), and how it is impacted by the presence or absence of an explicit breaking of translational invariance.
Finally, we close with a few words on the nature of the ground state.Models of the type we are considering are known to support hyperscaling violating geometries, which can be obtained by suitable choices of scalar couplings and potential.Such backgrounds are well-known to develop both striped [42,43,44,45,46] and superfluid instabilities (see [47] for general criteria).It would be interesting to generalize our study of intertwined orders to this case.We leave these questions to future work.
h (1)  xx + h (1)  yy + 1 2 2 , (A.7) 5 .(A.10) Here the five source terms C i (i = 1, ..., 5) depend on the leading order modes (w, b t ).We haven't explicitly written down their form for simplicity, as they are not essential to our arguments.Note that (A.10) is a constraint equation, constraining e.g.h (1) tt .Moreover, the tt component of Einstein's equation (A.7) is implied by the remaining equations.Thus, what we have can be reduced to four independent differential equations: three second order equations for (a xx , h yy ) and one first order equation for h (1) tt .

A.4 Set (IV)
This case comprises six equations of motion for the four modes (a 2 , (A.12) (A.16) As in the previous case, the six source terms C i (i = 1, ..., 6) depend on the two leading order perturbations (w, b t ).The precise form of the source terms is not important for our arguments and is not included for simplicity.Again, the tt component of Einstein's equation (A.12) can be obtained from the remaining equations.Using the rr and rx components of (A.15) and (A.16), we find that h (2) tt satisfies an algebraic equation.Finally, we obtain three independent second order equations for (a

B. Thermodynamics in Grand Canonical Ensemble
Recall that in order to obtain a renormalized action, we must supplement it with appropriate boundary terms.For the model we are considering in Section 5, such terms take the form where γ µν is the induced metric at the AdS boundary and K µν is the extrinsic curvature defined by the outward pointing normal vector to the boundary n µ .According to the holographic dictionary, the stress-energy tensor of the dual field theory is Notice that we are interested in the one point functions with respect to the metric Plugging in (5.9), we obtain Working in the grand canonical ensemble, the local free energy density Ω(x) is given by Ω(x) = T tt (x) − T s(x) − µ ρ A (x), (B.9 where the thermal entropy density s is given by where χ s is identified as the source of the dual scalar operator O χ .Therefore, the one point function for O χ is computed by (B.14) In the grand canonical ensemble the two free thermodynamic variables are the chemical potential µ and temperature T .The system (5.4)As one can check, they only depend on the dimensionless combination r h /µ, or equivalently T /µ.

C. Accuracy of Numerics
In this appendix we check the quality of our numerics.We start by checking the convergence of ξ 2 = ξ µ ξ µ , where the DeTurck vector field ξ has been defined in (5.6).
The number of grid points in the z and x directions are N z and N x , respectively.As one can see in Figure 16, the maximum value of ξ 2 converges towards zero so that indeed we have a solution to Einstein's equations.We also find an approximatively exponential rate of convergence with the number of grid points, which is expected by the pseudospectral collocation method.Next, we perform a check of the first law of thermodynamics.Note that we are working in the grand canonical ensemble with the free energy given by (B.9).Together with the first law,

Figure 1 :
Figure1: Critical temperature versus wave number for the onset of striped instabilities for the massless case (solid line) and massive case (dashed line).The peaks are located at (k ≈ 1.130µ, T ≈ 0.01649µ) for the solid curve and at (k ≈ 1.048µ, T ≈ 0.01281µ) for the dashed one.The remaining parameters are m 2 = −8, L = 1/2, a = 4, κ = q A = 1.

Figure 2 :
Figure 2: Critical temperature as a function of wave number, for different choices of coupling c.From top to bottom the coupling is increased.We have m 2 = −8, m 2 v = 0, L = 1/2, a = 4, κ = q A = 1.

Figure 3 :
Figure 3: The log-log plot of the VEVs near T c .Straight lines denote the expected scaling behaviors shown in (4.5), and dotted lines denote numerical values.We consider the model (5.1) with q B = 0.5 and k/µ = 1.The data at the top corresponds to O 1 and ρ B1 , the data in the middle to O 2 and ρ B2 , and at the bottom to ρ A1 .

Figure 4 :
Figure4: The left plot shows the VEV of the scalar condensate at x = 0 as a function of temperature.The right plot shows the difference in the average value of the free energy Ω between the striped phase and the normal phase.The parameters chosen here are k = µ = 1, q B = 0 and the critical temperature is T c = 0.0161.

Figure 5 :
Figure 5: Bulk solutions for the metric fields Q tt , Q xx , Q xz , Q zz for T = 0.01426 and k = 1, for the parameters chosen in (5.1) and q B = 0.The horizon is at z = 0 and the AdS boundary at z = 1.

2 Figure 6 :
Figure 6: Profile for the scalar field φ (left panel) and scalar VEV (right panel) for T = 0.01426 and k = 1, for the parameters chosen in (5.1) and q B = 0.

Figure 7 :
Figure 7: Profile for the gauge field α (left panel) and charge density ρ A (right panel) for T = 0.01426 and k = 1, for the parameters chosen in (5.1) and q B = 0.

2 Figure 8 :
Figure 8: Profile for the gauge field β (left panel) and charge density ρ B (right panel) for T = 0.01426 and k = 1, for the parameters chosen in (5.1) and q B = 0. Note that the average of ρ B is zero.

Figure 9 :
Figure9: The two charge densities for q B = 0, T = 0.01426 and k = 1.Here the frequency of ρ A is twice that of ρ B (and consequently of the scalar condensate).

2 Figure 10 : 2 Figure 11 :
Figure 10: Profile for the gauge field α (left panel) and charge density ρ A (right panel) for T = 0.01426 and k = 1, for the parameters chosen in (5.1) and q B = 0.5.

2 Figure 12 : 2 Figure 13 :
Figure12: The scalar condensate for T = 0.01426 and k = 1.The solid blue curve corresponds to q B = 0 and the dashed one is for q B = 0.5.Two horizontal lines denote their average values.

Figure 15 :
Figure 15: Bulk solutions for T = 0.00737µ and k = 1 in the presence of axions.The horizon is at z = 0 and the conformal boundary at z = 1.We have chosen the momentum relaxation strength τ /µ = 0.5 and the other parameters are the same as (5.1).
tt .Since A t corresponds to the charge density ρ A of the dual field theory, the mode a (2) t means ρ A is spatially modulated with wave number 2k.

ds 2 = 2 N T tt = r 3 hL 4 2 h+ 3 (
−dt 2 + dx 2 + dy 2 .(B.3) Substituting the expansion (5.9), we obtain the following components of T µν , 2κ q tt + q xx ) ,(B.4)with all other components vanishing.Moreover, using the relation(5.11),one can explicitly check that the stress-energy tensor is traceless and conserved,T µ µ = 0, ∂ ν T νµ = 0, (B.5)which is consistent with the expectation that it describes a CFT.The expectation values of the currents can be obtained by the variation of the on-shell action with respect to the gauge fields, ν (Z A A µν + Z AB B µν ), ν (Z B B µν + Z AB A µν ).(B.7)

3 h 2κ 2 N L 4 2 + L 2 2 µ 2 r 2 h 2 N r 2 hL 2 , ρ A = 1 2κ 2 N 12 )= χ s r + χ v r 2 +
x)Q yy (0, x).(B.10)For the normal solution (5.3), we then haveT tt = r , s(x) = 2π κ r h µ L 2 , (B.11)and the corresponding free energy density is given by Ω Notice that the UV expansion for χ in the r coordinate is in general given by χ O(r −3 ), (B.13) d T tt (x) = T ds(x) + µ dρ(x) , on a fixed k/µ branch of the striped solutions.As shown in Figure17, this condition is satisfied along the set of solutions we are considering in Figure4.

2 Figure 16 :Figure 17 :
Figure 16: ξ 2 as a function of N z = N x = N for T = 0.01426 and k = 1 with µ = 1.We consider the same parameters as in Figure 4.