Density wave instabilities of tilted fermionic dipoles in a multilayer geometry

We consider the density wave instability of fermionic dipoles aligned by an external field, and moving in equidistant layers at zero temperature. Using a conserving Hartree-Fock approximation, we show that correlations between dipoles in different layers significantly decrease the critical coupling strength for the formation of density waves when the distance between the layers is comparable to the inter-particle distance within each layer. This effect, which is strongest when the dipoles are oriented perpendicular to the planes, causes the density waves in neighboring layers to be in-phase for all orientations of the dipoles. We furthermore demonstrate that the effects of the interlayer interaction can be understood from a classical model. Finally, we show that the interlayer correlations are important for experimentally relevant dipolar molecules, including the chemically stable $^{23}$Na$^{40}$K and $^{40}$K$^{133}$Cs, where the density wave regime is within experimental reach.


Introduction
The trapping and cooling of molecules in their rotational and vibrational ground state is a new research direction within the field of ultracold atomic and molecular physics [1,2,3,4,5,6,7,8,9,10]. In contrast to the short-range isotropic interactions in typical cold atomic gases, molecules provide anisotropic potentials that typically have a long range dipolar part. This opens up a host of possibilities for exploring interesting physics [11,12,13] and also chemical reaction dynamics at low temperatures [14,15].
Chemical reaction losses can be large in three-dimensional (3D) samples [6] due to the attractive head-to-tail interaction of dipolar molecules.
However, recent experiments using optical lattices [9,10] have shown that a low-dimensional confinement of the system can effectively suppress the loss. Interesting many-body phases have been proposed in such settings, including s-and p-wave superfluid states in single-and bilayer setups [16,17,18,19,20,21,22,23], density-waves in homogeneous [24,25,26,27,28,29] and lattice systems [30,31,32], and non-trivial Fermi liquid behavior [33,34,35]. The long-range dipolar forces also opens up for a very rich spectrum of few-body bound state physics in 1D [36,37,38,39] and 2D setups [40,41,42,43,44]. The dipoles reside in layers in the xy-plane separated by the distance d. The angle θ, between the dipoles and the normal to the layers is in the xz-plane. The interaction between dipoles in the same layer, in adjacent layers, and in layers separated by a distance 2d is indicated by V 0 , V 1 , and V 2 respectively. (b) The phase boundary between the normal (left) and the striped phase (right) in the (g, θ) plane for strictly 2D layers with w = 0: a single layer (--), two layers separated by dk 0 F = 0.5 (----) and dk 0 F = 1 (-· -), and three layers separated by dk 0 F = 0.5 (----) and dk 0 F = 1 (-· -). The · · · · · · line gives the RPA (Random Phase Approximation) result for a single layer for comparison. The case of a single quasi-2D layer with wk 0 F = 0.1 is plotted by ----.
We consider the density wave (stripe) instabilities of fermionic dipoles at zero temperature in 2D layers. The dipole moments are aligned by an external field, and they are moving in equidistant layers as illustrated in figure 1(a). Calculating the density-density response function within a conserving Hartree-Fock approximation, we find the following. (I) The presence of several layers can decrease the critical coupling strength for stripe formation significantly. The effect is strongest for the dipoles oriented perpendicular to the planes whereas it vanishes when the dipoles are aligned in the planes, see figure 1(b). (II) The mechanism for this decrease is the formation of stripes on top of each other in adjacent layers. This in-phase stripe pattern is always energetically favorable, independent of the orientation of the dipole moment which is somewhat surprising. It can be understood from a purely classical calculation which also explains the angle dependence of the effect. (III)  The decrease in the critical coupling strength for stripe formation due to interlayer correlations is significant for experimentally relevant systems consisting of 7 Li 40 K, 23 Na 40 K, 40 K 133 Cs, 6 Li 87 Rb and 6 Li 133 Cs molecules. It is less important for 40 K 87 Rb and 6 Li 23 Na molecules which have smaller dipole moments.

