Spectral weight and spatially modulated instabilities in holographic superfluids

Free fermions form a Fermi surface, which results in non-zero spectral weight at low energy and finite wavevector. In this work, we find similar features in holographic phases dual to strongly coupled quantum superfluid matter. At zero temperature, the phases we consider exhibit semi-local criticality in the IR and all the charge is carried by the scalar condensate outside the black hole horizon. Depending on the value taken by the IR critical exponents, we find Fermi surfaces in the transverse sector, Fermi shells in the longitudinal sector or no spectral weight at all. When there is non-zero transverse spectral weight, the IR can be subject to an instability at finite wavevector, the endpoint of which is likely a spatially modulated phase.


Introduction
Fermi liquid theory, a model explaining why a system of weakly-interacting fermions can exhibit properties similar to a non-interacting Fermi gas, is a mainspring in condensed matter physics. Weak coupling plays a crucial role, as the effective degrees of freedom are long-lived quasiparticles and can be thought of as dressed electrons. Although this theory describes the normal state of many metals, there also exist so-called non-Fermi liquids, characterized by an anomalous temperature scaling of the resistivity, specific heat, and other quantities while in their normal states [1]. These materials are of significant technological interest, as they include the high-T c cuprates that become superconducting at unexpectedly high temperatures [2]. Many of these non-Fermi liquid features are believed to trace back to the lack of long-lived quasiparticles. It is now a major endeavor to move beyond Fermi liquid theory and develop an understanding of these strongly interacting materials.
While standard perturbation theory techniques fail in the strong coupling regime, it is possible to investigate strongly interacting field theories using holography [3]. Despite their many scaling anomalies, non-Fermi liquids are seen to retain a Fermi surface momentum space distribution at low energies. Thus an important step in understanding these strongly correlated materials holographically is to first understand the contexts in which Fermi surfaces appear in holographic settings.
At strong coupling, the quantity that most explicitly marks the presence of a Fermi surface is called the low energy spectral weight: It is denoted by σ(k) to remind us that this is the low energy limit of the real part of the spatially-resolved conductivity appearing in Ohm's law (not to be confused with the optical conductivity where the k → 0 limit is taken first). The diagnostic power of the spectral weight is that it directly counts the number of degrees of freedom at a given frequency and momentum. By "degrees of freedom" we mean the density of states into which a given state can transition. This counting property is easily seen from the spectral decomposition of the retarded Green's function: In the expression above, m and n are eigenstates of the system. Notice that there are three momenta in the problem: k (the perturbing momentum), k (the momentum of state n) and k (the momentum of state m). The counting comes from the explicit energy conserving delta function, and the momentum conserving delta function resulting from the inner product. Further, the appearance of the current operator J guarantees that only charged degrees of freedom are counted. As we will discuss in the main body of the article, the low energy limit of the spectral weight is determined holographically by the IR behaviour of the holographic dual, that is by the near horizon, low temperature limit of the solution to the classical equations of the motion. Spectral weight in certain holographic theories has been studied previously, notably in gravity duals exhibiting a hyperscaling violating near horizon geometry (that is, geometry defined by the dynamical critical exponent z and the hyperscaling violating exponent θ) [4][5][6]. These geometries describe a flow between a UV relativistic fixed point (with z = 1) to an IR fixed point with emergent non-relativistic symmetry (z = 1). Hyperscaling violation corresponds to an effective reduction of the spatial dimensionality of the IR fixed point. For the hyperscaling violating geometries, it was found in [7] that spectral weight is always exponentially suppressed. The following year a particular limit of the hyperscaling violating geometry was studied, in which the quantity η ≡ −θ/z is held fixed as z → ∞ [8]. The phase of matter dual to this theory is called the semi-local quantum liquid for its finite spatial correlation length but infinite correlation time [9]. In both the hyperscaling violating and semi-local cases [7,8], the bulk matter content was given by the Einstein-Maxwell-dilaton theory [4,5], and all electric flux is sourced by the charged black hole horizon (since there is no explicit charged matter in the bulk). According to the holographic dictionary, this electric flux is in turn the source for the field theory current density. This is the object of interest, allowing us to calculate the spectral weight. For the semi-local quantum liquid, nonzero low energy spectral weight was found to exist in the range 0 < η < 2, for momenta less than a critical momentum: This is the signature of a (smeared out) Fermi surface. The infinite spectral weight at low momentum is an artifact of the zero temperature calculation. At nonzero temperature, this becomes finite (see [8] for further discussion). The dual charge density in the bulk gravity theory is behind the horizon, and we conclude that these charges are responsible for the low energy spectral weight observed holographically in the boundary field theory.
In this work, we are interested in the low energy spectral weight of holographic superfluids. A holographic model of superconductivity was first established in [10][11][12], and is reviewed in [3,13]. The approach taken in [12] to study these theories was to consider an effective bulk action of the form Here we have an explicit charged scalar in the bulk (outside the horizon) which can undergo spontaneous condensation below a critical temperature, and due to Schwinger pair production from the superconducting instability, all electric flux reaching the boundary is sourced by the charged condensate rather than the black hole horizon [13]. For a generic choice of a quadratic or quartic potential for the scalar potential V , the IR geometry is another copy of Anti de Sitter or a Lifshitz spacetime [14,15]. The general results of [7] indicate that in both cases, the low-energy spectral weight is exponentially suppressed as the exponent z takes a finite value. That is, the low-energy spectral weight does not seem to depend on the opening of a superconducting gap at all, contrary to weak coupling intuition. Rather, what appears to govern its absence are the scaling properties of the IR spacetime. The main goal of this paper is to explore what happens in the cases where z → +∞, for which [7,8] found non-zero low energy spectral weight. To be precise, weakly-coupled intuition dictates the spectral weight at finite momentum should vanish after the U(1) symmetry is spontaneously broken and the charged condensate has formed in the bulk, accompanying the expulsion of all the charge outside the horizon.
To explore this question, we will consider a slightly more general model [16,17]: This is an effective theory that captures the physics of superconductivity: the scalar has already undergone Bose-Einstein condensation, as evidenced by the massive photon. The scalar field φ should be thought of as the modulus of the original complex scalar field, while its phase has been integrated out. This model allows for several interesting features, depending on the IR behaviour of the action couplings V (φ), Z(φ) and W (φ). First, the IR geometries now include semi-local fixed points with z → ∞, even in the presence of a condensate. Second, these superfluid semi-local geometries can come with either charge in the bulk exclusively, or with both charge behind the horizon and in the bulk [16]. Here we will focus on the former case.
Our main result is that we do find nonzero low energy spectral weight at finite momentum in the boundary field theory, in spite of the fact that a condensate has formed and all the bulk charge now sits outside the horizon. The presence or absence of low energy spectral weight appears to depend on the scaling properties in the IR and not on the spontaneous breaking of the U(1) symmetry. This raises several interesting questions: 1) How should we interpret low energy spectral weight that exists independently of horizon charge? 2) What other degrees of freedom could this weight represent? 3) To what extent do bulk charge distribution properties represent those of the boundary charge?
In what follows, we will outline the calculation of the zero temperature low energy spectral weight for holographic superfluids in two decoupled sets of variables: the transverse and longitudinal channels. For the transverse channel we solve and decouple the bulk equations of motion directly, and for the longitudinal channel we demonstrate a scaling argument that allows us to infer the low-frequency dependence of the retarded Green's function without the need of a full bulk calculation. We end by a discussion of these results and suggest directions for future work.

