{'fraction_non_alphanumeric': 0.03361539353193462, 'fraction_numerical': 0.01336888064603377, 'mean_word_length': 4.303688524590164, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Density Waves in Layered Systems with Fermionic Polar Molecules

A layered system of two-dimensional planes containing fermionic polar molecules can potentially realize a number of exotic quantum many-body states. Among the predictions, are density-wave instabilities driven by the anisotropic part of the dipole-dipole interaction in a single layer. However, in typical multilayer setups it is reasonable to expect that the onset and properties of a density-wave are modified by adjacent layers. Here we show that this is indeed the case. For multiple layers the critical strength for the density-wave instability decreases with the number of layers. The effect depends on density and is more pronounced in the low density regime. The lowest solution of the instability corresponds to the density waves in the different layers being in-phase, whereas higher solutions have one or several adjancet layers that are out of phase. The parameter regime needed to explore this instability is within reach of current experiments.


Introduction
After the great successes of cold atomic gas physics using neutral atoms with short-range interaction [1,2], many groups have now set their goals on obtaining ultracold samples of polar molecules that have an anisotropic longrange interaction [3,4,5,6,7,8,9,10]. These can, however, lead to strong losses and the design of experimental geometries that reduce these effects are now becoming a reality. In particular, the use of two-dimensional geometries can reduce losses and at the same time very interesting many-body phases in both single-and multilayer configurations have been proposed [11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]. One such proposal concerns the potential instability of a single two-dimensional layer with polar fermions toward the formation of density-waves as the polarization of the molecules with respect to the layer plane is varied [28,29]. However, the systems of current experimental interest are not single-layer [10], and the effect of adjacent layers is therefore of concern.
Using linear response within the the random-phase approximation, we consider how interlayer interactions influence the density-waves instability and how the critical strength is modified by interlayer terms. In order to estimate the effects of exchange terms, we use many-body local field factors. This approach has been successfully applied to electron systems. We find that the instability is enhanced by the presence of in-phase density-waves in neighboring layers. The effect depends on the density of fermions in each layer and is most pronounced in the low density limit where the critical value is inversely proportional to the number of layers. The latter effect is largely insensitive to the inclusion of exchange terms, Fermi surface deformation, or changes in the effective mass. The density-wave instability will therefore occupy a larger region of the zero-temperature phase diagram for a multilayered system as compared to a single layer system.

Linear Response and Effective Interaction
We consider a multilayer system of fermions with dipole moment D and mass m confined in planes parallel to the xy-plane and separated by the distance d. In the direction normal to the planes, all dipoles reside in the lowest quantum level which we take to be a Gaussian of width w, i.e. φ(z) ∝ exp(−z 2 /2w 2 ). The dipole moments D are aligned by an external field forming an angle θ with respect to the normal of the planes and with a projection onto the planes which is parallel to the x-axis. The experimental setup is illustrated in Fig. 1. Two dipoles separated by r interact with the potential V (r) = D 2 (1 − 3 cos 2 θ rd )/r 3 where θ rd is the angle between D and r. We assume that the layers all have the same density n of fermions.
To obtain the instabilities of the multilayered system we use linear response theory and the random-phase approximation (RPA) as was done for the case of a single layer in [28,29]. Within the RPA framework, the densitywave instability occurs at the poles of the density-density response function. To treat several layers we extend the θ (a) The experimental setup and the density waves corresponding to the lowest (a), the next lowest (b), and the highest (c) eigenmodes for the three layer case. The dipole moments form the angle θ with respect to the normal of the planes. Their projection onto the planes is parallel to the wave fronts.
RPA to a multilayer (or multicomponent) system. We can write a general density fluctuation in response to an external potential, φ ex , in momentum (q) and frequency (ω) space as where δρ is a vector quantity containing the disturbances in each layer as entries. Likewise, χ(q, ω) is now in general a matrix of response function with entries χ ij (q, ω).
The interactions between the various layers produce an induced potential which we write as where the matrix V (q) ij = V ij (q) contains the interaction between layers i and j. We approximate the system response to that of a non-interacting Fermi gas responding to both the external and the induced disturbance. We then have where χ 0 (q, ω) ij = δ ij χ 0 i (q, ω) is the matrix of response functions of the non-interacting system which is of course diagonal. Combining Eqs. (1) and (3) we arrive at the following matrix equation the response function In the case of a single layer this equation reduces to the standard RPA expression for the density-density response function. Here we are interesed in density-wave instabilities in the static limit ω = 0 and we have to determine the singularities of χ(q). By inversion, we see that these occur when and this is the equation that we will solve below. We assume here that the density in each layer is the same, so that the non-interacting response functions are all the same, i.e. χ 0 i (q) = χ 0 (q), and are given by where k = 2 k 2 /2m and f is the Fermi distribution.
In the two-dimensional case of interest here we have the explicit expression [30] where k F is the Fermi momentum and θ(x) is the Heaviside step-function. For simplicity, we ignore any Fermi surface deformation due to the dipolar interaction [28]. We will briefly comment on the influence of such effects in Sec. 3.5.

Exchange Corrections
The RPA analysis above neglects the role of exchange interactions. In the single-component Fermi system we consider here, the exchange effect can be significant. As an example, we note that for a momentum-independent potential, the exchange correction would completely cancel the direct term in a Hartree-Fock calculation. The effective dipolar interaction that we discuss in the next section depends, however, linearly on momentum. The effects of exchange can be included via the Hartree-Fock RPA approximation. This unfortunately involves a nonlocal interaction making the resulting numerics somewhat involved. We will not pursue such calculations here, but rather follow the simpler local field factor approach that has been very successful for the electron liquid [31]. It attempts to include the intrinsically non-local effects of the exchange term through the introduction of an effective local 'exchange' potential, in similar spirit to the highly successful density-functional method. This approach has been applied to two-dimensional double-layer electron systems (see for example Ref. [32]) which is a system closely analogous to the one studied here.
In the multilayer setup considered here, we must be careful when including exchange corrections in the correct places. Since we assume that there is no tunneling between the layers, the layer index is effectively a spin coordinate, and we therefore have no exchange corrections for the interlayer interaction. This means that only the diagonal terms in Eq. (4) have to be modified. We write a diagonal entry in the form where G(q) is the momentum-dependent local field factor. There are various more or less sophisticated ways to calculate this factor through self-consistent numerical methods [31]. As we are only interested in estimating the effects of exchange correlations on the density wave instability, we will follow a more intuitive approach originally introduced by Hubbard [33]. The Pauli principle introduces the so-called 'exchangehole' in Fermi systems. For large q, i.e. short length scales, the exchange-hole cancels the direct interaction and G(q) → 1 for q → ∞. For q → 0, i.e. for long distance, the exchange effect should not play a role and in turn G(q) → 0. Between these limits, the detailed functional form of G(q) of course depends on the particular form of the bare potential V 0 (q). As we are only interested in the qualitative effects of the exchange correlations, it is sufficient to use the simple function G(q) = 2 π tan −1 (q/s) which interpolates between the q = 0 and q → ∞ limits above. Here s is the natural scale in the problem at hand; we take s = 2k F .

Effective Dipolar Interaction
The direct dipole-dipole interaction has an intra-and an interlayer part in our multilayered setup. The Fourier transform of the former can be written [34] where q = |q| and α is the azimuthal angle between the wave vector q = (q x , q y ) and the projection of D onto the plane which is parallel to the x-axis. P 2 (x) is the second Legendre polynomial, and we have defined the function 2)] exp(x 2 /2) with erf(x) the error function. To obtain this formula, the z-direction confining the dipoles in the layers have been integrated out.
The interesting angular dependence of the intralayer interaction is contained in the function ξ(θ, α) = cos 2 θ − sin 2 θ cos 2 α. This function provides the anisotropy in momentum space which is absent at θ = 0 when the dipoles are oriented perpendicularly to the layer. For w d, the interlayer interaction can be written as [17] This approximation holds very well for small w and deviates less than 10% for w = 0.2d. As argued in [28,29], the most unstable direction is found at α = π/2. This is a configuration where the density-wave is perpendicular to the x-axis in order to reduce the side-by-side repulsion of the dipoles while optimizing the attraction from the head-to-tail setup, see Fig. 1. In this case we have ξ = cos 2 θ. The first term in Eq. (9) which is constant in momentum space can be discarded since we are working with a single-component Fermi system [28,29]. As discussed in [28], the critical value in a single-layer has some dependence on θ. Here we are interesting in the effects of multiple layers and we thus fix θ at cos 2 θ = 1/3, but our results can be easily mapped to a different angle through the substitution With the choices above, the intralayer interaction becomes whereas the interlayer interaction is simply multiplied by a factor of 1 3 . For wq 1, Eq. (11) reproduces the potential used for the single layer case in [28,29]. As we are mostly concerned with the effects of multiple layers, we will also assume wq 1 for simplicity, i.e. we assume that V 0 (q) is linear in q. This linear momentum dependence of the intralayer potential was used in [29] to argue that the density-wave instability must occur at some dipole strength always. We note that the most unstable mode is expected to be at q = 2k F (neglecting the effects of Fermi surface deformation). For consistency, we must therefore have that 2k F w 1. In terms of the density of a single layer, n, this condition reads w √ 16πn 1. Thus, either the density must be small or the transverse confinement strong. In terms of typical physical scales in experiments, we have Using experimentally relevant values [10] d = 0.5µm and assuming w/d = 0.1, we find that 2k F w < 1 for densities n < 8 · 10 8 cm −2 which is fulfilled by current experiments.