System
The system consists of identical fermionic dipoles of mass m moving in equidistant layers separated by a distance d. Their dipole moment p is aligned forming the angle θ with respect to the normal of the layers (z-axis) with their projection onto the planes defining the x-axis, see figure 1. The layers are formed by a deep 1D optical lattice so that the dipoles in layer l (l an integer) reside in the lowest Wannier state in the z-direction. This is to a good approximation a Gaussian ϕ l (z) = exp[−(z − ld) 2 /2w 2 ]π −1/4 w −1/2 , where w is the width of the layer which is centered at ld. We neglect any trapping potential in the xy-plane so that the transverse states are labelled by the momentum k = (k x , k y ) (we use units where = k B = 1). In this basis the grand canonical Hamiltonian readŝ where we have absorbed the harmonic oscillator energy of the z-direction in the chemical potential µ. Here,ĉ k,l removes a dipole in layer l with momentum k. The interaction between two dipoles separated by r is where θ r is the angle between r and the dipole moment p, and D 2 = p 2 /4πε 0 for electric dipoles. The effective interaction V l,l (q) between dipoles in layer l and l is obtained by integrating the interaction V 3D ( r) over the Gaussians ϕ l (z)ϕ l (z) combined with a 2D Fourier transform. This yields [46] V for the intralayer interaction where P 2 (x) = (3x 2 − 1)/2 is the second Legendre polynomial, while The constant term in equation (3) corresponds to a δ(r) interaction which plays no role since we consider identical fermions. The real part of the interlayer interaction is only dependent on the difference l in layer numbers and is given by [47,48] V l (q) = −2πD 2 ξ(θ, ϕ)qe −dlq (5) for w d. This approximation deviates less than 10% from the exact expression for w ≤ d/5. In equation (5), we have only given the real part of the potential since the imaginary part is zero for momenta along the y-direction which are the ones relevant for stripe formation, as discussed in section 3.
The strength of the interaction is parametrized by the dimensionless number where k 0 F = √ 4πn 2D is defined from the 2D density. Likewise, a dimensionless measure for the layer separation is given by dk 0 F . The ratio of the layer distance to the typical inter-particle distance in a layer is dk 0 F / √ 4π.