Einstein-Maxwell-dilaton with massive vector
In this work we consider the Einstein-Maxwell-dilaton theory with a massive vector, which spontaneously breaks the U (1) gauge invariance of the theory The equations of motion of the theory in terms of a general Ansatz for the metric and matter fields are given in Appendix A. The low energy spectral weight is determined by the IR behaviour of the solution to the classical equations, so we focus only on solutions describing the near horizon region of spacetime at zero temperature. They have the following runaway behavior deep inside the bulk This is a necessary condition to depart from the scale-invariant solutions considered in [14,15]. It is known [16][17][18][19] that such a solution will emerge from the equations of motion of if we allow the coefficient functions to take the following form in the IR, loosely motivated by top-down string theory realization of the model (2.1): 1 For the rest of the analysis we are free to set Z 0 = 1, as this is simply a rescaling of the gauge field. More precisely, the scaling solutions will take the general form The logarithmically running scalar manifestly realizes our condition (2.2). The metric and gauge field are covariant under rigid rescalings of the coordinates t → λ z t, (r, x, y) → λ(r, x, y). Upon turning on a small temperature, the entropy density is seen to scale as s ∼ T (2−θ)/z , which exemplifies how the spatial dimensionality of the boundary field theory has been effectively modified to d ef f = 2 − θ and characterizes hyperscaling violation in the IR field theory. In the bulk, the critical exponent ζ encodes departure from scale invariance for the gauge field. From the point of view of the dual theory, it governs the IR scaling of the electric conductivity [17,26]. It also contributes anomalously to the dimension of the charge density operator in the IR [17,26,27] (though one should note that when the current is not conserved, its dimension is not protected and the usual arguments preventing an anomalous dimension [28] do not apply). The parameters of the solution (2.4) are fixed in terms of parameters in the action, so (θ, z, ζ, L, κ) are actually functions of ( , γ, δ, W 0 , V 0 ). In this paper, we will not study the general class of scaling solutions (2.4) (for details see [16,17]) but instead focus on their semi-local limit z → +∞, θ → +∞, ζ → +∞ with η = −θ/z andζ = ζ/z kept finite. In the remainder of the paper, we drop the tilde.
Two cases may be distinguished, depending on whether the bulk charge sits entirely outside the horizon or also behind it. In this work, we will focus on the first possibility and refer to [16] for more details on the second.
The η geometries are exact solutions to the background Einstein equations provided the constraint = γ − δ is enforced [16,17]: Note that these solutions do not exist for all values of γ, δ and W 0 . Instead, choosing two fixes the third. This places constraints on the effective actions which admit these solutions. We have also chosen to set the IR 'AdS radius' to 1 by fixing V 0 . This is of course not necessary (and actually not desirable to construct the full flow to a UV AdS 4 ).
It is important that we only consider a physically-consistent parameter space. The constraints we have to take into account are: the null energy condition (NEC), the reality of all metric/scalar/gauge field coefficients, negativity of V 0 and positivity of W 0 (since this is a mass term and the charge squared of the scalar [16]). We also require that the specific heat is positive. According to our convention, V 0 < 0 because the V (Φ) term in the action replaces the cosmological constant term −2Λ, and Λ < 0 for Anti-de Sitter space. Furthermore, from the metric (2.5) we see that in order to have a well-defined IR we must have η(η + 2) > 0. This constraint is echoed in the NEC, which requires that the gravitational field generated by the stress tensor is attractive. Specifically, the NEC gives [29]: Combining all of our constraints, we find From these conditions, it follows that the IR is always r → +∞ in these coordinates and the electric flux, which scales as r −ζ−η , always vanishes there.
Finally, we should also be careful that there exists the right number of irrelevant deformations around the IR geometry. These correspond to the scaling dimensions of operators. The method to compute the radial deformations and the identification of the IR scaling dimensions is described in detail in [16,17]. The deformations all sum to 1 + η, which is the dimension of the free energy density in the IR. There are three pairs of modes. Two of them are simply (1 + η, 0). The other pair is the one of interest and should be used to connect the IR geometry to the UV AdS geometry. It reads We can see that β + > 0, so this is always a relevant mode. We need to impose that β − is irrelevant, which leads to the reduced parameter space 11) This is shown in Figure 1.
These solutions are called cohesive as the electric flux they source vanishes in the IR: the horizon is neutral and all the charge is generated by the charged condensate in the bulk.