Instability Conditions
We first recapitulate the findings for a single layer in the RPA neglecting exchange. The instability equation is Assuming that the instability occurs first at q = 2k F , we have χ 0 = − m 2π 2 and V 0 (2k F ) = −4πD 2 k F /3. We thus have the relation where D c is the critical dipole strength [29]. We define a dimensionless measure of the strength and density; g := 2mD 2 k F /3 2 . We thus have the critical value g 0 = 1.
To include exchange, we need to make the substitution To highlight the effects of the multilayer setup, we now proceed to discuss the bi-and trilayer cases without the 1−G(q) factors and defer the discussion of exchange corrections to Sec. 3.4.

The Bilayer
For the case of two adjacent layers we get the following algebraic equation from Eq. (5) We note immediately that if we set V 1 (q) = 0 we recover the usual RPA condition for density instabilities. It is also clear at this point that the bilayer will have a smaller D c than the single layer above since the V 1 (q) term is negative. If one considers the interlayer interaction from the point of view of induced interactions this is no surprise as such interaction are usually attractive at lowest order. At q = 2k F we can solve the equation above and find the lowest critical value for a density-wave instability in a bilayer For k F d = 1, we get roughly a 12 percent reduction, whereas for a lower density of k F d = 0.5 the difference is 27 percent. The other solution to the bilayer equation is Solving for the corresponding zero eigenmodes of χ(q, 0) −1 from Eq. (4), we find Here, δρ i = δρ i (2k F ) is the density fluctuation in layer i. We see that the density waves in the two layers are inphase for the lower solution g b and out of phase for theg b solution. Thus, the instability is enhanced by the density waves in neighboring layers being in-phase gaining more attractive head-to-tail energy and minimizing the side-byside repulsion. Likewise, when the density-waves are outof-phase the instability is suppressed.

