Gravity dual of spin and charge density waves

At high enough charge density, the homogeneous state of the D3-D7' model is unstable to fluctuations at nonzero momentum. We investigate the end point of this instability, finding a spatially modulated ground state, which is a charge and spin density wave. We analyze the phase structure of the model as a function of chemical potential and magnetic field and find the phase transition from the homogeneous state to be first order, with a second-order critical point at zero magnetic field.


Introduction
The ground states of many condensed matter systems feature spontaneously broken spatial symmetries. Modulated states, such as charge and spin density waves, break translation symmetry and some or all of the symmetries of the underlying lattice. Systems that exhibit these phases, such as doped Mott insulators, are typically characterized by competing types of order and often have complicated phase diagrams. The pseudogap regime of cuprate superconductors, for example, features a variety of symmetry-breaking phases.
Gauge/gravity duality has proven to be a powerful tool for investigating stronglycoupled systems, such as those arising in many interesting condensed matter systems. Although most of the focus has been on computationally simpler homogeneous and isotropic systems, in recent years symmetry-breaking instabilities and states with broken symmetry have received increasing attention. One motivation to consider spatially inhomogeneous systems is that breaking spatial symmetries can affect inter-esting physical quantities. For example, translational symmetry implies momentum conservation which automatically leads to a divergent DC conductivity. To obtain physically sensible results, translational symmetry must be broken.
In many cases, one is driven to investigate inhomogeneous states because the naive homogeneous state turns out not to be the ground state. Investigations have found a range of examples in which homogeneous states suffer from instabilities [20][21][22], and a lot of work has gone into the difficult problem of finding the inhomogeneous end states of such instabilities [23][24][25][26][27][28][29][30][31].
In this paper, we will find the ground state of the D3-D7' probe-brane model. The homogeneous phase of this system has been analyzed in great detail, both in the metallic, black hole-embedding phase [36][37][38] and in the quantum Hall, Minkowskiembedding phase [39][40][41]. 1 A fluctuation analysis [37] showed that, above a critical charge density, the homogeneous black hole phase is unstable. The tachyonic modes responsible for this instability have nonzero momentum, suggesting the system decays to an inhomogeneous ground state. This instability, however, is mitigated by a nonzero mass for the fermions or an external magnetic field. This perturbative result motivated us to search for spatially dependent solutions and to investigate whether the expectations from the fluctuation analysis hold at the full nonlinear level.
We find explicit numerical solutions for configurations which exhibit spontaneous modulation in one spatial dimension. The bulk fields are entirely coupled; as a result, the corresponding inhomogeneous state is both a spin and charge density wave. The spatial frequency of this striped phase is closely related to the momenta of the tachyons of the homogeneous phase. Comparing with the homogeneous phase, this striped phase is thermodynamically preferred at sufficiently large charge density and small magnetic field. In general, the phase transition is first order, although there is a second-order critical point at zero magnetic field.
In the next section, we will review the construction of the D3-D7' model. Section 3 will describe the pseudospectral method we employ to find our modulated solution. In Section 4, we will describe the solutions, the transition from the homogeneous phase, and the phase diagram. We finish in Section 5 with a discussion of open questions and directions for future research.