Linear response and the Hartree-Fock approximation
To analyze the instabilities of the homogeneous phase towards the formation of stripes, we consider the retarded density-density response function χ R ij (r − r , t − t ) = −iθ(t − t ) ρ i (r, t),ρ j (r , t ) where (i, j) denotes the layers and r = (x, y). The density operator for layer i isρ i (r) =ψ † i (r)ψ i (r) withψ i (r) the field operator for the dipoles in layer i. An instability toward the formation of a density wave with wave number q shows up as a zero frequency (ω = 0) pole of the Fourier transformed response function χ R ij (q, ω). We calculate the retarded response function using diagrammatic perturbation theory in Matsubara space χ ij (q c , iω n ) with ω n = 2nπT (n integer) a bosonic Matsubara frequency. The retarded function χ R ij (q c , ω) is then obtained by analytical continuation iω n → ω + i0 + in the usual way [49]. As illustrated in figure 3, χ ij (q) with q = (q, iω n ) can be written as where Π i (k, q) = G i (k + q)G i (k) is the particle-hole propagator with G i (k) the single particle Green's function for the dipoles. The scattering matrix Γ ij (k, k , q), which describes a particle with momentum k + q scattering on a hole with momentum k to produce a particle with momentum k + q and a hole with momentum k , obeys the Bethe-Salpeter equation Γ ij (k, k , q) = I ij (k, k , q) + l,k I il (k, k , q)Π l (k , q)Γ lj (k , k , q) (8) as shown in figure 3. Here I il (k, k , q) includes all scattering processes which are irreducible with respect to the particle-hole propagator Π l (k, q).
To proceed, we apply the self-consistent Hartree-Fock approximation illustrated in figure 4. All Green's functions are interacting in this approximation and the vertex I il (k, k , q) is given by the lowest order direct and exchange interactions. Writing the Green's function as G −1 l (k) = ik n − ε lk with ik n = (2n + 1)πT a fermionic Matsubara frequency and ε lk = k 2 /2m − µ − Σ l (k), we have the usual Hartree-Fock expression for the self-energy Here f lk = [exp(βε lk ) + 1] −1 is the occupation of the k momentum state in layer l.
In the calculations, we take β = 1/k B T → ∞ appropriate for the zero temperature case. It follows from (5) that the Hartree shift V l,j (0) for the energy of a dipole in layer l due to the interaction with the dipoles in layer j = l vanishes for thin layers, i.e. for w d. This is a consequence of the fact that the dipole-dipole interaction integrates to zero over a plane as shown in Appendix C. In this paper, we keep w d and since the density is the same in all layers, the self-energy and the chemical potential are both layer independent. We calculate the single particle Green's functions numerically, obtaining self-consistency through an iterative procedure. The scattering matrix Γ ij (k, k , q) is then determined from equation (8) using Note that exchange is only included for interactions within the same layer as dipoles in different layers are distinguishable. Diagrammatically, this approximation corresponds to the summation of bubbles containing intralayer ladder interactions, which are connected to each other by inter-and intralayer interactions.
Since the irreducible interaction vertex I i,j is independent of frequency in this approximation, the particle-hole scattering matrix is independent of the internal frequencies ik n , ik n and we can perform these frequency summations in equation (8), which then simplifies to Equation (11) corresponds to a series of inhomogeneous Fredholm equations of the second kind in the first variable. From equation (7) it follows that apart from the poles of Π describing the particle-hole continuum, the density-density response function χ and the particle-hole scattering matrix Γ share the same poles describing collective modes. Thus, we determine the critical coupling strength by searching for a zero frequency pole at a non-zero momentum q of the matrix Γ ij (k, k , q), which is analytically continued to real frequency by iq n → ω + i0 + . When θ > 0, the anisotropic interaction favors dipoles that are aligned along the x-axis. This corresponds to a density wave with a wave vector pointing in the perpendicular direction ϕ c = π/2. Since the density wave is formed by particle-hole excitations, we expect the length of the wave vector to be 2k F (ϕ c ) to minimize the kinetic energy cost. This indeed follows directly from RPA calculations of the densitydensity response function [24,25], and we expect it to hold even when exchange correlations are included.
The self-consistent Hartree-Fock calculation of the response function is demanding numerically, and we describe in Appendix A how it is implemented numerically. The payoff is that the approximation is conserving in the sense of Kadanoff-Baym [50]. We shall furthermore see that the exchange correlations have large effects on the critical coupling strength for the formation of stripes.

Numerical results
We now present numerical results for the critical coupling strength g c as a function of dipole angle θ, layer separation dk 0 F , and layer thickness wk 0 F for fixed momentum in the direction ϕ c = π/2 with magnitude q = 2k F (ϕ c ). We shall for concreteness consider the cases of one, two, and three layers.