The Trilayer
The case of three layers produces the algebraic equation where we have introduced the notation V 2 (q) for the interlayer potential of the two outer layers that are a distance 2d apart. Note that V 2 (q) differs from V 1 (q) by a factor of exp(−qd) and we thus expect it to be a much smaller quantity than V 1 (q) at 2k F . In light of the above, we therefore first consider the simpler case of V 2 (q) = 0, i.e. we include only nearestneighbor interactions. This means that the equation for the instability factorizes and the condition becomes that either 1 − χ 0 (q)V 0 (q) = 0 or Clearly the latter condition produces a lower critical value and we find for the trilayer with V 2 (q) = 0 that where the asterisk indicates that we include nearest-neighbor interactions only. As expected the trilayer has a reduced critical value. A naive guess for the trilayer might be to multiply the interlayer strength by a factor of 2 and then consider it as a bilayer problem. However, our result demonstrates that the enhancement is only by a factor of √ 2. As it turns out the trilayer equation has a rather simple analytic solution. The roots are where a = exp(−2k F d) and we have listed them in order of increasing magnitude. The top solution is always less than one and decreases with k F d (we denote it by g t in the following) whereas the others are always larger than one and increase with k F d. The corresponding eigenmodes are where layer 2 is the one in the middle, see Fig. 1. Again, the lowest solution corresponds to the density-waves in the different layers being in-phase with amplitude now being the largest for the layer in the middle. The second eigenmode has the middle layer out of phase and of larger magnitude than the outer layers. This is the same situation as theg b solution for the bilayer, only now the out of phase effect is more costly. This is also reflected in the fact that this solution is always larger thang b for any k F d, whereas the opposite holds for the lowest solution, i.e. g t < g b . The last solution is an eigenmode with density waves of the outer layers out of phase and no amplitude change in the middle layer (within the RPA). We sketch in Fig. 1 the density waves for the three eigenmodes. The physically relevant mode for the instability is of course the one with the lowest critical value corresponding to the density waves in the planes being in-phase. In Fig. 2 we plot the (lowest) bilayer and trilayer critical values at which the density-wave instability appears for k F d ≤ 2. The single-layer critical value, g 0 = 1, is approached asymptotically, however, for the range plotted the multilayer cases are all below that value by at least a few percent. As expected the bilayer is always above the trilayer value. For k F d → 0, the critical value undergoes the largest reduction which is a factor of two for the bilayer, while for the trilayer it is a factor of 3 (we return to this fact below). We also compare the trilayer with and without the interaction of the two outer layers, V 2 (q). When excluding the term, we see a larger critical value, g * t , for all k F d than when taking it into account in g t . The additional attraction of V 2 (q) thus reduces the critical value as one would expect.

