# Stability and excitations of a bilayer of strongly correlated dipolar Bosons

## Abstract

We study correlation effects and excitations in a dipolar Bose gas bilayer which is modeled by a one-dimensional double well trap that determines the width of an individual layer, the distance between the two layers, and the height of the barrier between them. For the ground state calculations we use the hypernetted–chain Euler Lagrange method and for the calculation of the excitations we use the correlated basis function method. We observe instabilities both for wide, well-separated layers dominated by intra-layer attraction of the dipoles, and for narrow layers that are close to each other dominated by inter-layer attraction. The behavior of the pair distribution function leads to the interpretation that the monomer phase becomes unstable when pairing of two dipoles becomes energetically favorable between or within layers, respectively. In both cases we observe a tendency towards “rotonization”, i.e. the appearance of a soft mode with finite momentum in the excitation spectrum. The dynamic structure function is not simply characterized by a single excitation mode, but has a non-trivial multi-peak structure that is not captured by the Bijl-Feynman approximation. The dipole-dipole interaction between different layers leads to additional damping compared to the damping obtained for uncoupled layers.

###### pacs:

03.75.Hh, 03.75.Kk, 67.85.De, 67.40.Db## I Introduction

Experimental advances in achieving Bose-Einstein condensation (BEC) of atoms with large magnetic moment (Cr Lahaye et al. (2007, 2008), Dy Lu et al. (2011), Er, Aikawa et al. (2012)) and in generating quantum gases of heteronuclear molecules (KRb Ni et al. (2008); Ospelkaus et al. (2008); Ni et al. (2009, 2010), LiCs Deiglmayr et al. (2008), LiK Voigt et al. (2009), RbCs Sage et al. (2005); Takekoshi et al. (2012)) by Feshbach association have lead to a growing interest in effects caused by the dipole-dipole interaction (DDI). The shape of a trapped dipolar condensate Goral et al. (2000); Ronen et al. (2007) and its stability against collapse Lahaye et al. (2008); Koch et al. (2008) have been investigated. The dynamics O’Dell et al. (2003); Santos et al. (2003); Mazzanti et al. (2009); Hufnagl et al. (2011); Macia et al. (2012) has been studied theoretically, but recently also experimentally Bismut et al. (2012). The generation of novel phases with topological order using polar molecules has been proposed Micheli et al. (2006). A recent review of the field can be found in Ref. Baranov et al., 2012. The strength of the DDI can be characterized by the dipolar length , where is the mass of the dipolar atom or molecule and is proportional to the square of the dipole moment. For magnetic moments, is usually much smaller than for electric dipole moments of molecules, which can range to thousands of Å. Since achieving BEC with heteronuclear molecules is much harder than with homonuclear molecules, Er is a promising candiate to reach a much stronger DDI regime with purely magnetic dipole moments fer ().

The two-dimensional limit of a dipolar Bose gas (DBG) polarized perpendicularly to the plane has been studied extensively by quantum Monte Carlo methods Buchler et al. (2007); Astrakharchik et al. (2007); Mazzanti et al. (2009); Filinov et al. (2010); Filinov and Bonitz (2012) for a wide range of dimensionless densities , including high densities where the dipole-dipole repulsion leads to such strong in-plane correlations that the excitation spectrum exhibits a roton similar to the roton in superfluid He, and even higher densities where the ground state of the 2D DBG is a triangular crystal. Such large may soon be in experimental reach because can be very large, as mentioned above. For the 2D DBG with tilted polarization a stripe phase has been predicted recently for sufficiently large tilt angle Macia et al. (2011, 2012). Also the more complicated case of a quasi-2D layer was studied, i.e. of a DBG in a one-dimensional trap . The finite extent in this direction allows pairs of particles to explore the anisotropy of the DDI. Already by using the mean field approach of the Gross-Pitaevskii equation, “rotonization” of a quasi-2D layer of a polarized DBG was found to occur if the strength of the DDI surpasses a critical value with respect to the short-range repulsionSantos et al. (2003); Wilson et al. (2008). The roton in this case is not a signature of the repulsive correlations as in the 2D limit for high densities, but a signature of the attractive correlations for head-to-tail configurations of pairs of dipoles. We have shown this conclusively in Ref. Hufnagl et al., 2011, where even a cross-over between these two kinds of rotons was demonstrated by variation of the trap length. While the effect of attractive correlations on the excitation spectrum can be qualitatively described by mean field methods that optimize only the ground state density by the Ritz’ variational principle, repulsive correlations leading to He-like rotons require the optimization of at least density and pair density. This is exactly what the family of hypernetted chain Euler-Lagrange (HNC-EL) methods does. The HNC-EL method was therefore used in Ref. Hufnagl et al., 2011 in order to investigate both kinds of rotons. A summary of how HNC-EL works is given in the next section, all details about the method can be found in Ref. Krotscheck et al. (1985).