Single layer
We first focus on the case of a single layer with vanishing thickness, i.e. w = 0. In figure 1(b), we plot as a solid line the phase boundary between the normal phase (left) and the striped phase (right) for a single layer. The boundary has an intriguing non-monotonous behavior with a maximum critical coupling strength g c for θ π/4. For comparison we also plot as a dotted line the result of a RPA calculation using the interacting Green's functions [24,25].
For small dipole angles θ, the RPA result underestimates the critical coupling strength significantly as compared to the conserving HF approximation. This demonstrates that exchange correlations suppress the formation of stripes. For larger angles, the shape of the phase boundary differs qualitatively from that obtained from the RPA calculation which predicts g s ∝ cos(θ) 2 (neglecting the effects of Fermi surface deformation).
The dependency on θ for the RPA result (green, dotted line of figure 1(b)) can be understood purely from the fact that the repulsive part of the interaction decreases as cos 2 θ as the dipoles are tilted towards the layer. For small θ, the exchange correlations suppress stripe formation. As θ crosses the "magic angle" cos −1 (1/ √ 3), the spatial (or momentum) average of the interaction goes from being repulsive to being attractive. Since exchange correlations enter through an average over the momentum transfer in the term V (q) − V (k − k ), see (10) and the k -sum in (8), this means that the effects of the exchange term vanishes right at the magic angle as can be seen from figure 1(b). For larger θ, exchange correlations enhance the stripe instability and for θ → π/2 the instability is entirely driven by the exchange term since the direct interaction vanishes. This is the qualitative origin of the maximum in the critical coupling strength for an angle close to cos −1 (1/ √ 3). Our result for the phase boundary agrees within 7% with that of ref. [28], and the critical coupling strength for θ = 0 agrees within 5% of that reported in ref. [27]. As explained in Appendix A, we have taken care to use many k-space points in the integration around the singular points in (11), and we estimate our results to be numerically accurate within 1%. On the other hand, our results differ substantially from those obtained using a self-consistently determined local field factor [29].
We also plot in figure 1(b) the phase boundary for a quasi-2D single layer system where the finite depth of the 1D optical lattice gives a width of wk 0 F = 0.1 for the Gaussian transverse wavefunctions. This softening of the z-direction reduces the repulsive part of the effective dipole-dipole interaction (3), and as a result the critical coupling strength for stripe formation is increased by up to 18%. For larger wk 0 F , one has to take into account higher states in the z-direction as investigated for the case of θ = 0 in Ref. [27].
Note that we for simplicity have not included the region of p-wave superfluidity for θ π/4 [16] in the phase diagram, nor the region of collapse due to a negative compressibility for large angles and coupling strengths [16,25,29].

Several layers
In figure 1(b), we also plot the phase boundary in the case of two and three layers separated by the distances dk 0 F = 0.5 and dk 0 F = 1. This illustrates a main result of this paper: the presence of neighbouring layers reduces the critical coupling strength for stripe formation significantly when the layer distance d is comparable to the distance between particles within each layers. As expected, the effect decreases with increasing layer separation as follows directly from the exponential decay of the interlayer interaction (5). This is illustrated further in figure 5(a) where the critical coupling strength for the bi-layer case is compared to that of the single layer case as a function of layer separation and dipole angle. For small layer separation, the critical coupling strength can in fact be shown to scale as 1/N for N 1 layers, since the exchange correlations within each layer can be neglected in this limit, so that the problem reduces to that of dipoles with N internal degrees of freedom moving in a single layer [26]. From figure 5(a) we furthermore conclude that the effects of neighboring layers is strongest for small dipole angles. This effect has a simple classical interpretation as we shall demonstrate below. Figure 5(b) shows the critical coupling strength for the bi-layer case for w = 0. It increases with increasing layer separation dk 0 F whereas it has an interesting angular dependence as a result of the interplay between the angular dependence of the interlayer and intralayer interaction. Note that for well-separated layers with w d, a non-zero layer thickness only changes the intra-layer interaction, and it therefore affects the critical coupling strength in a way similar to that of a single layer.
It should be noted that s-wave interlayer superfluidity is likely to occur in the multilayer system. This has been discussed for perpendicular dipoles in for both bilayer [19,20,21] and multilayer cases (restricted to nearest-neighbour pairing) [62].
For θ = 0, the interlayer superfluidity is expected for any value of g when densitywaves are ignored. At sufficiently large g, we will likely have a competition of the two phases and a model containing both is necessary to infer which phase dominates or whether there can be coexistence. The superfluid could have p-wave symmetry for larger angles θ which may be more favorable to coexistence. This will be explored in future studies.