Multiple Layers
For more than three layers we expect similar behavior as seen above, i.e. a critical strength that decreases with decreasing k F d. In the limit k F d → 0 we can in fact find the exact solution for g N for any number of layers, N . Here the matrix in Eq. (5) simplifies considerable since it has 1 − g in all diagonal and −g in all the non-diagonal entries. It is easy to verify that a vector consisting of ones in every entry is an eigenvector of this matrix with eigenvalue 1−N g. We thus conclude that the system is unstable towards the formation of in-phase density waves in all planes for the critical coupling strength g N = 1/N . The limit of very small k F d should therefore approach this simple value. This limit is clearly seen for the lowest critical value in the bi-and trilayer cases above with the corresponding eigenmodes approaching one in all entries. If we take this limit by reducing d while keeping k F constant, we see that for large N the critical value approaches zero as the layers come closer. Here the (non-interacting) system is equivalent to that of N spins moving in two dimensions. The response function Eq. (6) is then multiplied by a factor of N which reduces the critical value by a factor 1/N .

Exchange effects
The exchange correction has to be included in the diagonal terms of the response function only as discussed in the previous section. This means that the eigenvectors corresponding to the critical couplings are the same irrespective of whether the exchange effect is included or not. In the large k F d limit, the off-diagonal terms of the interaction are negligible, and we thus obtain the critical value g = 1/(1−G(2k F )). Without exchange the single-layer result is recovered, i.e. g = 1. Using the value G(2k F ) = 1/2 as estimated in the previous section, the limit is a factor of two larger. These arguments make it clear that the effects of exchange are more pronounced in limit of large k F d where the intralayer correlations dominate. However, the interlayer correlations dominate for small k F d which means that the exchange effects are insignificant in this limit. This means that one of our main results, the 1/N scaling of the critical coupling strength for k F d 1, still holds when exchange is included.
In Fig. 3 we show numerical solutions for the lowest critical values for N = 2, 3, 10, 20, and 30, when neglecting (lower full (blue) lines) and including exchange (upper dashed (black) lines). The expected decrease of g N with N is clearly seen both with and without exchange corrections. For example, at k F d = 1 the critical value is 0.88 for N = 2, whereas for N = 30 it is 0.77 when neglecting exchange. The numbers are 1.57 for N = 2 and 1.23 for N = 30 when including exchange. This trend continues for higher N . Note that all the limits discussed above are clearly confirmed by the numerics. In particular the results with and without including the local field factor to account for exchange approach each other as k F d becomes small and when N grows.
The corresponding eigenmodes all have the density waves in the layers in-phase. There are also other solutions with larger critical coupling strengths as for the bi-and trilayer cases. The eigenmodes can be analyzed in similar fashion and one finds that the solutions can be organized according to the number of adjacent layers that are out of phase with each other with the lowest solution (plotted in Fig. 3) fully in-phase across all layers. We speculate that these higher modes correspond to collective modes in the striped phase. This will be examined in the future.