Computation of the spectral weight
We now perturb the background fields and determine the linearized equations of motion. If we choose the plane wave perturbations to be in the x-direction, as δX(r, x, t) = δX(r)e i(kx−ωt) , then the perturbed modes naturally decouple into two categories: those even or odd under the parity transformation y → −y. The set of odd modes are called transverse {δA y , δg xy , δg yt , δg yr } and the even modes are called longitudinal {δA t , δA x , δA r , δg tt , δg rt , δg xt , δg xr , δg xx , δg yy , δg rr , δΦ}. Not all of the perturbations mentioned above are independent, and we must combine them to form gauge invariant variables. The gauge group here is just diffeomorphism invariance, and so gauge invariant variables are those that are invariant under a coordinate transformation X α → X α + ξ α . Under such a transformation, the metric changes according to the Lie derivative Thus for a gauge invariant variable, the Lie derivative must vanish. Note that the choice of variables is not unique, though some choices are wiser than others as a tool for decoupling the equations of motion. Since there are four degrees of freedom in ξ α , we can pick a gauge such that four metric perturbations are zero. A popular choice is the so-called radial gauge δg µr = 0. However, it will turn out that postponing this gauge choice will allow us to more easily decouple the perturbed equations of motion in the transverse channel. The vector potential also transforms according to the Lie derivative: The vector ξ is contracted, so this transformation gives us the freedom to set only one of the variables δA µ equal to zero. A further discussion of Lie derivatives and gauge invariant variables is given in the Appendix.