Review of the D3-D7' model
We begin by setting up the notation and recalling the model under study. We are interested in a holographic quantum liquid modeled by a probe D7-brane embedded in a black D3-brane background. The intersection of the D3-D7' system is (2+1)dimensional, supersymmetry is completely broken, and the low-energy excitations of the field theory dual are purely fermionic. The model thus makes a compelling case for condensed matter applications.
The ten-dimensional background metric reads: where the thermal factor is h = 1 − r T r 4 . The metric on the internal sphere we write as an S 2 × S 2 fibered over an interval: with the ranges for angles ψ ∈ [0, π/2], θ, α ∈ [0, π], and φ, β ∈ [0, 2π]. Recall also that the Ramond-Ramond four-form is C (4) txyz = −r 4 /L 4 AdS , and the curvature radius is related to the 't Hooft coupling as In the following we will work in units where L AdS = 1.
We will embed the D7-brane such that it spans t, x, y Minkowski directions, is extended in the holographic radial direction r, and wraps both of the two-spheres. We are interested in solutions which depend on x and r, such that the embedding functions are z = z(x, r) and ψ = ψ(x, r). Moreover, to creep towards more realistic physical models, we wish to turn on all the external components of the gauge potential along the brane worldvolume: a t,x,y,r (x, r) = 0. Since the D7-brane would otherwise be unstable towards slipping off the internal space, we also turn on internal components of the gauge field, such that their field strengths are: The quantized constants f 1 and f 2 are proportional to the number of D5-brane fluxes diluted on the internal two-spheres. The dynamical fields we have are ψ, z, a 0 , a x , a y , a r . Since we are only interested in time-independent solutions, we can work in the gauge a r = 0.
The action, which consists of a Dirac-Born-Infeld term and a Chern-Simons term, reads: where the prime denotes differentiation with respect to r, and The overall constant factor of the action is N = 4π 2 T 8 V 1,1 , where V 1,1 is the spacetime volume in the homogeneous t and y directions. Notice that a x only appears once, as a single term ha 2 x in A, which means that there is a trivial constant (both in x and r) solution for it and it decouples from the rest of the dynamics; we can thus drop a x in what follows. Before proceeding to solve the equations of motion that follow from the above action, let us make a coordinate transformation and in particular absorb one scale.
The model has several parameters. These correspond to some boundary values (r → ∞) of the bulk fields ψ, z, a 0 , a y , as well as the fluxes f 1 and f 2 , together with the parameter r T of the background which is proportional to the temperature. Not all the parameters are independent, however. It was shown in [36] that f 1 and f 2 fix ψ(r → ∞) ≡ ψ ∞ via the relation The fluxes also determine the subleading exponent of the field ψ ∼ ψ ∞ + Ar ∆ ± : where we have used (11) to eliminate f 2 . We note that, in order for the model to be perturbatively stable, one needs large enough flux to have real exponents ∆ ± . In this paper we will only consider the case f 1 = f 2 = 1/ √ 2 which implies that ψ ∞ = π/4 and also that the anomalous mass dimension of the fermions is vanishing: where m ψ is proportional to the quark mass and c ψ is proportional to the condensate. For simplicity, we will consider massless fermions in this paper, m ψ = 0.
We can freely set the boundary values of z and a y to zero since neither z nor a y enter in the equations of motion without derivatives and they do not also have IR constraints. This is unlike for the case of a 0 , whose IR value has to vanish at the horizon, leaving the UV boundary value a 0 (r → ∞) ≡ µ as a physical parameter of the theory. We will specify the UV boundary conditions explicitly for the different fields in App. A.1.
We have a freedom to scale out one parameter of the model and for this we choose to pick the temperature. We will introduce a new compact radial variable To get rid off explicit dependence of r T factors everywhere, we introduce notation with hats as follows:x Essentially, all the parameters are read in the units of the horizon radius.
The action in the new radial coordinate reads where 2 G = f 2 1 + 4 cos 4 ψ f 2 2 + 4 sin 4 ψ (18) and The new overall factor is N T = r 2 T N . The prime now denotes differentiation with respect to u.
The equations of motion that we will study in this article are the Euler-Lagrange equations for the fields ψ = ψ(x, u),ẑ =ẑ(x, u),â 0 =â 0 (x, u), andâ y =â y (x, u), which follow from (17). To streamline the discussion and because they are quite lengthy, we will not write the equations of motion down explicitly.