In this work we extend our previous investigation of a single quasi-2D DBG layer Hufnagl et al. (2010, 2011) to a bilayer. The coupling between layers via the long-ranged DDI and the possibility of pairing of dipoles on different layers to form dimers (and more generally of bound dipoles in layers) has been investigated previously Wang (2007); Trefzger et al. (2009). Superfluidity of fermionic bilayers was studied in Ref. Baranov et al. (2011). The bound and scattering states of just two dipoles on different layers have been studied Volosniev et al. (2011); Armstrong et al. (2010) as well.

Our bilayer is realized by a one-dimensional double well potential

(1) |

A DBG in this trap is homogenous and infinitely large in – and –direction and finite in the confinement direction . The dipole moments are aligned along the –direction, therefore the dipole-dipole interaction (DDI) potential takes the form

(2) |

where is the angle between the dipoles and measured from the –axis, and . To stabilize the system against collapse Lahaye et al. (2008) we add a hard core repulsion, that is modeled by a potential. The Hamiltonian describing this many–body system, in the reduced length and energy units, and respectively, looks as follows:

Correspondingly, we give all values for energy, length, wave number, and density in units of , , , and ; hence all quantities in tables and figures are dimensionless.

In previous work O’Dell et al. (2003); Santos et al. (2003); Hufnagl et al. (2010, 2011) it was established that a translationally invariant single layer of a DBG in a one dimensional harmonic trap can become unstable due to the attractive part of the interaction. The pair distribution calculated in Ref. Hufnagl et al., 2011 shows that this instability can be understood as a dimerization, where two dipoles can form a bound state; such weakly bound dipoles would not be stable and indeed experiments with Cr Lahaye et al. (2008) and Er Aikawa et al. (2012) observe a collapse of the BEC for traps that are too wide in the polarization direction. Similarly, two coupled layers can become unstable not only due to the attractive interaction within a layer, but also due to the attraction between dipoles in different layers. This latter “instability” actually indicates the dimerization of dipoles in different layers. As long as the barrier between the layers is high enough that the two bound dipoles remain in their respective layer, such a dimerized phase would be stable.

We have studied the DBG in the trap potential (1) for various potential parameters , , and that control the barrier height between the wells, their separation, and their individual width. We changed the parameters such that we can study the transition from two broad, but well-separated layers to two thin, but close layers. We thereby go from a limit that is dominated by intra-layer attraction to a limit dominated by inter-layer attraction. Both limits are characterized by the appearance of a soft mode (roton) with a respective typical parallel wave number . In the first limit, is the oscillator length of the approximately harmonic well felt by each layer; in the second limit, is the distance between the layers. The six combinations of potential parameters and that we used are listed in the first two columns of table 1; was fixed to . The corresponding trap potentials are plotted in the lower panel of Fig. 1, where we scale it by in order to show all six potentials in the same figure.

We note that, although the average of the DDI over the whole plane vanishes, two dipoles on different two-dimensional planes separated by a distance always dimerize, i.e. form a weakly bound state due to the attractive head-to-tail well of the DDI, regardless of the value of Simon (1976). Hence in the zero density limit, the ground state is always dimerized. At finite density, our calculations in the 2D limit show that dimerization is suppressed by many-body effects unp (). In other words, increasing the density in the two layers stabilizes the monomer phase.

## Ii Ground State

Ground state properties of Bose gases can be calculated using various methods, such as the Gross–Pitaevskii method Gross (1961); Pitaevskii (1961); Dalfovo et al. (1999) and quantum Monte Carlo methods Ceperley (1995); Boronat (1998); Sarsa et al. (2000); Boninsegni et al. (2006). The Gross–Pitaevskii method is widely applied for dilute systems, where the correlations between particles are sufficiently weak such that the interaction between them can be approximated by an effective mean potential felt by each particle. Quantum Monte Carlo on the other hand is also suited for strongly interacting systems, but is computationally demanding. For our calculations of the ground state properties we use the hypernetted–chain Euler Lagrange (HNC-EL) method, which is a variational method suitable for strongly correlated systems Krotscheck et al. (1985), but with lower computational demands than QMC. Starting point is a Jastrow–Feenberg ansatz for the many-body wavefunction

(4) |

which is optimized by solving the Euler-Lagrange equations numerically

Here is the one-body density and is the pair distribution function. and are special cases of the -body density reduced from the full -body probablity of the (normalized) full wave function

Due to translational invariance in and -direction, for the present layer geometry all two-body functions, such as the pair distribution function, depend on three variables: the modulus of the projection of on the plane, , and the two -components and . Hence we effectively have a pair distribution function .

