Density-wave instability and collective modes in a bilayer system of antiparallel dipoles

We consider a bilayer of dipolar particles in which the polarization of dipoles is perpendicular to the planes, in the antiparallel configuration. Using accurate static structure factor S ( q ) data from hypernetted-chain (HNC) and Fermi HNC calculations, respectively for an isolated layer of dipolar bosons and dipolar fermions, we construct effective screened intralayer interactions. Adopting the random-phase approximation for interlayer interactions, we investigate the instability of these homogeneous bilayer systems towards the formation of density waves by studying the poles of the density–density response function. We have also studied the collective modes of these systems and find that the dispersion of their antisymmetric collective mode signals the emergence of the density wave instability as well.


Introduction
Layered structures combine interesting physics of the low-dimensional systems with the additional tunability coming from the interlayer interactions and tunneling. In condensed matter systems, several unique phenomenon such as Coulomb-drag effect [1], formation of indirect excitons and their eventual condensation [2], density wave instabilities (DWI) and Wigner crystallization [3] and fractional quantum Hall effect [4,5] in layered structures have been the subject of an immense interest in the past few decades. The isolation of graphene and other layered van der Walls materials [6] in recent years have enormously raised this enthusiasm.
Ultracold atomic and molecular systems, on the other hand, with their impressive controllability have become natural simulators of the condensed matter and many-body theories. In particular experimental progress in trapping and cooling atoms with large magnetic moments and polar molecules, opened up a new and interesting area of exploring quantum many-body systems with large and anisotropic dipole-dipole interactions [7][8][9]. Both polar molecules [10][11][12][13][14] and atoms with large permanent magnetic moments [15][16][17][18] have been trapped and cooled down. Very recently, the droplet crystal phase of atomic dysprosium Bose-Einstein condensate (BEC) have been directly observed by Kadau et al [19].
In bulk geometries, the attractive part of the dipole-dipole interaction could in principle lead to instabilities, as in Bose-Einstein condensate collapse [20] or chemical reactions between particles [9]. Therefore, it is usually favorable to confine the dipolar gases into quasi-two or one-dimensional geometries, and use an external electric or magnetic field (depending on the nature of dipoles) to polarize all dipoles in the same direction. As might be expected, layered structures are also a configuration of great interest which one can tune the attractive interactions and pairing between different layers without the fear of having chemical reactions.
While the stripe or density-wave phase is naturally expected in an isolated two-dimensional (2D) system of tilted dipolar bosons [21,22] and fermions [23][24][25][26][27][28][29] due to the anisotropy of the dipole-dipole interaction, this instability has been the subject of much dispute in the limit of perpendicular dipoles, where the inter-particle interaction is isotropic. While mean-field approximation [23] as well as density-functional theory (DFT) [30] and Singwi-Tosi-Land-Sjölander (STLS) [25] calculations all predict stripe phase formation at relatively low interaction strength for 2D dipolar fermions, quantum Monte Carlo (QMC) simulations suggest that the stripe phase never becomes energetically favorable, up to the liquid-to-solid phase transition for perpendicular bosons [22] and fermions [31].
In double-layer structures, both bosonic and fermionic systems have attracted a lot of attention so far. The ground state properties and instabilities of fermionic bilayers have been studied within the Hartree-Fock [32][33][34] as well as STLS methods [35]. The QMC simulations [36], as well as DFT calculations [37], have been employed to study the crossover from BEC to Bardeen-Cooper-Schrieffer state too. For bosonic bilayers, on the other hand, Hufnagl and Zillich [38] have used the hypernetted-chain (HNC) approximation to calculate the ground-state quantities of a bilayer system of tilted dipolar bosons. Then using the correlated basis function (CBF) method they obtained its dynamical properties. It has been also suggested that a bilayer system of dipolar bosons becomes a self-bound fluid when the polarization of dipoles in two layers is opposite [39]. More recently, the competition between single-dipole and dimer condensation in a bilayer of perpendicular dipolar bosons with parallel polarization, i.e. the same direction of polarization in both layers, have been investigated using QMC method by Macia et al [40]. They have observed that the pair superfluidity dominates over the singleparticle superfluidity at very strong interlayer couplings, i.e. when the separation between two layers is much smaller than the average intralayer distance between two dipoles. The dynamical properties of the dipolar bosonic bilayer in the atomic and pair superfluid regimes have been looked at using QMC and CBF techniques [41]. The correlation effects in a bosonic bilayer have been extensively studied using QMC simulation for the ground state properties as well as the stochastic reconstruction method and method of moments for the dynamical properties by Filinov [42].
Our aim in this work is to study symmetric bilayers with the equal number of identical dipoles in each layer, whose moments are aligned perpendicular to the 2D-planes, over a wide range of the inter-layer and intra-layer couplings. We investigate bosonic and fermionic dipolar systems on equal footing, but consider only antiparallel polarization of dipoles in two layers (see, figure 1 for a schematic illustration). Perpendicular alignment of dipoles makes both the intralayer and interlayer interactions isotropic. While the bare intralayer interaction is purely repulsive, the bare interlayer interaction could be either repulsive or attractive, depending on the in-plane separation of two dipoles. In our antiparallel configuration, the interlayer interaction is repulsive at small inplane separations and becomes weakly attractive at large separations (see, equations (1) and (2), below). We should note that in bilayers with a parallel polarization of dipoles in two layers, the dominant interlayer interaction is attractive. At small layer spacings, this in principle leads to the pairing between dipoles of two adjacent layers. This problem has been extensively studied for both bosonic [40,41] and fermionic systems [36,43]. In this work, instead, our focus is on bilayers with the antiparallel polarization of dipoles. In this configuration the pairing is either absent or extremely weak [44,45] and therefore is not expected to affect the strong interlayer screening at small layer spacings [46]. We investigate the possibility of the instability of a homogeneous fluid towards the formation of inhomogeneous densities, i.e. density waves. For this purpose, we look at the poles of the static density-density response function. The effective intralayer interactions are obtained from an accurate HNC and Fermi HNC (FHNC) results for the static structure factor of an isolated 2D layer of bosons [47] and fermions [48], respectively, combined with the fluctuation-dissipation theorem. We have treated the interlayer interactions within the random-phase approximation (RPA) [49]. A similar study of the instability of a homogeneous liquid with respect to the inhomogeneous phase of charge density wave has been studied in a variety of quantum charged systems ranging from single-layer electron gas [50] to electronelectron and electron-hole double-layers [51][52][53][54], charged Bose systems [55] and superlattices [56,57].
We also find semi-analytic expressions for the full dispersions of in-phase and out-of-phase collective modes (i.e. zero-sound modes) from the poles of the dynamical density-density response function. For both bosonic Figure 1. Schematic illustration of a bilayer system of dipoles with the antiparallel polarization of dipolar moments in two layers. d refers to the layer spacing between two layers. and fermionic bilayers, the signatures of the emergence of DWI show up in the dispersions of these collective modes.
The rest of this paper is organized as follows. In section 2, we introduce the density-density response function of our system and explain how effective intralayer interaction could be obtained from the static structure factor. In section 3, we look at the density wave instability for bilayer systems of dipolar bosons and dipolar fermions. In section 4 we calculate the collective modes of these bilayer structures and investigate their dispersions in the vicinity of the DWI. In section 5, we summarize and conclude our main findings. Finally, in the appendix, we report some analytic results for bilayer of dipolar fermions using the so-called Hubbard local field factor (LFF) for the effective intralayer potential.