Effective mass and Fermi surface deformation
Finally, we briefly address the question of influence of effective mass and Fermi surface deformation caused by the dipolar interaction. For the single-layer case, these corrections have been calculated in Ref. [28]; it was found that these terms pushes the critical value up by about 25% for cos 2 θ = 1/3 and about 18% for θ = 0. In our setup this factor has to be included as a prefactor of χ 0 in Eq. (4), i.e. it effectively amounts a redefinition of our g. The neglect of effective mass and deformation effects means that our results represent lower bounds. Note again that the eigenvectors for the unstable modes are unaltered by these corrections.

Phase Diagram
In the multilayer setup, the interaction parameter, U = mD 2 / 2 d, is a convenient dimensionless measure for the strength of interactions in the system. In Fig. (4) we show the zero temperature phase diagram in the (U, k F d) plane for N = 2, 10, and 30. The more layers, the earlier one expects to enter the density-wave regime as before. We also see that one can probe the phase diagram by changing either the dipole moment or the density of fermions. Changing d is also an option. This is, however, somewhat harder as U ∝ 1/d and the lines of constant U k F d = 3g/2 are very similar to the lines shown in Fig. (4). We note that the inclusion of the exchange term causes an interesting plateau of the critical values for large N at k F d ∼ 1. This implies that there can be a large region with U 2 and k F d 1.5 where the system is not unstable towards the formation of density waves. This is valuable for the study of other phases like superfluidity which persists to small U [18].  The regime of validity of the RPA approach augmented by the local field factor when applied to dipolar systems can be related to the corresponding situation for the electron liquid. In the latter case the RPA is known to provide reasonable results in the high density limit while it performs poorly at low densities where the Coulomb to kinetic energy ratio, r s , becomes large [31]. However, for dipolar systems the interaction dominates in the high density limit whereas the low density limit is weakly interacting. We thus expect the RPA to be accurate for low densities and weak dipolar strengths, i.e. when g 1. This is precisely the case for the large N limit which is our main interest in this work.

Competing Phases and Finite Temperature
The zero-temperature phase diagram for density-wave instabilities presented above needs to be considered in the light of other possible ground-states of the layered dipolar system. In the case of a single layer and in the weak couping limit, a p-wave superfluid state was proposed [18]. Likewise, a region of negative compressibility leading to collapse of the system was found [18,28], although this happens outside the parameter regime considered here. For several layers, the system can become superfluid with the Coopers pairs formed between dipoles residing in different layers [24,25,26,27]. In the strong-coupling limit, a single layer of dipoles can also form a Wigner crystal with a symmetry which depends on how the dipoles are aligned with respect to the plane [14,15,16,35]. The presence of bound states in single and multilayer configuration of both fermionic and bosonic dipoles has also been given a great deal of attention recently [20,36,37,38,39,40,41]. For strong coupling, chains of multiple dipoles in bound states could be the relevant degree of freedom in the system and the effective interaction of such constituents should determine the ground-state, and could be very different for odd fermionic chains as opposed to even bosonic ones. How the phase diagram of single-and multilayer system at zero temperature maps out is an extremely interesting topic for future research.
At finite temperature one expects the physics to be governed by the Berezinskii-Kosterlitz-Thouless transition (BKT) [42]. In the bilayer case the BKT physics is contained in the pairing order parameter in the weak-coupling limit or in a condensate of bosonic dimers in the strongcoupling limit [26]. For multiple layers similar dimerized phases are expected that are governed by the BKT transition [24]. The universal relation for the critical temperature scales with the superfluid density As the latter is proportional to the total density for strong coupling, the low density regimes can be difficult to access. We speculate that the interlayer interactions could help stabilize the low-temperature phases of the system and in turn easier to access experimentally as compared to a single layer. Again this is a topic for future research.

Conclusions
We have considered the density-wave instability of dipolar fermionic polar molecules confined to a stack of twodimensional layers. As the number of layers increases we find a reduction of the critical strength to enter the densitywave regime at all densities. The corresponding density waves are in-phase in all the planes. In the low density limit the critical strength even approaches zero as the number of layers grow.
We thank M. M. Parish for numerous discussion and for providing valuable references.