All calculations are done for a total area density of , i.e. for each of the two layers. The HNC equations can be formulated in terms of Mayer cluster diagrams known from classical statistical mechanics Hansen and McDonald (1986). The exact solution of the HNC equations would require the calculation of a class of diagrams called elementary diagrams that cannot be summed exactly. Elementary diagrams are especially important at high densities, while they can be neglected a lower densities. Furthermore, we note that we use a Jastrow–Feenberg ansatz that does not extend to triplet correlations, . From He we know that both elementary diagrams and triplet correlations are important for quantitative agreement with experiment and quantum Monte Carlo simulations Krotscheck (1986). For the area density used in this work, we have checked the influence of the elementary diagrams and the triplet correlations in the 2D limit. In this limit correlations are stronger than for quasi-2D geometries, so the 2D limit gives a conservative estimate of their importance. We found that they improve the accuracy of the static structure function (see below) by less then 2%. Therefore we neglect elementary diagrams and triplet correlations. In Ref. Hufnagl et al. (2010) we compared results for a 2D system of aligned dipoles to QMC calculations Astrakharchik et al. (2007) for the density . The agreement between the HNC-EL results and the QMC calculations for the static structure function is very good. The total energy for both calculations differs by about 2.8%, if elementary diagrams and triplet correlations are neglected. For the even smaller area density considered here, we expect the accuracy to be even better.

In table 1 we show the six parameter combinations that we choose for the trap potential (1) along with the tunnel splitting for a single particle and the chemical potential of the many-body system. is measured with respect to the single-particle ground state energy, i.e. with respect to a non-interacting Bose gas in the same trap potential. As expected, the chemical potential increases if we decrease the thickness of the individual layers, which is due to the intra-layer repulsion of both the DDI and the short-range interaction. With increasing parameter , the two trap wells are not only getting closer, but, with our choice of parameter combinations, also the tunnel splitting increases. As mentioned above we gradually move from thick layers that are widely separated to thin layers that are close to each other, while keeping the total area density fixed at . The results for the density profiles of the DBG in the trap potentials can be seen in the upper panel of Fig. 1. The corresponding trap potentials are shown in the lower panel using the same line style and color (online). is scaled by the inverse of the trap parameter such that all six potentials can be shown using the same scale. Each potential is offset such that the ground state energy of a single particle is zero.

In the limit of zero density or in the non-interacting limit, the density profile is given by the square of the ground state solution, , to the one-body Schrödinger equation . How closely the density of the one-body problem approximates the density of the many-body problem depends on the area density and the strength of the interactions, but also on the strength of the trap potential. For very tight confinement, the eigenenergies of above the ground state energy (or above the first two modes in case of a double well) have energies so high that their contributions to the -body ground state can be neglected. The weaker the confinement, the more and will differ from each other. This is indeed what we find for the three weaker traps with , and . The comparison in Fig. 2 between and shows that the interactions lead to a wider density , that does not agree anymore with the one-body assumption . For the three more confined geometries, and are almost indistinguishable (not shown in Fig. 2). Note that this does not mean that a DDI in a tight trap is well described by the one-body Hamiltonian ; the opposite is true, in-plane correlations are stronger in a tight trap Hufnagl et al. (2011).

For the layer geometry we define a static structure function as

where is the parallel wave number, is the pair distribution function introduced above and is any wave vector in the -plane with wave number . Note that depends only on , while we integrate over the and dependence of . In Fig. 3 we show the static structure function as a function of for the six traps studied here. As we go from well-separated thick layers to close thin layers we observe a peak in in both limits, whereas the peak vanishes in between. It is natural to assume that for wide layers the peak is caused by correlations due to the intra–layer attraction of the dipoles whereas at for a small layer distance it is caused by correlations due to the inter–layer attraction of the dipoles. However, does not contain enough information to distinguish between the these two mechanisms.

In order to gain information about intra- and inter-layer correlations, we define the partially integrated pair densities and ,

(6) | |||||

(7) |

The prefactors are chosen such that for , thus and can be regarded as intra- and inter-layer pair distribution functions. They are the normalized probabilities to find two dipoles in the same layer and in opposite layers at a parallel distance , respectively, regardless of their -coordinate within the layer. and are shown in the lower and upper panel of Fig. 4 for all six traps. For wide, but well-separated layers there are strong intra-layer correlations at , whereas the inter-layer correlations are vanishingly small. This means that two particles in the same layer have a very high probability for head-to-tail configurations, with no parallel separation. As we decrease the thickness of each layer, these intra-layer correlations vanish. At the same time we decrease the distance between layers, thereby increasing the inter-layer correlations. For the smallest distance, two particles in different layers are strongly correlated and have a high probability for head-to-tail configurations. Since the layer is thin, particles in the same layer have a vanishing probability for zero parallel separation because of the DDI and the short-range repulsion. In both limits of two independent wide traps and two close narrow traps the respective strong positive correlations of and suggest a tendency towards dimer formation, where two dipoles align head-to-tail either within a layer or across two layers.