Correlations between stripes in different layers
To examine further the effects of neighboring layers, we now analyze the ω = 0 collective mode and the associated density oscillations at the critical coupling strength g c . We first focus on the case of two layers. The density oscillations induced by an external perturbation l d 2 rV ext l (r, t)ρ l (r) are within linear response given as At the critical coupling strength g c , one eigenvalue of the density response matrix in equation (12) diverges. We find that the mode which first diverges always is symmetric in the layer index l, independent of the dipole angle, except for θ = π/2 where the modes are degenerate. Close to the critical coupling 1 − g/g c 1, the density-density response matrix has the pole structure as we demonstrate in detail in Appendix B. Here, χ 0 is the Lindhard function including Hartree-Fock shifts of the single particle energies. Thus, the density instability in the two layers is in-phase and the stripes will be on top of each other as illustrated in figure (1). The mode which is anti-symmetric in the layer index becomes unstable at a larger coupling strength. This can be understood as a splitting of the eigen-mode for a single layer into a symmetric and anti-symmetric mode. The same result holds for more than two layers: The mode with the stripes in phase always becomes unstable first, irrespective of the value of θ, except for θ = π/2 where the modes are degenerate. This is illustrated in figure 6(a) which depicts the critical coupling strength for the even and odd modes for the cases of two and three layers. The modes with higher g c were calculated using the self-energy at the value for the lowest modes for simplicity. This changes the obtained values slightly but retains the relative ordering of the modes. For the case of three layers, there is an additional mode with no density fluctuations in the middle layer. It goes unstable for a coupling strength in between the values for the even and odd modes. This agrees with the results found within the RPA [26] ‡. ‡ The odd mode and the mode with no density modulations in the middle layer was mistakenly swapped around in the text of ref. [26].

Classical model
We now demonstrate that the effects of neighboring layers can be understood from simple classical considerations. First, it is easy to show, see (C.2), that the classical interaction energy of a single dipole with an infinite layer of dipoles with homogenous density is zero. Second, we analyze the interaction between stripes by calculating the classical interaction energy between a single dipole and stripes of increased/decreased density of dipoles in a layer separated by a distance d. As explained in detail in Appendix C, a straightforward calculation gives that this interaction energy is where the − is for the in-phase case where one stripe in the plane is directly above the single dipole, and + is for the anti-phase case where the projection of the dipole onto the plane is in between two stripes. The linear density of dipoles within a stripe is given by γ, and λ c is the distance between the stripes within the layer. Equation (14) is plotted in figure 6(b). This demonstrates that the in-phase/out-of-phase configuration of the stripes has a negative/positive interaction energy thereby explaining why they decrease/increase the critical coupling strength g c for stripe formation as compared to the single layer case. The cos 2 θ dependence in (14) furthermore explains why the effect decreases for increasing angle, with the modes becoming degenerate for θ = π/2 as seen in figure 6(a). It is reassuring that our full quantum mechanical calculation which is rather numerically involved, agrees with this classical analysis.

Experimental realisations
We now examine the importance of the interlayer correlations discussed above for typical experiments. In an experiment where the planar confinement is caused by a 1D optical lattice formed by counter propagation lasers with wavelength λ, the layer distance will be d = λ/2. Figure 2 shows the critical dipole strength times the square root of the mass, √ mp, as a function of density for one, two, and three layers and the dipoles perpendicular to the planes for a typical wavelength of λ = 1064 nm. For comparison, in [9] the JILA group reports a peak density √ n 2D = 0.58 µm −1 or dk 0 F = 1.1 in an experiment with 40 K 87 Rb molecules in the rotational and vibrational ground state. In the figure, the horisontal lines are the permanent electrical dipole moment times the square root of the mass for several experimentally relevant fermionic dipolar molecules. Note that the effective dipole moment in the trap is somewhat smaller, since the dipoles are not aligned perfectly. It is however of the same order of magnitude as the permanent dipole moment; in [5] the JILA group reports a maximum value for the average dipole moment in the laboratory frame of about 40% of the permanent moment.
From the figure, we see that the molecules 7 Li 40 K, 23 Na 40 K, and 40 K 137 Cs and moreover 6 Li 87 Rb and 6 Li 133 Cs, which lie outside the figure, have such large values of √ mp, that one will observe stripe formation already in the regime of relatively low density where multilayer effects are important. This demonstrates the experimental relevance of the results discussed in this paper. The 40 K 87 Rb and 6 Li 23 Na molecules on the other hand have small values of √ mp, and the density required to observe stripe formation is so high that interlayer correlations are not important. The molecules containing Li are all chemically unstable [51] in the sense that the reaction 2YX → Y 2 + X 2 is exothermic for Y=Li and X any other alkali metal. Molecules of 23 Na 40 K and 40 K 137 Cs are however chemically stable in this sense, so they are prime candidates for studying the density wave instability.