Transverse channel
In the transverse channel the modes are:

(3.4)
Note that it is sometimes more useful to work with perturbed variables with one index raised.
The resulting equations of motion (prior to substituting in the background quantities) are The perturbations we have introduced are not gauge invariant variables, and so our first task is to find linear combinations of the fields {δA y , δg y x , δg y t , δg y r } that are gauge invariant. We introduce the vector field ξ = e −i(kx−ωt) ξ y (r)∂ y (3.9) and calculate the Lie derivative with respect to this field. We have chosen ξ so that it has one nonzero component, the y-component. This is because we are working in the transverse channel, where the modes are odd under the parity transformation y → −y.
An explicit calculation is given in the appendix. The answers are: In general, beyond finding vanishing combinations of these quantities, one must also take care that no longitudinal modes { δA t , δg rt , δg yy , ...} are generated in this process.
In this case, that issue does not arise. We choose the combinations: The factors of r included above allow us to decouple the equations of motion using the technique of "master variables" [30]. In this approach, one finds a linear combination of the invariant variables ψ 1 and ψ 2 that will automatically decouple the equations. A brief discussion of how to find such variables is given in the Appendix (though no general algorithm is yet known). The answer is: In terms of these fields, we obtain two decoupled equations of motion: with 2ν ± = 5 + 2η + η 2 + 4k 2 ± 4X. (3.15) It is quite remarkable that the transverse perturbation equations can still be decoupled after the introduction of the mass term for the vector. As we will see below, we will not be so lucky in the longitudinal sector. The decoupled equations are solved by Bessel functions We must choose only the ingoing modes at the horizon (r → ∞), since we know that matter falls into the black hole and does not come out. The solution satisfying these infalling boundary conditions is This is the near-horizon solution. We are ultimately interested in the retarded Green's function in the UV, G R (ω, k). The imaginary part is the so-called spectral weight of a given operator, which counts the number of degrees of freedom that overlap with that operator, ImG R O A O A (ω, k). In general, to achieve this one must solve the equations of motion throughout the entire bulk, rather than just in the near horizon limit like we did. However, as shown in [31], for the geometries we are considering it is possible to map the IR Green's function G R (ω, k) directly onto the UV one G R (ω, k) through a matching procedure, provided we are only interested in the low frequency behaviour of the UV Green's function. The result is: where the d s are real constants and the sum runs over the fields involved. In our case this is just G R + and G R − . Returning to the IR solution, the Green's function is obtained by taking the limit of our near horizon solution as r → 0. This corresponds to a solution far from the horizon: We are now in a position to access the small ω dependence of the UV Green's function It is the imaginary part of this quantity that gives the low energy limit of the spectral weight, which is the quantity we wish to compute. Since the constants in are all real, and because the ν − exponent dominates over the ν + exponent at low energies (cf. (3.15)), we obtain ImG R ∝ ω 2ν − . (3.22) We note at this point that ν ± are both real in the allowed parameter space (2.11). Thus we find the following low energy spectral weight: This is a similar conclusion to the one reached in [8]. At zero temperature, this system enjoys low energy spectral weight over a finite range of momenta. This result runs contrary to weak coupling intuition, where we do not expect low energy degrees of freedom at finite momentum when the charges manifestly form a condensate. At strong coupling, gapless degrees of freedom can survive thanks to the scaling geometry in the IR. Notice, however, that k vanishes for the special value η = 2 (independent of the value of ζ). For η > 2, there is never any low-energy spectral weight. It is worth recalling that without a condensate, k also vanishes for η = 2 [8]. The full parameter space where the transverse spectral weight diverges at low frequency is depicted in Figure 1. We can also see that ζ > −2 in order for the low energy transverse spectral weight to diverge.
There is an alternate method of obtaining the ω dependence of the Green's function without the necessity of decoupling the equations of motion. We will demonstrate this technique here for the transverse channel, and will then take advantage of it in the computationally complicated longitudinal channel.
In this method, we endow all perturbations with a scaling behavior: δA y (r) = a 0 r a 1 , δg ty (r) = t 0 r t 1 , δg xy (r) = x 0 r x 1 , δg ry (r) = 0. (3.25) Upon substituting this Ansatz into the four perturbed equations of motion, we notice that by fixing t 1 and x 1 in terms of a 1 , 2 the equations take the schematic form A(a 0 , t 0 , x 0 ) + B(a 0 , t 0 , x 0 )ω 2 r 2 . Namely, all the r and ω dependence is relegated in the second term. This means we can set ω = 0 and simply scale r to reinstate it. This is in agreement to the exact solution (3.17) we found using the master variables. By manipulating the equations, we can solve for two of (a 0 , x 0 , t 0 ) in terms of a third. Finally, the remaining equation can be solved to get an expression for a 1 in terms of η and ζ, yielding 3 a 1 = 1 2 2ζ + η + 1 ± 5 + η 2 + 2η + 4k 2 ± 4X , (3.26) where X is defined above in (3.13). Note that the radical is the same as that obtained in (3.15), but here we avoided using the full machinery of gauge invariant variables and decoupling the equations of motion.