What happens, if we would drive the system to even larger correlation peaks in or ? The instability with respect to dimerization manifests itself as a numerical instability of the HNC-EL equations. Unlike other approximations, the HNC-EL equations have the benefit that they do not produce a solution, if a ground state of an assumed variational form does not exist. In the present case, the Jastrow-Feenberg ansatz (4) does not allow for the dimerization that our above analysis of intra- and inter-layer pair distributions clearly suggests. Since the ground state we try to compute does not exist, our iterative procedure to solve the HNC-EL equations does not converge. In order to actually compute the properties of the dimerized phase, one would have to optimize a variational ansatz that is flexible enough to allow dimerization, or alternatively perform quantum Monte Carlo simulations, as e.g. in Ref. Macia et al. (2012).

## Iii Excitations

### iii.1 Bijl-Feynman modes

Owing to the translational invariance, excitations can be characterized by a parallel (i.e. in-plane) wave number . For a given there are in principle an infinite number of excitations, that are indexed by a perpendicular quantum number associated with out-of-plane motion. Especially for narrow double-well traps these modes have a much higher energy than the lowest two modes in the interesting regime of wave numbers , therefore we restrict our discussion to the two lowest modes and the appearance of a soft mode. In Fig 5 we show the first and the second excitation mode, and in Bijl-Feynman approximation as a function of for the six layer geometries for which we studied the ground state above. The first and the second excitation mode are almost degenerate for a large distance between layers and considerably split for a small distance. Two completely independent DBG layers would of course have two-fold degenerate excitation energies. Even for the most separated layers, the Feynman dispersion is not truly degenerate, but split for small wave numbers. In order to illustrate this, we show the energy difference between the two lowest Feynman energies as a function of in Fig. 6 for the two traps closest to instability. The lifting of the degeneracy can only be due to the DDI that is long ranged and hence couples even well separated layers. We can estimate the typical range of parallel wave numbers for which excitations are most strongly affected. We assume a circular density wave , of wave number , in one layer. A particle at in the other layer will feel a particularly strong dipole force if is such that where is the radius where the DDI changes from attractive to repulsive. is given by , where is the angle of the attractive cone of the dipole interaction, , which gives . From this we get an estimate for the wave vector at which we observe the strongest dipole coupling, which is . If we estimate as the distance between the two maxima of the density profiles shown in the top panel of Fig. 1, we obtain and for the closest and most separated layers, respectively. This simple estimate agrees reasonably well with the maximum energy splitting of the Feynman spectrum at and in Fig. 6. Note that for the DDI averages out, leading to zero DDI-induced splitting for , which is what we observe for well-separated layers. For the closest layers the splitting at is large, however, which is caused by our short-range repulsion model which at such small can be felt between different layers. One could decrease without compromising the stability against intra-layer dimerization, but we preferred to tune only the external trapping potential while keeping the interaction parameters fixed. Furthermore, the tunnel splitting is not small anymore for the closest layers, adding to the splitting caused by the short-range repulsion.

In order to test our conclusions regarding inter-layer coupling for two layers of finite thickness, we also performed calculations for the limit of two 2D layers. In this case the interaction within the same layer is purely repulsive and the attractive part of the interaction is completely missing. Positive correlations are possible only for the inter-layer pair distribution . The 2D results for the two lowest Bijl-Feynman energy dispersions are shown as symbols in Fig. 5. For the wide layers that are far apart, the quasi-2D and 2D results differ substantially (top left panel), which demonstrates that the bending of the dispersion towards forming a roton is an intra-layer effect. As we make each layer narrower, the quasi-2D and 2D results become almost identical. This means that the intra–layer attraction plays less of a role and the roton formation is truly an inter–layer effect. Note that in the 2D limit, the splitting between the two lowest modes vanishes for in agreement with the above argument that the DDI averages out when integrated over the whole 2D plane.

### iii.2 Dynamic structure function from CBF-BW

Calculations of the excitations in the 2D limit of single layers have shown Mazzanti et al. (2009) that the Feynman approximation is adequate for the dispersion relation only at very low densities, but correlation effects become more important as the density is increased and fluctuations of pair correlations must be taken into account. Pair correlation fluctuations are accounted for in the correlated basis function - Brillouin-Wigner (CBF-BW) formalism Krotscheck (1985). The CBF-BW method not only improves the accuracy of the excitation energies, it also describes damping via decay of collective modes into lower energy modes. We will see that the DDI coupling between layers leads to even larger deviations of qualitative features of the excitation spectra in the Feynman approximation.

The CBF-BW method was adapted to layer/film geometries in Ref. Clements et al. (1996) and applied to superfluid He films Clements et al. (1996); Apaja and Krotscheck (2003a, b) and recently to single layers of a DBG Hufnagl et al. (2011). The CBF-BW method has been demonstrated to yield excitation energies much closer to the experimental results than the Feynman approximation, even for such a strongly correlated system as He. Further improvement had been achieved for bulk He by including fluctuations of triplet correlations Campbell and Krotscheck (2010). The added complexity, however, precludes an application to inhomogeneous systems.