Solution by pseudospectral method
We will now briefly sketch the techniques we used to solve the equations of motion, which are a coupled, nonlinear, second-order PDE system. A more detailed description can be found in Appendix A.
We employed a pseudospectral method (see, e.g., [45]), in which the difficult partial differential equations are rendered into a more tractable system of algebraic equations. Most of the runs could be performed on a personal computer, but we ended up using a supercomputing cluster for longer runs in order to achieve sufficient accuracy to analyze the phase structure.
The first step is to determine the correct set of boundary conditions. The boundary conditions imposed must be consistent with the equations of motion and must fix all constants of integration. We are interested in the spontaneous breaking of translation invariance rather than breaking it explicitly. Therefore, the nonnormalizable terms in the UV remain as unmodulated free parameters. The normalizable terms will be, in general, modulated, and their values are fixed by requiring regularity in the IR.
In the UV, at u = 0, we fix the unmodulated sources as follows (see Appendix A.1 for more detail): whereμ is the (reduced) chemical potential for the fermions of massm ψ , andb is the (reduced) external magnetic field: 3m At the horizon u → 1, we demand that the solutions are regular and finite. One particular requirement for the regularity is thatâ 0 is constant in the IR but in fact we need to require this constant to vanish: which ensures thatâ is also a well-defined one-form as the thermal circle pinches off.
Finally, we need to specify the boundary conditions in the x-direction. We are looking for spatially modulated solutions, so it is natural to impose periodic boundary conditions; that is, ψ(x) = ψ(x +L), etc, where L = r TL is the spatial period. However, it is not obvious a priori what is the correct periodL (or equivalently frequencyk = 2π L ) to choose. Our approach is to find solutions for a range of spatial periods, and by comparing their free energies, determine which wavelength is preferred by the system.
Having fixed the boundary conditions, we proceed to turn the system of differential equations of motion into a system of algebraic equations. To do this, we implement the standard pseudospectral method based on a Fourier series having N x terms in thex-direction and an expansion in the Chebyshev polynomials having N u terms in the u-direction. The values of ψ,ẑ,â 0 , andâ y at the O(N x N u ) collocation points are chosen as the variables, and the derivatives are computed in the pseudospectral approximation, i.e., from the Fourier and Chebyshev series with O(N x N u ) terms which exactly match with the variables at the collocation points. We then insert these numbers into the equations of motion. By demanding that the equations are satisfied at the collocation points, we obtain O(N x N u ) algebraic equations for the values of the functions which we then solve numerically. See Appendix A for more details. The accuracy converges exponentially with increasing N x and N u . We tested the code up to N x = 38 and N u = 40, which yielded an accuracy of at least 10 −9 for all functions.