Longitudinal channel
In the longitudinal sector, we could not succeed in decoupling the field equations in terms of gauge invariant variables. We could not even express the equations solely in terms of gauge invariant variables. This is likely due to the fact that a δa t perturbation is turned on, which directly affects the right-hand side of the t component of Maxwell equations via the mass term for A µ . Nonetheless, we can still work out the scaling dimensions of the IR operators according to the method outlined at the end of the previous section. In more details, the constraints coming from Einstein equations can be used to solve for δa x and δg tx . Upon replacing and substituting for the longitudinal fluctuations an Ansatz similar to (3.25), the remaining equations all take the form A(a 0 , t 0 , x 0 ) + B(a 0 , t 0 , x 0 )ω 2 r 2 , as in the transverse sector. Applying similar manipulations to this set of equations, after setting ω = 0, we find the following IR scaling dimensions: We can check that the dominant exponent at low frequencies is always ν − . Setting W 0 = 0 (or equivalently ζ = −η), we recover the result reported in [8]. In this limit, [8] found no low energy longitudinal spectral weight, that is 2ν − − 1 is always positive. There are several interesting differences when W 0 = 0. First, while ν 0 and ν + are always real in the allowed parameter space (2.11), there is a region where ν − can be complex for a range of k: (3.28) and k − < k < k + , where (ν − (k ± )) 2 = 0. Effectively, this restricts ζ < 0, see Figure  2. This suggests the homogeneous superfluid phase is unstable and that the endpoint of this instability is a spatially modulated, superfluid phase. To confirm its presence, we would need to construct the superfluid geometry at finite temperature and look for normalizable modes at non-zero k, along the lines of [32]. Spatially modulated holographic superfluids have been constructed before in various other setups [33][34][35][36][37][38][39]. In this region of parameter space, there is always a corresponding Fermi surface in the transverse sector.
Also in Figure 2, we show in the brown region the parameter space (η, ζ) where 2ν − − 1 is negative for k − < k < k + (where 2ν − (k ± ) − 1 = 0). When k < k − or k > k + , then 2ν − − 1 > 0. This suggests we do not have a Fermi surface and corresponding Fermi sea for k ≤ k F , but rather a Fermi shell for k − < k < k + . This shell may result from two nested Fermi surfaces, one of 'particles' and one of 'holes'. When the longitudinal sector exhibits this shell, there is always a corresponding Fermi surface in the transverse sector. Let us note here that Fermi shells have also appeared in top-down truncations, both in N = 4 SYM [40] and ABJM theory [41]. There, the situation is clearer as the field content of the dual field theory is known and there are indeed two