The CBF-BW excitation energies are conveniently obtained by following the linear response approach that yields the density-density response operator and – via the fluctuation-dissipation theorem Forster (1975) – the dynamic structure function , where is the frequency of a small external perturbation. The derivation of can be found in Ref. Clements et al., 1996. If we project onto plane waves

we obtain the inelastic cross section for a perturbation imparting the momentum to the system. For a given , a peak in at an energy indicates an excitation of energy . Peaks can have zero linewidth, if decay of an excitation is kinematically forbidden, or finite linewidth otherwise. Translation invariance in the -plane implies that the projection of on the -plane is a good quantum number. A perturbation transferring a parallel momentum and energy to the system probes the dispersion relation of the collective excitations, that we have calculated above in the simpler Feynman approximation. One might think that, since only the parallel component of matters for measuring the dispersion relation, we can restrict ourselves to a parallel , with a vanishing perpendicular component . However, the corresponding dynamic structure function probes only excitation modes of even symmetry with respect to the -plane. Since we are interested not just in the lowest (even) mode but also in the second (odd) mode, we will show also for wave vectors which have an angle with the -plane. A purely perpendicular perturbation () could be implemented by fluctuations of the trapping potential (1) itself, but such a perturbation does not probe the dispersion relation.

#### Parallel Momentum Transfer

In Fig.7 we show for the six different traps shown in Fig.1, the parameters given in table 1 and for wave vectors that are parallel to the -plane, i.e. and hence . is represented in Fig. 7 by mapping to a gray scale. The power of makes sure that also broad, but low-intensity features can be seen well. Full lines track peaks of of zero linewidth, i.e. which are proportional to a -function. The resulting line is an undamped dispersion relation. For sufficiently large wave number , the dispersion relation merges with the gray area, where damping by decay of an excitation into two lower-energy excitations is kinematically possible (i.e. energy and momentum are conserved). The Bijl-Feynman spectrum is shown as dotted lines for comparison, including also higher modes. For wider, well separated traps, the Bijl-Feynman dispersion agrees quite well with the CBF-BW result – even for the widest trap where a roton starts to form due to the intra-layer instability (top left panel). Of course, the Bijl-Feynman approximation does not account for damping. As we confine the two layers more strongly by increasing both trap potential parameters and , the dipole coupling between films leads to a splitting of the Bijl-Feynman energies, as discussed above. has a much richer structure that is poorly represented by the Bijl-Feynman spectrum. On the one hand the density increases as the trap tightens (see Fig.1), and the Bijl-Feynman approximation becomes worse at higher density. On the other hand, the DDI between layers leads, in addition to a splitting of excitation energies, also to more decay channels.

We demonstrate the importance of inter-layer DDI coupling by switching it off for the two traps resulting in the closest layers ( and ). This is simply achieved by setting to zero if and have opposite signs. In Figs. 8 and 10 we show with the full DDI in the left panels and without inter-layer DDI in the right panels. For (Fig. 8), the lack of inter-layer DDI almost completely decouples the two layers, leading to an almost degenerate Bijl-Feynman spectrum. What we get is the dynamic structure function of a single layer, which has been studied in Ref. Hufnagl et al., 2011. The inter-layer DDI leads to significant additional damping for higher energies, seen by the wider peak in the energy regime where the dispersion becomes approximately quadratic. Note that even without inter-layer DDI, there is a bend in the dispersion for , which shows it is not so much caused by inter-layer coupling, but by the intra-layer repulsion of the DDI that, for much higher area density, results in the type of roton studied in Ref. Mazzanti et al., 2009 in the 2D limit.

While a full map is necessary to track the dispersion relation, the detailed line shapes of the various peaks are best seen by plotting slices of for fixed values of . The top panel of Fig. 9 shows, for trap parameters and , a slice of at , which is slightly below the value of where the sharp dispersion curve merges into the damping regime and thus becomes broad (see full map in Fig. 8). The full line and dashed line are the results for with and without inter-layer DDI, respectively. The corresponding excitation energies in Bijl-Feynman approximation are indicated by arrows. The vertical lines are the undamped peaks of the sharp dispersion. We see that has only one broadened peak without inter-layer DDI, while the inclusion of the inter-layer DDI leads to two broadened peaks. We stress again that for parallel momentum transfer, only probes the lower, even mode, hence the two broad peaks are not due to the splitting of a degenerate eigenmode ( for non-parallel momentum transfer is presented below).