Striped solutions without a magnetic field
Having set up the numerical method for solving a system of coupled PDEs, we can vary the parametersμ,L, andb and look for striped solutions. We consider massless fermions, som ψ = 0. For simplicity, we begin with the restricted case of zero magnetic field and generalize tob = 0 in Sec. 4.2.
A representative striped solution is shown in Fig. 1. All the fields are modulated: ψ andâ y have periodicity equal to the periodicity imposed on the spaceL soln =L but are π/2 out of phase with each other;â 0 andẑ are modulated with frequency 2k. Because of the charge density and wrapped internal flux, the bulk fieldsâ 0 andẑ have nontrivial homogeneous profiles; in Fig. 2, we have subtracted out this homogeneous background to highlight the modulation. By construction, the modulation goes to zero near the boundary u → 0, indicating that the translation symmetry is being spontaneously broken. To determine the preferred modulation frequency, we compare the free energies of solutions of varyingL, withμ held fixed. As we are working at fixedμ, that is, in the grand canonical ensemble, the relevant free energy is the grand canonical potential, which is defined by the Euclidean action (17) evaluated on a solution to the equations of motion: This divergent expression (31) is regulated using the standard holographic renormalization; the required counterterms are explicitly written in [36]. Note that the range of thex integral in the action isL. In order to compare solutions with different periodicitiesL, the quantities we need to compare are the average grand potential densities Ω/L. Bear in mind that, while we are artificially fixing the spatial period-icityL by hand here, in the actual system the periodicity is determined dynamically as the system minimizes its energy.
In Fig. 3, we plot Ω/L as a function ofk = 2π L for a fixedμ. The various curves on the plot correspond to different ratios between the periodicity ofx, which isL, and the periodicity of the solution, denoted byL sol . For the solution shown in Fig. 1, the two periods are equal, but in general there may be any integer number of modulations of the solution within the spatial period, i.e.,L = nL sol for any positive integer n. In particular, if a field configuration with modulation periodL sol is a solution when L =L sol , it will necessarily also be a solution forL = nL sol .
As can be seen from Fig. 3, there is a minimum energy solution at nonzerok =k 0 ; this is the frequency preferred by the system for a givenμ and is the frequency the system will take whenL is not fixed by hand. Forcing the system to have a higher frequencyk >k 0 costs energy, and the amplitude of the modulation decreases, as shown in Fig. 4. For sufficiently largek =k max , the modulation vanishes, and the In the opposite limit, if we make the spatial periodicity largeL 2π k 0 , the system will try to have as many oscillations as necessary sok is as close tok 0 as possible. Ask is decreased belowk 0 , the energy similarly grows and the amplitude decreases, eventually reaching a minimumk =k min .
Rather than merging with the homogeneous solution, the n = 1 branch of solutions merges atk min with a branch of solutions which interpolates between the n = 1 and n = 3 solutions (shown in red dots in Fig. 3). Interestingly, along this interpolating branch of solutions is a local minimum of Ω/L, indicating the presence of a metastable solution.
We now allowμ to vary, and compute Ω/L as a function of bothμ andk. In Fig. 5, we plot the difference between the modulated and homogeneous free energies ∆Ω/L = (Ω mod − Ω hom )/L. Forμ above the critical chemical potentialμ c , there is a range ofk where ∆Ω/L < 0, indicating that the modulated solution is preferred. The preferred frequencyk 0 (μ) is given by thek that minimizes ∆Ω(μ)/L and sets the frequency of the striped phase. As shown in Fig. 5, the striped phase appears at µ c with a nonzero critical frequencyk c =k 0 (μ c ), andk 0 (μ) increases withμ. The transition between the homogeneous and the striped phases occurs atμ c . Aŝ µ decreases toμ c , ∆Ω(k 0 ,μ)/L → 0 smoothly, as illustrated in Fig. 6(left). Atμ c , the charge densityd = − ∂Ω ∂μ of the two phases are equal, as shown in Fig. 6(right), indicating that the phase transition is continuous. We can fit to very good accuracy the difference in densities of the two phases ∆d =d mod −d hom nearμ c to a linear function ofμ. We conclude from this that the phase transition is second-order. We can further investigate the nature of the phase transition by analyzing how the modulation disappears as the critical point is approached. Fig. 7 shows the amplitude of modulation, which is well fit near the critical point by √μ −μ c for the spin density wave (ĉ ψ ) and by a linear function for the charge density wave (d). This phase transition exactly matches the expectations obtained from the fluctuation analysis of the homogeneous phase. As shown in [37], at a particular chemical potentialμ inst , one quasinormal mode develops a positive imaginary dispersion at nonzero momentumk inst ; see Sec. 4.4 of [37], 4 particularly Fig. 4. The chemical potential and momentum at which the instability occurs match the critical chemical potential and spatial frequency computed here from the explicit modulated solutions: µ inst =μ c andk inst =k c .
There is another lesson to be drawn. Comparing the magnitudes of the modulations of different fields, for example as depicted in Fig. 7, we see that the amplitude of the spin wave dominates over that of the charge density wave. This conforms nicely with the fluctuation analysis [37]. When the fermion mass is set to zerom ψ = 0 at zero magnetic field strengthb = 0, the fluctuation equations decouple. The tachyonic mode responsible for the modulated instability resides entirely in the sector which only mixes the embedding scalar δψ fluctuation and the transverse gauge field fluctuation δâ y , thus suggesting only a spin density wave. Comparison of the plots in Fig. 7 shows, however, that coupling to the charge density wave is present at second order in the fluctuation analysis. It is also interesting to note that the nonlinearities do not essentially ramp up the charge density wave, even further away from the critical point.