Experimental issues and detection
Experiments with polar molecules operate at finite temperatures. The JILA experiments with KRb have reported temperatures down to T = 220 nK or T /T F = 1.4 [52]. This is close to but not quite in the degenerate regime, so a decrease in temperature is most likely needed to see many of the predicted phases. Furthermore, in the low-dimensional setups created by optical lattices, heating can occur as the lattice is turned on, demanding that the initial temperature be even lower to reach critical temperatures.
In the strict 2D limit, the molecules are completely confined in the transverse direction. This means that their motion is reduced to a strict planar geometry. In addition, we do not allow any tunneling between the layers. The layer index can be considered an effective spin label. In this case, the non-zero temperature phase transition to ordered states such as the density wave is governed by the Berezinskii-Kosterlitz-Thouless (BKT) transition [53,54], and no true long-range order occurs. One consequence is that the BKT transition temperature is below the transition temperature obtained from mean-field theory.
However, studies of one-dimensional arrays of tubes with two-component fermions that can undergo a pairing transition indicate that weak tunneling between the tubes can stabilize long-range order [55]. We imagine that this could also work for our multilayer system, i.e. by allowing weak tunneling between the layers the densitywave ordering becomes long-range and the transition temperature may be increased from the BKT value to the mean-field prediction. This intriguing possibility will be explored in future work.
For detection of the density wave state in the multilayered setup a number of different techniques are possible. A quantum non-demolition measurement [38,56,57,58,59,60,61] could be used to detect the large density fluctuations close to the phase transition. Alternatively, light scattering experiments proposed to detect dimerized pairing in multilayers [62] could be adapted to the density wave case.