As is increased further, the sharp peak loses more and more spectral weight and eventually becomes damped. This case is shown in the lower panel of Fig. 9, where and the sharp dispersion has vanished for both the coupled and uncoupled bilayer, see Fig. 8. Again there is only a single peak without inter-layer DDI and two peaks with inter-layer DDI. The lower peak is caused by the DDI coupling while the higher one is only shifted slightly with respect to its position without DDI coupling. Note that the DDI coupling approximately doubles the width of the higher peak, hence reduces the lifetime of the associated excitation by about a factor of two. Thus, as one can expect, the dipole-dipole coupling between layers leads to faster decay of excitations compared to uncoupled layers.

Finally, in Fig. 10 we compare with and without inter-layer DDI for even closer layers ( and ). The bending is now significantly enhanced by the inter-layer DDI. In CBF-BW approximation, the dispersion (blue line) has a small slope at , i.e. the system is close to “rotonization”. Note that the residual splitting of the dispersion without inter-layer DDI is due to tunneling and the short-range repulsion as mentioned above.

#### Non-Parallel Momentum Transfer

Parallel momentum transfer only probes those excitations which are even with respect to inversion at the -plane, because a perturbation independent of is even and therefore cannot excite odd modes. In order to probe odd modes, we study for wave vectors with an arbitrary angle with respect to the -plane. Fig. 11 shows for , for trap parameters and . We plot as a function of , not , since only is a good quantum number that is meaningful for characterizing the excitation spectrum. Unlike in all previous figures of , we now add an artificial small imaginary part to the energy which slightly broadens all features of . The rationale behind this broadening is that it makes the spectral weight of peaks with zero intrinsic linewidth visible.

The case was shown already in Fig.7 (without artificial damping), and is shown here again for better comparison with . For indeed only the lowest, even mode is visible. As is increased, a second mode becomes visible and gains weight. For , both modes can be seen equally clear in where they appear as a narrow dark trace. Note that the second, odd mode is damped for small . This is very different from the low- behavior of the lowest mode (sound mode) which is not damped because of its negative curvature. The damping of the second mode can be seen as a broadening (in addition to the artificial broadening) for . We introduce the damping limit which is the energy above which an excitation of parallel wave number can decay into two modes with perpendicular wave number and . is given by

where, due to the limitations of the CBF-BW approximation, are the excitation energies in Bijl-Feynman approximation, not the excitation energies following from CBF-BW itself (inclusion of triplet correlations has been shown for homogeneous systems to lead to a self-consistent formulation of the self-energy, see Ref. Campbell and Krotscheck (2010)). Even modes can decay into combinations where is even and vice versa for odd modes. Since we are interested only in decays of the lowest two modes, we obtain three decay limits, which fulfill . They are shown in Fig.11 as dashed lines. The lowest mode can decay into and and the second mode can decay into . The effect of these respective limits are clearly seen in Fig.11. Damping indeed sets in as the dispersion relation of the mode crosses the damping limit with a symmetry appropriate for the mode.

Also visible in Fig.11 are interference patterns that lead to a modulation of the intensity of as and thus is increased. These are simply due to the perpendicular wave number being in phase or out of phase with even or odd modes. For example, a value of , where is a measure for the distance between the layers, leads to a cancellation of the intensity for even modes, but to a maximal intensity for odd modes.

## Iv Discussion and Conclusion

In this work we generalized our previous studies Hufnagl et al. (2011) of dipolar Bose gas layers from a single layer in a harmonic trap to double-well traps which results in a bilayer geometry. As in our previous work, the bilayer calculations are based on the HNC-EL method for the many-body ground state and on the CBF-BW method for excitations. Dipolar bilayers have a richer structure than single layers owing to the inter-layer dipole coupling. The possibility of head-to-tail pairing of two dipoles on different layers leads to similar rotonization effects in the non-paired (monomer) phase as previously predicted in single layers. We restricted ourselves to the calculation of ground state correlations of the monomer phase, as well as its excitation spectrum, including damping due to decay of excitations into two lower excitations.

We systematically varied the double-well trap parameters between two close, but thin layers and two well separated, but wide layers, while keeping the total area density fixed at a modest . The two end points of the range of trap parameters are marked by instabilities of the monomer phase. Either if layers are too close or if one layer is too wide, inter-layer or intra-layer dimerization occurs, respectively. The latter kind of dimers are not stable and would quickly collapse via 3-body collisions, but inter-layer dimers are stable, given a sufficiently high double well barrier. The propensity to pairing was clearly seen in the monomer pair distribution functions or , which are the normalized probabilities to find two dipoles in different or the same layers, respectively, at a parallel distance . We showed that, at , grows a peak for small interlayer distance, while grows one if each single layer is sufficiently wide. In both cases, the peak of is a precursor to the pairing of two dipoles in head-to-tail orientation.