Density-density response function and the effective interactions
We consider two identical 2D planes of dipoles, separated by the distance d. No tunneling is allowed between two layers. Therefore, layers are coupled together only through the dipole-dipole interaction. All dipoles are assumed to be polarized perpendicular to the planes, but the relative direction of this polarization is assumed to be antiparallel in two layers (see, figure 1). The bare intralayer and interlayer interactions respectively read [39] where C dd is the dipole-dipole coupling constant, r is the in-plane distance between two dipoles. After Fourier transformation one finds [33,58] Here erfc(x) is the complementary error function and w is the short distance cut-off introduced to heal the divergence of Fourier transform of the intralayer interaction.
In this work we are interested in the density-wave instabilities and collective density modes of this double layer structure. For this we begin with the following matrix equation for the density fluctuations [49] is the external potential applied to layer j and c w ( ) q, ij is the density-density response function, whose matrix form reads

eff
Here is the non-interacting density-density response function, and Eigenvalues of the density-density response matrix c ŵ , , s d are the symmetric and antisymmetric components of the effective potentials.
The non-interacting density-density response function w P( ) q, of a 2D system is analytically well known. In the case of 2D bosons, it reads  w e w e P = + - where n is the particle density in each layer and is the single-particle energy of dipoles of mass m. The full analytic form of for fermions is slightly more complicated and could be found e.g., in [49].
The exact form of the effective potentials are not known, and one has to resort to some approximations. In the celebrated RPA [49], the effective potentials are replaced with their bare values. But as the effects of exchange and correlation become more significant with increasing interaction strength, naturally the RPA which entirely discards these effects needs to be improved at strong couplings. On top of this, as the bare intralayer potential in q-space (3), has an artificial cut-off dependence, a simple application of RPA appears to be not very appropriate for dipolar systems even at weak couplings [34]. In order to overcome both of these problems, we use the fluctuation-dissipation theorem to find an approximate expression for the effective interlayer potential [48]. At zero temperature the fluctuation-dissipation theorem reads [49] Here, ( ) S q is the static structure factor of an isolated 2D dipolar liquid, which can be obtained very accurately both for bosons and fermions e.g., from QMC simulations [31,[59][60][61] or HNC [47] and FHNC [48] calculations. Therefore, the idea here would be to invert equation (9), and find an approximate expression for the static effective interaction in terms of the static structure factor ( ) S q . This is in principle possible if one ignores the dynamical effects, i.e. replaces w ( ) W q, s with a static and real function ( ) W q s . As the effects of exchange and correlation are already included in the static structure factors, these effects will be transferred into the effective intralayer potentials, at least at the static level.
For bosons the integral over ω in equation (9) could be performed analytically, to result in Whereas in the fermionic case, the complicated form of the exact prevents an analytic solution to equation (9), however resorting to the so called 'mean-spherical approximation' (MSA) for the density-density response function is the non-interacting static structure factor of a spin polarized 2D system of fermions [49], again an analytic solution of the frequency integral in the fluctuation-dissipation relation gives Note that, the MSA expressions for the non-interacting density-density response function (11), and the effective interaction of the fermionic system (12), reduce to the corresponding ones of the bosonic system with  ( ) S q 1 0 , which is indeed the correct static structure factor for non-interacting bosons. As already mentioned, the effects of exchange and correlation, entirely ignored in the RPA, are partly included in equations (10) and(12) through the interacting static structure factor ( ) S q . In figure 2 the HNC and FHNC results for the static structure factor of a single layer of dipolar bosons and fermions together with the effective intralayer interactions obtained from equations (10) and(12) are illustrated for several interaction strengths l = k r 0 0 . Here, is the characteristic length scale for dipoles, and p = k n 4 0 . Note that k 0 is indeed the Fermi wave vector k F of each layer for fermionic bilayers, but it is simply a measure of the density for bosonic bilayers.
In the following, we set the interlayer part of the effective interaction to the bare interlayer interaction , as an accurate knowledge of the interlayer static structure factors over a wide range of parameters for both bosons and fermions, is not yet available in the literature. Such an approximation is equivalent to RPA and we surmise it will be adequate for large enough layer separations. We should note that all the properties of these bilayer systems are governed by two dimensionless parameters, namely the intralayer coupling constant λ, and the dimensionless layer spacing k d 0 .

Density-wave instabilities
Density-wave instabilities could be obtained from the poles of the density-density response function given in equation (7) in the static limit, or equivalently from the solution of In fact, for a given system parameters such as the particle density n and layer spacing d, if equation (13) satisfies a solution with a specific wave vector q c , then the homogenous fluid becomes unstable towards the spontaneous formation of density modulations with the wavelength l p = q 2 c c . In the following, we investigate such an instability first for a bilayer system of dipolar bosons and then for a bilayer system of dipolar fermions.

Bosonic bilayers
In the static limit, the non-interacting density-density response of equation (8) reduces to which together with equation (13), gives Now, using the bare interlayer potential(4) and the effective intralayer potential of (10) in equation (15) we find As the static structure factor is positive by definition, the above expression with positive sign will not have any solution which means that no density-wave singularity in the in-phase channel (i.e., c + ) appears. On the other hand, in the out-of-phase channel (i.e., c -) one can find instabilities for suitable values of the interaction strength and layer spacing from the solutions of equation (16) with the negative sign. This means that the maxima and minima of the modulated density in two layers would be shifted by l 2 c with respect to each other. Numerical investigation of equation (16) reveals that (see, figure 3) for  l 1 the density wave instability at a finite wave vector develops below a critical layer spacing d c . At smaller intralayer couplings, the homogenous superfluid phase remains stable up to zero layer separations.
We note that for an isolated single layer, one has = ( ) W q 0 d , and the criteria for the density wave instability becomes which evidently has no solution at any finite q. Therefore, within the approximations we use here, no density wave instability is expected to happen in an isolated 2D system with purely repulsive dipolar interaction. This is in agreement with the QMC findings [22].

Fermionic bilayers
The non-interacting density-density response function of a 2D system of fermions in the static limit is [49] n P = --Q --  The phase diagram in figure 5 illustrates our numerical solution of equation (20). Similar to the bosonic bilayer, instability emerges only in the out-of-phase channel. The main observation here is that at a fixed density, the critical layer spacing for the formation of density waves in fermionic bilayers is slightly larger than the bosonic ones, and no DWI develops at  l 0.5. Figure 6 shows the static density-density response functions of a bilayer system of fermions A similar behavior to the bosonic system is observed also here. The antisymmetric component signals the emergence of DWI as the layer spacing approaches its critical value.

Collective modes
In this section, we turn to the discussion of the collective modes of dipolar bilayers. In symmetric bilayers and in the absence of tunneling between two layers, two collective density modes are simply the in-phase and out-ofphase oscillations of the particle density in two layers. The dispersion of these collective modes could be obtained from the singularities of the density-density response functions c w  ( ) q, at finite frequencies, or equivalently from the zeros of     (10), the dispersion of collective modes read Note that the first term on the right-hand-side of this equation is the Bijl-Feynman excitation spectrum of a single layer dipolar Bose liquid [47]. In the long-wavelength limit, as the static structure factor vanishes linearly is the sound velocity of bosonic system with ¢ = = ( ) ( ) | S Sq q 0 d d q 0 . Unlike the charged boson bilayer [62], both collective modes of a bilayer system of dipolar bosons have acoustic nature. The reason we find same sound velocity for both collective modes relies on the fact that we are using the bare interlayer potential which vanishes linearly at small q and hence does not contribute to the sound velocity (see the second term inside the bracket in equation (24)). One would expect deviations from this simple picture at small d, where the interlayer coupling is strong, but at larger layer spacings both sound velocities should approach the same value. Indeed, this has been verified in [39] for a bilayer of dipolar bosons with antiparallel polarization.
Using the numerical results for the static structure factor from [47] in equation (24), the full dispersions of the collective modes could be readily obtained. The results for w  ( ) q and single-layer collective mode are presented in figure 7 for a fixed value of the coupling constant λ and for different values of the layer separation d. We find that as the layer separation approaches the critical spacing for the density wave instability, the antisymmetric mode w -( ) q touches zero and becomes soft. This occurs at the same q-value that the static density-density response function c -( ) q diverges (see Figure 4). The energy of antisymmetric collective mode at smaller layer separations becomes negative, which is an indication of homogenous liquid phase becoming unstable.

Fermionic bilayers
In finding the collective modes of the fermionic system, one should solve the complex equation is the Fermi velocity, and this solution is valid as long as dispersions lie outside the PHC i.e, F . In the long wavelength limit, using the fact that the fermionic static structure factors are also linear at long wavelength, and therefore the intralayer effective interaction ( ) W q s is finite at q=0, we find is the fermionic sound velocity, and is related to the slopes of both interacting and noninteracting structure factors at the origin through equation (12). As v s F, is always larger than the Fermi velocity v F , the zero sound waves are undamped at the long wavelength for any coupling constant and layer spacing. Interestingly, similar to the dipolar bosonic bilayers, in-phase and out-of-phase collective modes are both linear at long wavelength. Again, the degeneracy of both modes at small q should be valid only at large layer spacings. At smaller separations, the exchange-correlation effects in the effective interlayer interaction will split these two modes. Whether the lower branch will still survive the Landau damping at long wavelengths or not, requires further investigations with a more careful treatment of both intralayer an interlayer correlations.
In figure 8 we show the behavior of collective modes of the fermionic bilayer system w  ( ) q at a fixed coupling parameter λ and for different values of the layer separation. The PHC is also shown in these figures. Collective excitations are well defined only outside this continuum. Note that as the layer separation becomes smaller than the critical value (see, the right panel in figure 8), an unphysical low energy branch at > q k 2 F appears below the PHC.

Summary and conclusions
In summary, we have investigated the instability of a homogenous bilayer system of perpendicular dipolar bosons and dipolar fermions with the antiparallel polarization of two layers towards density waves. Accurate HNC results for the intralayer static structure factor of bosons and FHNC results for the intralayer static structure factor of fermions are used together with the fluctuation-dissipation theorem to extract the static intralayer effective potentials and the RPA is employed for the interlayer interaction. We have observed that for both fermionic and bosonic bilayers, below a threshold intralayer coupling strength λ, no density wave instability emerges. At higher couplings, DWI forms below a critical layer spacing d c . In a given λ the DWI in fermionic bilayers sets in at a larger layer spacing in comparison with the bosonic bilayers. We have predicted that a homogenous bilayer with antiparallel polarization of dipoles in two layers is unstable towards the formation of density waves when the layer separation d becomes comparable or smaller than the average inplane distance between particles k 1 0 , and both of these length scales are much smaller than the dipolar length scale r 0 . We would expect this regime to be readily accessible experimentally with polar molecules whose dipolar length could easily exceed several thousands of nanometers [9].
The full dispersions of the in-phase and out-of-phase zero-sound modes of the bilayer system have been calculated too. We observed that both modes are linear at the long-wavelength limit, independent of the statistics of the particles. Finally, we should note that in the limit of closely separated layers, improvements beyond the RPA in the effective interlayer potential, like the inclusion of exchange-correlation effects might be necessary. Dynamical effects and frequency dependence of the effective potentials would become important in the strongly correlated regime too.

Acknowledgments
This work is supported in part by the Turkish Academy of Sciences TUBA. SHA thanks the hospitality and support of Department of Physics at Bilkent University, during his visits.

Appendix. Hubbard LFF for fermionic dipolar bilayer
Although the semi-analytic expression (12) for the effective interaction of dipolar fermions is neat and carries most information of the exchange-correlation hole through the static structure factor, but it requires an accurate knowledge of the static structure factor in the first place. A systematic approach to go beyond the RPA in electron liquid systems have been through the inclusion of the LFF [49] w where ( ) V q is the bare interaction, and the LFF w ( ) G q, itself is defined thorough equation (A.1) to give the exact dynamical effective potential w ( ) W q, . In practice, one needs to approximate the LFF, and the first and simplest approximation is suggested by Hubbard, which for spin-polarized fermions reads [49] Now, using equations (3), we find where the vanishing cut-off (i.e.,  w 0) is understood on the right-hand side. The interlayer effective interaction is same as the bare one ( ) V q d , as interlayer LFF vanishes within the Hubbard approximation, due to the absence of exchange interaction between particles from different layers [49]. 58. This is pretty close to the λ-threshold obtained from the fluctuation-dissipation theorem in section 3.2 (see figure 5). In figure A1 we have compared the critical layer spacing for the DWI, obtained from the numerical solution of equation (A.4) with the approximate expression(A.6).
The full dispersions of the collective modes within the Hubbard approximation are still given by the general expressions (28)