Conclusions
Using a conserving Hartree-Fock approximation, we examined the density wave instability of aligned fermionic dipoles moving in equidistant planes. We found that while exchange correlations suppress the instability, it can be significantly enhanced by correlations between the layers. The inter-layer correlations exhibit an interesting dependence on the dipolar angle, and they result in the density waves in the different layers to be in-phase for all angles. We furthermore demonstrated that the physics of the interlayer correlations can be understood from a classical model. The density wave instability was shown to be experimentally accesible with realistic densities for experiments using 7 Li 40 K, 6 Li 87 Rb, 6 Li 133 Cs, 23 Na 40 K and 40 K 137 Cs molecules. For these molecules, interlayer correlations were furthermore predicted to decrease the critical coupling strength for density waves significantly.
The inversion symmetry of both inter-and intralayer potentials V (k) = V (−k) gives f (−k) = f (k) and ε(−k − q) = ε(k + q) which allows for a shift of the displaced Fermi sea to give for the static (ω = 0) scattering matrix q, k , q). The scattering matrix is then given by Γ(k, k , q) = 1 2 [Γ + (k, k , q) + Γ − (k, k , q)] and thus has poles if and only if at least one of the symmetrized Γ's has a pole (unless the poles of Γ + and Γ − happen to cancel).
The Fermi sea is deformed by the dipole-dipole interaction, so we write the sum as an integral in polar coordinates and for fixed ϕ parametrize the norm integral by x = k /k F (ϕ ). Letk be unit vector in the direction of k, then the integral becomes An approach to solving the integral equation (A.1) for Γ ± ij is to choose suitable abscissa k α and weights w α for points in the single Fermi sea and approximate the integral by a sum. Then (11) becomes a matrix equation where (i, j)'th block of the matrices describes the interaction on layer i from layer j. We suppress the common q-dependency and introduce the weights as a diagonal matrix W . In the double indices of layer number (roman letter) and k-grid point (Greek letter) the equation reads where the diagonal matrix χ f has entries χ f l,γ = kγ ε lkγ −ε lkγ +q . The irreducible interaction is finite, so the scattering matrix diverges when the matrix in the brackets becomes singular. This happens when the matrix I iα,lγ W lγ χ f lγ has an eigenvalue of 1 (see also Appendix B). For I − the direct interaction cancels and the layers decouple as there is no exchange between different layers. For the single layer we find that the matrix has only negative eigenvalues, so Γ − never contributes to the divergence.
For fixed ϕ c = π/2 and q = 2k F (ϕ c ) we vary g until the first eigenvalue crosses 1. The difference of the single particle energies of the particle and hole in the denominator of χ f gives rise to an integrable singularity at the edge of the integration region at ϕ = 3π/2, k = 2k F (ϕ c ), so we partition the ϕ interval and cast 60 points in the [−π/2+∆ϕ, 3π/2−∆ϕ] and 10 points in [3π/2−∆ϕ, 3π/2+∆ϕ]. For the x = k /k F (ϕ) variable we cast 10 points in [0, 1 − ∆x] and 30 points in [1 − ∆x, 1]. We choose the abscissa and weights in the four intervals by a Gauss-Legendre quadrature rule. For ∆ϕ = 0.05 and ∆x = 0.02, a tripling of the number of points in either interval gives a change of the critical coupling strength that is less than the absolute tolerance 2·10 −3 .
The matrix depends on g explicitly from the potential in I, but also indirectly through the self energies and Fermi function. The iterative procedure is as follows: For the current guess for g c , calculate the deformation f (ϕ) of the Fermi sea according to Appendix A.1, rescale the k-points and calculate the matrix I ± iα,lγ g W lγ χ f lγ , find the largest eigenvalue λ 1 and set the new guess to g = λ −1 . The iteration is terminated when the absolute change in g is less than 2 · 10 −3 .
Appendix B. Scattering matrix near the critical point As g → g c for fixed q = q c the matrix in the square brackets in (A.2) becomes singular and the response blows up. If we make a spectral decomposition of I ± iα,lγ W lγ χ f lγ , the eigenvector with eigenvalue 1 will dominate in the vicinity of the divergence. To see this (for simplicity write I = I ± iα,lγ and for the diagonal matrix W lγ χ f lγ = D), consider a g < g c , so that ID has eigenvalues λ i < 1∀i. We need to find the inverse (1 − ID) −1 , so we make an eigenvalue decomposition for ID, say ID = V AV −1 with A = diag(λ 1 , . . .). Then (ID) n = V A n V −1 → 0 for n → ∞. Now take the geometric series As g → g c , the first 1/(1 − λ 1 ) is divergent.
For g < g c we have the inversion in (B.2) above and now wish to examine the diagonalizing matrix V . The matrix I + D has the symmetry found by switching blocks by multiplying with C = 0 1 1 0 : C has C −1 = C and commutes with both I + and D since all 3 are (block) symmetric, thus it commutes with the product. C is it's own inverse, so it has eigenvalues ±1. Each has n linearly independent eigenvectors of the form v = [w T , ±w T ] for λ = ±1 respectively. Because the matrices commute, the eigenvectors of I + D can be chosen to have the same form and the numerics show that the eigenvector with the largest eigenvalue has the +-form. The diagonalizing matrix can be thus chosen to be of the block form where the n × n matrices W i are determined by the eigenvectors. As g → g c only the first 1/(1−λ 1 ) is divergent in the decomposition (B.2). To find the response functions, we need to sum all matrix elements in each block (7). Near the divergence all but the first entry in the diagonal matrix can be ignored, so where v 1 = 1 2 w 1 w 1 , and w 1 is the first column of W 1 from (B.4). From (B.4) we see, that the first row of V −1 has the structure y T We see that the blocks of the scattering matrix are all the same close to the divergence Γ ij = Γ c . The response function is given by (7), so χ i,j is found by multiplying the block matrix in (B.6) by the diagonal matrix χ f on both sides and then summing all matrix elements within each n × n block and finally adding χ 0 , the response function for non-interacting dipoles with HF single particle energy. The g-dependency of the I + D-matrix comes both directly from the dipole-dipole interaction in I + = gĨ + wherẽ I + is determined by the vectors k, k , q and from the self-energy in χ f . The critical g c is defined by the nonlinear eigenvalue equation g(Ĩ + )D(g)v = λ 1 v, so by expanding near g = g c we see 1 − λ ∝ g c − g, so after the summation, the 2 × 2 structure is χ 1,1 χ 1,2 χ 2,1 χ 2,2 = χ 0 0 0 χ 0 + The numerics shows that the second largest eigenvalue value has an eigenvector of the type v T 2 = [w T 2 , −w T 2 ]. Near the corresponding g, the inversion formula for {1 − (I + D)} does not work since (I + D) n does not go to zero. This is because the system has already gone unstable. If we ignore this fact, say if the highest eigenvalue actually was the one with this eigenvector, the rank one approximation to the inverse would be the outer product of the (n + 1)th column and row vectors of V and V −1 from (B.4), y T 2 = 1 2 [z T 2 , −z T 2 ] so the response near this g c2 would be χ 1,1 χ 1,2 χ 2,1 χ 2,2 = χ 0 0 0 χ 0 + (B.8) .