We presented calculation of the dynamic structure function in the CBF-BW approximation. for parallel momentum transfer probes only even modes, where we are mostly interested in the lowest one. typically consists of a lower, undamped peak (that vanishes for higher ) and two broad peaks that are due the inter-layer DDI coupling (without it, there is only one broad peak), which also enhances damping. The double-peak structure is not to be confused with the more trivial effect that each mode is split into two because the inter-layer DDI lifts its degeneracy. The rich structure of is not captured by the simple Bijl-Feynman approximation which would predict a single, undamped peak for the lowest mode. The intra- and inter-layer instabilities of the monomer phase are characterized by a bending of the dispersion relation of the lowest (intra-layer dimer) or lowest two (inter-layer dimer) excitation modes. This bending, that is less pronounced but still visible in the Bijl-Feynman approximation, indicates “rotonization”, which is well-studied for single layers. As in our previous work on single layers, we found that the iterative procedure to solve the nonlinear set of HNC-EL equations becomes unstable as we approach rotonization, i.e. as the dispersion relation starts to have a local minimum at finite . This leads to the conjecture that the ground state is only metastable when the excitation spectrum exhibits a roton, while the true ground state, i.e. the state of lowest energy is the (intra- or inter-layer) dimerized phase. For a proof of this conjecture, however, one would need to compare our monomer results with results for the dimerized phase to find out which state has the lowest energy. Finally, we also presented results for non-parallel momentum transfer, i.e. where the angle between and the plane of the layers is non-zero. The dynamic structure function depends on both the parallel and perpendicular components of , and . still is a good quantum number, while the non-zero allows to probe also odd modes, particularly the second excitation mode. Showing as function of the parallel wave number, the second mode becomes clearly visible for e.g. an angle . Unlike the lowest mode, the second mode is damped for small due to decay into two excitations of lower energy.

An interesting topic are the correlations of the monomer phase and its excitations generalized to layers. The long-ranged DDI coupling between different layers for example lifts an -fold degeneracy of the excitation spectrum and opens many possible decay channels for the resulting modes. Another direction is the study of “unbalanced” bilayers where the two layers have different area densities, or bilayers with different kinds of particles (e.g. different mass) on each layer. If for example the density in one layer is very low, the DDI interaction with the other layer would constitute a very well controlled model of an impurity particle moving in one layer, coupled to a bath of particles in the other layer. The mechanism how an impurity attains an effective mass could be investigated in a well-controlled fashion over a much wider range of densities and interaction strengths than in condensed matter systems.

###### Acknowledgements.

We are grateful for discussions with Eckhard Krotscheck and Vesa Apaja. We acknowledge financial support by the Austrian Science Fund FWF under grant No. #23535 and by the National Science Foundation under Grant No. NSF PHY11-25915.### References