Discussion
We find that holographic superfluids with a semi-locally critical IR can exhibit low energy spectral weight over a finite range of momenta in the transverse correlators, which is a signature of a (smeared) Fermi surface. This result is surprising, as the bulk charge density that is thought to be responsible for this spectral weight manifestly forms a condensate in our model. The formation of the condensate is not accompanied by a gapping out of the low energy degrees of freedom.
In the longitudinal sector, we find that there is no low-energy spectral weight in a large part of the parameter space. Interestingly, we do find a region where there exists a Fermi 'shell', which suggests the presence of two Fermi surfaces, one of 'holes' and one of 'electrons'. In yet another region of the parameter space, we find that the scaling dimension of the least irrelevant IR operator becomes complex for non-zero values of the wavevector, indicating an instability, the endpoint of which is likely to be a spatially modulated superfluid phase. Obviously, it would be very interesting to construct these new ground states. Both of these features are always accompanied by non-zero spectral weight in the transverse sector, although the converse is not true.
While we have focussed on the case where all of the bulk charge sits outside the horizon, superfluid η geometries with charged horizons can also be constructed [16]. Their spectral weight is computed by same scaling arguments as developed here, and would be characterized by the same ν(k) exponents as in [8]. The reason for this being that the effects of the mass term for the vector are subleading in the IR and so do not affect the IR Green's function at low frequency. Thus the same distribution of spectral weight as in [8] would be found: transverse spectral weight provided 0 < η < 2, but no longitudinal spectral weight.
Our results raise some interesting questions about spectral weight in holographic theories, such as: a) To what extent do charge distribution properties in the gravity theory represent those of the field theory? More precisely, it would be desirable to understand more clearly the relation between bulk charge carried by the horizon or a condensate, and the normal and superfluid charge densities on the boundary. Preliminary work in this direction appeared in [42]. b) How should one interpret spectral weight at finite wavevector when all the charge in the bulk is manifestly carried by a condensate? and c) Why does the presence of low energy spectral weight seem to depend more upon the background geometry than matter content?
To work toward answering these questions, it may be important to recall that the spectral weight found in the holographic superfluid case vanishes for specific values of the background parameters of the theory, namely η = 2. Interestingly, a recent work has uncovered an example of a theory such that, after being dimensionally reduced from a higher dimensional Supergravity theory, the Fermi surface present in fermionic correlators vanishes when the parameters in the holographic action take their Supergravity values [25]. Moreover these top-down constructions have the advantage over bottom-up ones in that they have a well-identified field theory dual.
Our results also connect with work on the presence of holographic Friedel oscillations in the spatially resolved static susceptibility χ(k), either for Reissner -Nordstrom AdS [43] or the U (1) 4 black hole of [44] with three of the charges equal and the fourth zero. In the former case, exponentially damped Friedel oscillations were uncovered at long distances, apparently at the wavevector where Re[ν − (k)] = 0 (with the expressions appropriate for an AdS×R 2 IR geometry). [45] found that, on the contrary, these oscillations do not survive at long distances in the U (1) 4 black hole, which happens to have an η = 1 IR. This is consistent with the results in [8] that there is no low energy longitudinal spectral weight. It would be worthwhile to compute the static susceptibility in our superfluid geometries and work out whether Friedel oscillations are present in the corner of parameter space where we do find longitudinal spectral weight.

A Field Ansatz and equations of motion
We begin with the general background fields ansatz: The background equations following from this ansatz are the Maxwell equation the scalar equation and the Einstein equations There is a background conserved quantity which gives the Smarr law when evaluated at the boundary and at the horizon.

(B.3)
To get rid of the factors r −η , it is common to work in one raised index, so that [L ξ g] y t = −iωξ y , [L ξ g] y r = ξ y , [L ξ g] y x = ikξ y .

(C.3)
The master variables will be some linear combination of the ψ's. A general combination would be φ ± = f 1± (r)ψ 1 (r) + f 2± (r)ψ 2 (r) + f 3± (r)ψ 1 (r) + f 4± (r)ψ 2 (r), (C.4) but we would like a simpler choice. A good practice is to only include two of the terms that occur in both coupled equations. In the present case, this rules out ψ 1 and ψ 2 , so we have φ ± = f 1± (r)ψ 1 (r) + f 2± (r)ψ 2 (r). (C.5) Now it remains to solve for the functions f 1± (r) and f 2± (r) that satisfy φ ± ∝ φ ± . We will see how this condition is also used to fix the exponents N 1 and N 2 . To that end, let's take a look at φ ± : φ ± = f 1± (r)ψ 1 (r)+2f 1± (r)ψ 1 (r)+f 1± (r)ψ 1 (r)+f 2± (r)ψ 2 (r)+2f 2± (r)ψ 2 (r)+f 2± (r)ψ 2 (r) (C.6) We can then use (C.2) and (C.3) two get rid of the second derivatives in favor of lower order terms. Schematically, then, φ ± becomes φ ± = g 1± (r)ψ 1 (r) + g 2± (r)ψ 1 (r) + g 3± (r)ψ 2 (r) + g 4± (r)ψ 2 (r). (C.7) To satisfy the condition φ ± ∝ φ ± , we see immediately that the coefficient functions of ψ 1 (r) and ψ 2 (r) should be set to zero. This yields two differential equations, one in f 1± (r) and the other in f 2± (r), and the exponents N 1 and N 2 are chosen so that these two functions are constants. From there it is easy to determine what these constants need to be in order to achieve φ ± ∝ φ ± . As we noted earlier, it is not always possible to find simple solutions to these differential equations, and as yet there are no criteria to determine from the beginning which coupled equations can be decoupled using master variables and which cannot.