Stripes in a magnetic field
We now turn the magnetic field back on and again solve the PDEs. As in theb = 0 case, we find solutions with all the fields modulated. The main qualitative difference is that thex →L 2 −x parity symmetry is now broken. For eachb and withμ fixed, we again minimize the free energy Ω/L to find the spatial frequencyk 0 of the striped phase. We find that decreasingb causesk 0 to increase. The energy difference ∆Ω/L with the homogeneous phase at fixedμ is plotted as a function ofb andk in Fig. 8.
We find that the magnetic field suppresses modulation. For a givenμ, modulation is thermodynamically preferred below a critical magnetic fieldb c , appearing at a critical frequencyk c > 0. Fig. 9 shows how the free energy of the striped phase compares with that of the homogeneous phase. Nearb c , the energy difference is a linear function ofb. The magnetizationM = − ∂Ω ∂b therefore jumps atb c , indicating that the the phase transition is first order, in contrast to theb = 0 case; compare Fig. 9 with Fig. 6(left).
By varying bothμ andb, we can map out a two-dimensional phase diagram, which is shown in Fig. 10. For eachμ, we compute the critical chemical potentialb c , which at nonzerob c is the location of a first-order phase transition. This line of first-order transitions ends atb c = 0 at a critical point, where the transition becomes continuous.
We can compare this phase diagram Fig. 10 with the results of the fluctuation analysis of the homogeneous phase at nonzerob, performed in [38]. Because the phase transition is in general first order, we would expect the homogeneous phase to change only from stable to metastable at the transition pointμ c ; the instability should appear atμ inst >μ c . In fact, this is what is found by comparing our results The blue curveb c (μ) is a curve of first-order phase transitions between the homogeneous phase (above the curve) and the striped phase (below the curve). This curve ends in a second-order critical point atb c = 0 andμ = 2.66, indicated in green. The homogeneous phase is perturbatively unstable below the red curve; between the two curves, it is metastable.
to those of in Sec. 4.2 and Fig. 5 of [38]. 5