- T. Lahaye, T. Koch, B. Fröhlich, M. Fattori, J. Metz, A. Griesmaier, S. Giovanazzi, and T. Pfau, Nature 448, 672 (2007).
- T. Lahaye, J. Metz, B. Fröhlich, T. Koch, M. Meister, A. Griesmaier, T. Pfau, H. Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. Lett. 101, 080401 (2008).
- M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Phys. Rev. Lett. 107, 190401 (2011).
- K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401 (2012).
- K.-K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Pe’er, B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne, D. S. Jin, and J. Ye, Science 322, 231 (2008).
- S. Ospelkaus, A. Pe’er, K.-K. Ni, J. J. Zirbel, B. Neyenhuis, S. Kotochigova, P. S. Julienne, J. Ye, and D. S. Jin, Nature Physics 4, 622 (2008).
- K. K. Ni, S. Ospelkaus, D. J. Nesbitt, J. Ye, and D. S. Jin, 11, 9626 (2009).
- K.-K. Ni, S. Ospelkaus, D. Wang, G. Quemener, B. Neyenhuis, M. H. G. de Miranda, J. L. Bohn, J. Ye, and D. S. Jin, Nature 464, 1324 (2010).
- J. Deiglmayr, A. Grochola, M. Repp, K. Mörtlbauer, C. Glück, J. Lange, O. Dulieu, R. Wester, and M. Weidemüller, Phys. Rev. Lett. 101, 133004/1 (2008).
- A.-C. Voigt, M. Taglieber, L. Costa, T. Aoki, W. Wieser, T. W. Hänsch, and K. Dieckmann, Phys. Rev. Lett. 102, 020405 (2009).
- J. M. Sage, S. Sainis, T. Bergeman, and D. DeMille, Phys. Rev. Lett. 94, 203001 (2005).
- T. Takekoshi, M. Debatin, R. Rameshan, F. Ferlaino, R. Grimm, H.-C. Nägerl, C. R. Le Sueur, J. M. Hutson, P. S. Julienne, S. Kotochigova, and E. Tiemann, Phys. Rev. A 85, 032506 (2012).
- K. Goral, K. Rzazewski, and T. Pfau, Phys. Rev. A 61, 051601 (2000).
- S. Ronen, D. C. E. Bortolotti, , and J. L. Bohn, Phys. Rev. Lett. 98, 030406 (2007).
- T. Koch, T. Lahaye, J. Metz, B. Fröhlich, A. Griesmaier, and T. Pfau, Nature Physics 4, 218 (2008).
- D. H. J. O’Dell, S. Giovanazzi, and G. Kurizki, Phys. Rev. Lett. 90, 110402/1 (2003).
- L. Santos, G. V. Shlyapnikov, and M. Lewenstein, Phys. Rev. Lett. 90, 250403/1 (2003).
- F. Mazzanti, R. E. Zillich, G. E. Astrakharchik, and J. Boronat, Phys. Rev. Lett. 102, 110405 (2009).
- D. Hufnagl, R. Kaltseis, V. Apaja, and R. E. Zillich, Phys. Rev. Lett. 107, 065303 (2011).
- A. Macia, D. Hufnagl, F. Mazzanti, J. Boronat, and R. E. Zillich, Phys. Rev. Lett. 109, 235307 (2012).
- G. Bismut, B. Laburthe-Tolra, E. Marechal, P. Pedri, O. Gorceix, and L. Vernac, Phys. Rev. Lett. 109, 155302 (2012).
- A. Micheli, G. K. Brennen, and P. Zoller, Nature Physics 2, 341 (2006).
- M. A. Baranov, M. Dalmonte, G. Pupillo, and P. Zoller, Chem. Rev. 112, 5012 (2012).
- F. Ferlaino, private communication.
- H. P. Buchler, E. Demler, M. Lukin, A. Micheli, N. Prokofev, G. Pupillo, and P. Zoller, Phys. Rev. Lett. 98, 060404/1 (2007).
- G. E. Astrakharchik, J. Boronat, I. L. Kurbakov, and Y. E. Lozovik, Phys. Rev. Lett. 98, 060405 (2007).
- A. Filinov, N. V. Prokofâev, and M. Bonitz, Phys. Rev. Lett. 105, 070401 (2010).
- A. Filinov and M. Bonitz, Phys. Rev. A 86, 043628 (2012).
- A. Macia, F. Mazzanti, J. Boronat, and R. E. Zillich, Phys. Rev. A 84, 033625 (2011).
- R. M. Wilson, S. Ronen, J. L. Bohn, and H. Pu, Phys. Rev. Lett. 100, 245302/1 (2008).
- E. Krotscheck, G.-X. Qian, and W. Kohn, Phys. Rev. B 31, 4245 (1985).
- D. Hufnagl, E. Krotscheck, and R. E. Zillich, J. of Low Temp. Phys. 158, 85 (2010).
- D.-W. Wang, Phys. Rev. Lett. 98, 060403 (2007).
- C. Trefzger, C. Menotti, and M. Lewenstein, Phys. Rev. Lett. 103, 035304 (2009).
- M. A. Baranov, A. Micheli, S. Ronen, and P. Zoller, Phys. Rev. A 83, 043602 (2011).
- A. G. Volosniev, N. T. Zinner, D. V. Fedorov, A. S. Jensen, and B. Wunsch, J. Phys. B 44, 125301 (2011).
- J. R. Armstrong, N. T. Zinner, D. V. Fedorov, and A. S. Jensen, Eur. Phys. Lett. 91, 16001 (2010).
- B. Simon, Annals of Physics 97, 279 (1976).
- Unpublished.
- E. P. Gross, Nuovo Cim. 20, 454 (1961).
- L. P. Pitaevskii, Sov. Phys. JETP 12, 451 (1961).
- F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999).
- D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
- J. Boronat, in Microscopic Quantum Many-Body Theories and Their Applications, Lecture Notes in Physics, Vol. 510, edited by J. Navarro and A. Polls (Springer, 1998) pp. 359–379.
- A. Sarsa, K. E. Schmidt, and W. R. Magro, J. Chem. Phys. 113, 1366 (2000).
- M. Boninsegni, N. V. Prokofev, and B. V. Svistunov, Phys. Rev. E 74, 036701 (2006).
- J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, 1986).
- E. Krotscheck, Phys. Rev. B 33, 3158 (1986).
- E. Krotscheck, Phys. Rev. B 31, 4258 (1985).
- B. E. Clements, E. Krotscheck, and C. J. Tymczak, Phys. Rev. B 53, 12253 (1996).
- V. Apaja and E. Krotscheck, Phys. Rev. Lett. 91, 225302 (2003a).
- V. Apaja and E. Krotscheck, Phys. Rev. B 67, 184304 (2003b).
- C. E. Campbell and E. Krotscheck, J. of Low Temp. Phys. 158, 226 (2010).
- D. Forster, Hydrodynamics, Fluctuations, Broken Symmetry, and Correlation Functions, Frontiers in Physics, Vol. 47 (W. A. Benjamin, 1975).