Appendix C. Classical calculations
First we show that a dipole in one layer does not feel a homogeneous distribution of dipoles in another (infinite) layer. The dipoles are in the xz-plane and aligned by the external field as p = D(sin θ E , 0, cos θ E ) and the two layers are separated by a distance d. The angle between the relative vector and the dipole orientation is cos 2 θ r = (p · r) 2 r 2 = (x sin θ E + d cos θ E ) 2 r 2 , (C.1) while dipole-dipole interaction is given by (2). The calculations is done in a polar coordinate system in the xy-plane, ρ 2 = x 2 + y 2 , so V plane = D 2 n 2D This shows that any homogeneous background density cancels, so the interaction with a secondary layer is only dependent on the deviation from the average density. A density lower than the background can be modelled by changing the sign of the potential.
Since the homogeneous background density does not contribute to the interlayer interaction, we can model the density wave phase by a periodic series of lines with alternating sign on the dipole-dipole interaction. The distance (wavelength of the density wave) between two lines with the same change in density is denoted λ c . The full quantum mechanical calculation (13) shows that the collective eigenmodes correspond to in-phase and anti-phase density modulations in the two layers, so we only consider these two extremes.
We calculate the interaction between a single dipole in a stripe in one layer and the stripes in the other layer by splitting the lines into two sub-series. The first has lines directly on top while the lines of the other series are shifted by λ c /2. For the in-phase density modulations the first sub-series has excess density while the second has decreased density. For the anti-phase configuration the signs of the interactions are switched, but the geometry is otherwise the same. Using the coordinate system from figure 1, the stripes are parallel to the x-axis. The relative coordinates between the dipole and a point in line n in the first sub-series and second sub-series is respectively. For convenience we define a 2 n = n 2 λ 2 c + d 2 or a 2 n = (n + 1 2 ) 2 λ 2 c + d 2 as the squared distance in the yz-plane in the two cases. The dipoles are assumed to be smeared out along the stripes with a linear dipole density of γ. So the interaction with line n is 1 − 3 cos(θ r ) 2 r 3 = 2D 2 γ cos 2 θ E a 2 n − 2d 2 a 4 n (C.4) Summing the contributions from all stripes in both sub-series gives the interaction between a single dipole in a stripe with all the stripes in the other layer V classical = ∓ 2π 2 D 2 γ cos 2 θ E λ 2 c csch 2 (πd/λ c ) + sech 2 (πd/λ c ) , where the − is for the in-phase and + is for the anti-phase configuration.