Discussion and outlook
In this paper, we presented the striped phase of the D3-D7' system indicated by the perturbative instability of the homogeneous state in [37,38] and found by numerically solving the coupled equations of motion. This modulated phase was found to be both a charge density and spin density wave, while the latter component was vastly more dominant. The phase structure presented in Fig. 10 is in excellent agreement with the previous perturbative analysis.
There have been several recent examples of spontaneously modulated ground states of a similar type as those studied here; in particular, [22,25,[27][28][29] analyzed bottom-up holographic models whose homogeneous states suffer from similar types of instability as the D3-D7' system. They find striped, cohomogeneity-one solutions, dual to charge and/or current density waves. Though these solutions seem to resemble those that we found in this work, our key player was the spin density wave and not so much the charge density wave component of the modulation. Also, the phase transition from the homogeneous to striped phases in those works were typically found to be of second order, though those systems were not studied in the presence of a magnetic field.
There are a number of interesting avenues left to explore in the future. The D3-D7' system features a number of parameters which have been fixed so far, such as the mass m ψ and the internal fluxes f 1 and f 2 , and which could be varied to potentially interesting effect. In particular, the alternative quantization of the bulk gauge field, studied in [40,41], corresponds to an SL(2, Z) mapping of the boundary CFT. An upcoming stability analysis [46] indicates that the choice of quantization strongly affects the phase structure, but a construction of modulated solutions will be needed for a complete understanding.
In the quantum Hall phase, corresponding to a Minkowski embedding of the D7brane, the lowest neutral excitation, for a range of parameters, is a magneto-roton, that is, a collective mode whose minimum energy is at nonzero momentum [39]. One could look for modulated Minkowski solutions, corresponding to a roton condensate. However, since the homogeneous solution was found to be perturbatively stable, such a modulated state would not likely be preferred energetically.
A compelling but technically difficult open question is whether the striped phase is in fact the ground state. Or, does continuous translation symmetry break completely, resulting in a two-dimensional lattice? In principle, one could solve the system of PDEs, allowing the fields to depend on both x and y. This presents a challenging numerical problem. 6 Another, possibly more tractable approach, might be to perform a fluctuation analysis of the striped phase to detect perturbative instabilities. Just as a tachyon with nonzero momentum signaled an instability of the homogeneous mode to forming stripes, a tachyonic mode with nonzero momentum along the stripes would imply an instability toward forming a lattice. Such a calculation would involve solving the linearized fluctuation equations on top of the numerical striped background.
Another computation remaining is the calculation of the conductivity of the striped phase. Translationally invariant phases typically feature infinite DC conductivity, though probe-brane models have finite DC conductivity because momentum can be dissipated into the large N number of bulk degrees of freedom. Still, the conductivity in the modulated direction is expected to have interesting properties.
There are two standard methods for computing conductivities in holographic probe-brane models. One can compute the retarded current-current correlator which relates to the conductivity via a Kubo formula. 7 In this case, such a calculation would be similar in difficulty to the linear stability analysis discussed above. Alternatively, the Karch-O'Bannon method [48] can be used to compute the DC conductivity. However, because the embedding is inhomogeneous and only known numerically, em-ploying this technique will not be straightforward.
Finally, there are many basic open questions to be addressed. Of particular importance is to what degree symmetry-breaking is generic. Symmetry-breaking instabilities appear to be a common feature of holographic systems. To what degree can the result of this model be generalized? Given the current nature of holographic modeling, the detailed features of any one construction is less relevant than generic features shared by a wide class of systems.
For example, the magnetic field tends to inhibit the symmetry-breaking instability in the D3-D7' system. A similar effect was also seen in the linear stability analysis of the closely related D2-D8' model [43], as well as in the Sakai-Sugimoto model [49]. We might be tempted to conjecture that magnetic fields tend to stabilize modulated instabilities more generally. This type of statement is about as much as we can currently aim for in holographic condensed matter physics. The instabilities in these models are all of the type identified in [50], so this may fail generalize to instabilities generated by different mechanisms. So, it would be interesting to investigate the stabilizing effects of a magnetic field in other holographic models, to test this conjecture.

A Finding inhomogeneous solutions numerically
In this appendix, we describe in more detail the pseudospectral method for numerically solving the equations of motion.

A.1 Boundary conditions
Let us start from the boundary conditions in the UV, where u → 0. After using the radial gauge condition and rotation symmetry in the x-y-plane, the dynamical fields are ψ,ẑ,â 0 , andâ y . One can solve the equations of motion as a series around u = 0, assuming unmodulated sources: a y (x, u) =ā y +bx +ĵ y (x)u + 1 48 The constantsm ψ ,μ, andb are parameters we are free to choose. The translation symmetry of the background transverse to the D7-brane means that the constants z andā y can be set to any value without affecting observables, so we can set them to zero without loss of generality. Therefore, we impose the following UV boundary conditions: These free parameters have clear physical interpretations. Theμ is the (reduced) chemical potential for the fermions of massm ψ , andb is the (reduced) external mag-netic field:m In the IR, at the horizon where u → 1, we require that the solutions are regular and finite. In particular, we notice that the action (17) contains, inÂ x of (20), a term 1 h ∂xâ 2 0 . Since h → 0 at the horizon, we need to require that ∂xâ 0 → 0 in order for the action to remain real. Therefore,â 0 takes a constant value on the horizon, which must be set to zero in order forâ to be a well-defined one-form at the horizon. This choice also implies that the boundary value ofâ 0 is precisely the (reduced) chemical potential. The value ofâ 0 is the only integration constant we need to fix in the IR, as the behavior of the various fields is otherwise determined by requiring regularity. It turns out, however, that imposing explicitly the leading-order constraint, which follow from expanding the equations of motion as series at u = 1, improves the stability of the code. We therefore require explicitly that our solution satisfies botĥ a 0 (x, 1) = 0 (44) as well as the three leading-order conditions from the equations of motion for the fields ψ,ẑ, andâ y , which are lengthy expressions and not reproduced here.
Note that the action (17) is covariant under the reflectionx → −x ifâ y is odd and the other fields even inx. 8 This property is reflected in the UV expansions (32), (33), (34), and (35) (if we set to zero the "trivial" constantā y ). Whenb = 0, there is also another symmetry: the action is covariant under the reflection if ψ − π/4 is odd and the other fields are even inx. There does not appear to be a specific reason solutions must obey this parity behavior. But, any inhomogeneous solutions we found respected the first symmetry and the solutions withb = 0 respected both symmetries when reflected with respect to some value ofx. By translation invariance we required that the first symmetry appears in reflections with respect tox = 0. Then the solutions withb = 0 have the second symmetry in reflections with respect tox =L/4.

A.2 Implementation of the pseudospectral method
Having fixed the boundary conditions, we implement the pseudospectral method based on a Fourier series in the x-direction and an expansion in the Chebyshev polynomials in the u-direction. We assume that the desired inhomogeneous solution respects the first parity symmetry with respect tox = 0. 9 Therefore, we choose the collocation points from a grid of N x + 1 evenly spaced points ranging fromx = 0 tox =L/2 in the x-direction and from the Gauss-Lobatto grid with N u points ranging from u = 0 to u = 1 in the u-direction. The problem is then reduced to finding the values of all the functions ψ,ẑ,â 0 , andâ y at the collocation points.
The solution is found by requiring that the boundary conditions and equations of motion are satisfied at the collocation points. Notice that the UV boundary conditions (36)- (40), except for the second one fixing the quark mass, and the IR boundary condition (44) directly fix the values of the functions at the boundary points, thus reducing the number of variables. For the other conditions, we need the values of the derivatives at collocation points. These are computed from the values of the functions in the pseudospectral approximation. In more detail, we require that the values of the functions match with the expansions at the collocation points. Here T j are Chebyshev polynomials, defined in the range [−1, 1]. This fixes the coefficients ψ j,k ,ẑ j,k ,â 0,j,k , andâ y,j,k in terms of the values of the functions. The derivatives are then computed from these expansions.
Using the expressions for the derivatives, the remaining boundary conditions and equations of motions can be turned into algebraic conditions for the values of the functions. The UV and IR boundary conditions are evaluated at the collocation points at u = 0 and at u = 1, respectively, and the four equations of motion are evaluated at the other collocation points with 0 < u < 1, so that the total number of conditions matches with the number of remaining variables. The resulting algebraic equations are then solved numerically by using the Newton method.

A.3 Initial configurations for the Newton method
The last step in the numerical method, i.e., solving the algebraic conditions discussed above, is a rather nontrivial task for our system. The difficulties arise because unlike in most AdS/CMT models where inhomogeneities have been studied in the literature, the square roots in our brane action give rise to nonanalytic behavior. As it turns out, the inhomogeneous solutions typically lie close to the branch points of the square roots and tiny deformations of these solutions can make the action complex. Consequently, the initial guess for the Newton method must be chosen very carefully, as otherwise the square root factors become complex at some intermediate step, which causes the method to diverge. Naturally, the main idea for constructing initial conditions is to start at low N x and N u , and iteratively increase the number of grid points using the previous solution as initial guess for the next step. In our case, however, the inhomogeneous solution may be absent at very low number of grid points (or more precisely, the solution has turned complex), which limits the range of N x and N u where the iterative procedure can be applied.
In order to overcome these difficulties, we used a separate code for generating the initial data, which also relied on a spectral method. We started with a very small grid (3-4 points in each direction) where the system could be solved by an initial Ansatz motivated by the fluctuation analysis [37,38]. The existence of the solution was guaranteed by adding a small random noise in the locations of the collocation points and by running the code in loop until it found a valid real root. The size of the grid was then gradually increased, keeping the random noise, until the grid was about ten points in each direction. After this, the obtained solution could be used as an initial configuration for the actual code, and grid size could be iteratively increased further without issues.

A.4 Convergence and accuracy
We tested the convergence of the code at a single set of parameter values (μ = 4, L = π/2,m = 0, andb = 0) up to N x = 38 and N u = 40. The expected roughly exponential convergence was found. At the largest grid size, the relative accuracy of a y was about 10 −9 , that of ψ andâ 0 was about 10 −10 , and that ofẑ was about 10 −11 . Most of the data in this article was computed using N x = 24 and N u = 26. At this grid size, the accuracy ofâ y was smaller than 10 −6 , that of ψ andâ 0 was smaller than 10 −7 , and that ofẑ was smaller than 10 −8 . This accuracy is more than enough for all analysis presented here.