A G Volosniev

N T Zinner

D V Fedorov

A S Jensen

B Wunsch

Bound dimers in bilayers of cold polar molecules

Keywords: numbers: 67.85.-d, 36.20.-r, 05.30.-d

Bound dimers in bilayers of cold polar molecules

a G Volosniev, N T Zinner, D V Fedorov, a S Jensen,

Introduction

A strong experimental drive in the field of polar atoms and molecules has realized controllable samples in the rotational and vibrational ground-state that are close to quantum degeneracy [1,2,3,4,5,6]. These heteronuclear systems have a number of very interesting properties due to the long-range and anisotropic dipole-dipole force which can give rise to highly non-trivial many-body states in both the weak-and strongcoupling regime [7,8]. The attractive head-to-tail configuration can, however, lead to strong chemical reactions [6] or many-body collapse of the system [9], and confinement in optical lattices has been suggested as a means of avoiding this problem [10]. These confined one-or two-dimensional geometries have led to a number of predictions of novel few-and many-body states [10,11,12,13,14,15,16,17,18,19,20,21,22], and very recently the first experimental implementation of a multilayered stack of pancakes containing fermionic polar molecules was reported [23].

Here we consider the case of two adjacent layers. However, even in this seemingly simple case there is a competition of intra-and interlayer interactions which can vary between repulsion and attraction as one changes the polarization angle of the dipole moments with respect to the layers. In the present paper we will be concerned with few-body states with one particle in each layer in order to describe the simplest complex in such a system in detail. The case of dipoles oriented perpendicular to the layers was considered from the few-body bound state and scattering point of view in previous works [24,25,26,27]. At the so-called 'magic' angle where the intralayer repulsion vanishes in a one-dimensional trap the few-body bound state structure was also discussed [28,29,30].

To our knowledge, the full two-body bound-state problem as a function of the polarization angle and the dipole moment has not been studied previously. This problem is highly non-trivial due to (i) the anisotropy and (ii) the vanishing integral over space of the potential for arbitrary polarization angle. The problem is a specification of a more general problem of two particles in two dimensions interacting via anisotropic potentials where the net volume is negative or zero. The more general problem has been addressed only in a recent letter [31], where the details of analytical and numerical methods are not elaborated.

In this paper, we present in section 2 analytic results for the dipole-dipole potential. They are specialization of the derivation in Appendix A for an arbitrary potential. In section 3, we describe the novel numerical method based on stochastic variation along with the computed results for the dipole-dipole potential. We present energies, wave functions, and expectation values of relevant operators as functions of the potential strength and the polarization angle. One of our main results, analytical as well as numerical, is that the bilayer system has a bound state for any polarization angle and any value of the dipole moment. We also calculate a partial-wave decomposition that characterizes the geometric structure of the wave function which indicates the likely symmetries of the corresponding many-body problem. In section 4 we present a first application of our results in a many-body context. We consider the limit of strong coupling where the system forms bound bosonic dimers that can potentially form a (quasi)-condensate. Finally in section 5 we brifly summarize and conclude.

Analytic results

The general setup we consider consists of two particles confined to two spatial dimensions with a pair potential V ( r) depending on the relative coordinate r. In polar coordinates, r = (x, y) = (r cos ϕ, r sin ϕ), we have the Schrödinger equation:

- 1 s ∂ ∂s s ∂ ∂s - 1 s 2 ∂ 2 ∂ϕ 2 + λ V (s, ϕ) Ψ = α 2 Ψ , ( 1
)
where Ψ is the wave function, µ the reduced mass, λ is the dimensionless strength, λ V = 2µd 2 V / 2 and α 2 = 2µd 2 E/ 2 , V is the potential and E the energy, d the unit of length, and s = r/d, is the reduced coordinate.

The investigation of possible bound states was briefly sketched analytically in [31] for an arbitrary potential in two spatial dimensions (2D). The result is that any cylindrical or non-cylindrical potential has a least one bound state provided the volume of the potential is negative or zero. We give more details of a similar derivation in Appendix A. The difficulties are centered around exceedingly weak potentials where no ordinary perturbation treatment is applicable because no unperturbed solution is available. The binding energy approaches zero as the potential vanishes and the continuum is approached. Thus the limiting energy is zero and the corresponding wave function is uniformly distributed over all space.

We specialize to the dipole-dipole potential arising for the system of two polarized molecules of mass M confined to two parallel planes separated by a distance d as shown in figure 1. The corresponding dipole-dipole potential, V , projected to this two-dimensional geometry is

V (r, ϕ) = D 2 r 2 + d 2 -3(r cos ϕ cos θ + d sin θ) 2 (r 2 + d 2 ) 5/2 , ( 2
)
where D is the dipole moment ‡, and θ denotes the polarization angle measured from the layer plane to the z-axis which intersects the two layers at right angles. The potential in equation ( 2) is found in the ideal limit of zero-width layers. The dipole polarization is measured such that for θ = π/2, the dipoles are oriented perpendicular to the layers as in [24,26,27]. One can take corrections to the zerowidth layer limit into account by integrating out a gaussian in the transverse direction. However, the corrections are second-order in the width, w, and we neglect them as we here are interested in the w d limit. We have to solve the 2D Schrödinger equation in equation (1) with the potential in equation (2). The reduced mass is µ = M/2 and the dimensionless dipolar strength is given by U = MD 2 /( 2 d), which is a measure of the ratio of potential to kinetic energy. We will also consider the case where U < 0 which is also physically realizable as explained below. In the notation above we have λ = U and we will use these notations interchangeably in order to emphasize the generality of the analytic approach presented in Appendix A. This potential is invariant under reflection in the x-axis and has the peculiar property that d xd yV (x, y) = 0 for any θ. In particular, it does not fulfil the Landau criterion for bound states in two dimensions [START_REF] Landau | Quantum Mechanics[END_REF], which states that a bound state always exists for d xd yV (x, y) < 0. An early existance proof was given in [START_REF] Simon | [END_REF] using a method that is not well-suited for expansions in the strengh of the potential. A discussion of such an expansion appeared in [27] but only for case where θ = π/2 and cylindrical symmetry holds. Here we are interested in the appearance and properties of bound states for arbitrary θ. A partial-wave decomposition of the potential in the basis {1, cos ϕ, cos(2ϕ)} (which are the only non-zero terms) leads to

λ V (s, ϕ) = V 0 (s) + V 1 (s) cos ϕ + V 2 (s) cos(2ϕ) ( 3 ) V 0 (s) = U [3 sin 2 θ -1][s 2 /2 -1] (s 2 + 1) 5/2 , ( 4
) V 1 (s) = -3U s sin(2θ) (s 2 + 1) 5/2 , (5)
V 2 (s) = - 3 2 U s 2 cos 2 θ (s 2 + 1) 5/2 , (6)
which we will refer to as monopole, dipole, and quadrupole terms, respectively. The monopole potential V 0 , has in itself zero net volume and it vanishes identically for θ = θ c where sin 2 θ c = 1/3. The dipole term only vanishes for θ = 0 and π/2, whereas the quadrupole term is finite except at θ = π/2. Thus for θ > θ c and U > 0, the monopole term has an inner attractive pocket and a repulsive barrier outside s = √ 2, and vice versa for θ < θ c . For U < 0 the story is reversed. We expect the monopole term to be most important for the system properties, at least when it is non-vanishing away from θ = θ c . However, the monopole term is, except for the factor of (3 sin 2 θ -1)/2, identical to the full potential at π/2, i.e. we know from previous work that it always supports bound states [26,27,[START_REF] Simon | [END_REF]. We also know that the configuration with an attractive inner pocket and a repulsive outer barrier leads to considerably stronger binding than in the reversed case [26]. We will see this explicitly in the energies presented below.

It is very important to notice that the angle θ c is different from the magic angle, θ * c , where the potential of two dipoles moving in one dimension vanishes (determined by cos 2 θ * c = 1/3) [8]. This demonstrates an important difference between one-and twodimensional dipolar systems. We will address this fact in more detail when we discuss many-body physics below.

First we specialize the analytical results in Appendix A derived for general interactions to the dipolar potential. A partial wave expansion is employed in analogy to that of the decomposition in equation ( 3). For the dipole potential in equation ( 2) the resulting energy expression from equation (A.25) then becomes

E = - 4 2 Md 2 exp -2γ - 2(1 + UB 1 + U 2 B 2 ) U 2 (A 0 + UA 1 + U 2 A 2 ) , (7)
where γ is Euler's constant and the coefficients, A 0 , A 1 and B 1 , B 2 , are defined by

A 0 = 1 4 M 2 c + 1 8 sin 2 2θ + 1 32 cos 4 θ , (8)
A 1 = + 0.0053 sin 2 2θ cos 2 θ -0.0033 sin 2θ cos 4 θ -0.0019 cos 6 θ -M c (0.0349 sin 2 (2θ)

+ 0.0054 cos 4 (θ) + 0.0156M c cos 2 θ + 0.0343M 2 c ) ,

B 1 = -1.204M c - 1 16 cos 2 θ , (10)
B 2 = 0.8382M c (M c + 0.0667 cos 2 θ) ( 1 1 )
-0.0037 sin 2 (2θ) + 0.0894 cos 4 θ ,

M c = 3 2 sin 2 (θ) - 1 2 . ( 12
)
The expression for A 2 is much more elaborate consisting of more than a hundred terms each given as double, triple, and quadrupole integrals over well defined functions. We refrain from showing them all here. Note that since the spatial integral of the potential vanished, the term A 2 has to be considered in an expansion to second order in U (see also (A.25)).

The only non-zero matrix elements for the dipole potential in equation ( 2) are V m,m , V m,m±1 , V m,m±2 . This implies that the m = 0, 1, and 2 partial waves are sufficient to get all contributions to the orders given on the right hand side of equation (7). Higher partial waves beyond m = 2 only contribute to the wave function for the dipole potential through orders of U that are higher than those in equation 7.

The scattering length, a, is usually considered to be a measure of the most crucial model independent property of any potential. In the present case it is a function of strength and polarization angle, and to any order related to the energy as in [34], i.e. through the general equation

E = -2 2 µa 2 exp(-2γ) . (13)
Any accuracy of E is then directly transferable to a through equation (13). The energy can be calculated to any order in powers of U , as second order in equation (7). Then the scattering length becomes

a d = exp (1 + UB 1 + U 2 B 2 ) U 2 (A 0 + UA 1 + U 2 A 2 ) , (14)
where the dependence on strength and polarization angle now is explicit.

The energy, and scattering length, very close to threshold is exponential in U -2 , as seen in eqs.( 7) and ( 14), and determined by the polarization angle through A 0 . The first order terms, (A 1 , B 1 ), in U exhibits the difference in approach to threshold for the different signs of the strength, U . The second order terms, (A 2 , B 2 ), are necessary to get the correct U -independent pre-exponential factor in the energy. Here B 2 is both much simpler and more significant than A 2 which consists of sums of a large number of contributions expressed as definite integrals.

The expressions simplify substantially for θ = π/2. In [35] the energy is calculated for θ = π/2 to the order given in equation (7), including the A 2 term, in agreement with our result. The computation of the energy in [27] for θ = π/2 deviates from our results in equation (7) in the first order correction.

Numerical procedure

The potential is in general anisotropic and the wave equation is not easy to solve by discretization or integration. We therefore turn to the stochastic variational approach using gaussian wave functions which has been successfully applied to other interactions [START_REF] Suzuki | Stochastic Variational Approach to Quantum-Mechanical Few-Body Problems[END_REF]. However, in the limit of weak binding the wave functions become very small and spatially extended without structure at large distances. The special method to achieve convergence with a fair amount of gaussians is described in this section, and the results for energies and wave functions presented in the next section.

Method for weakly bound systems

For numerical calculations we employ the correlated Gaussian method which has been sucessfully used in a range of few-body problems in atomic physics [START_REF] Thøgersen | [END_REF]38,39]. The wave function Ψ(x, y) is found through the variational principle as a linear combination of basis-functions G i (x, y),

Ψ(x, y) = N basis i=1 c i G i (x, y), (15)
where N basis is the size of the basis, c i are the linear variational parameters, and the basis-functions are chosen in the form of shifted correlated Gaussians,

G i (x, y) = e -(q-s i ) T A i (q-s i ) , (16)
where the superscript T denotes transposition, q ≡ (x, y) T is the column-vector of the coordinates, and where the elements of the symmetric positive correlation matrices A i and the shift vectors s i are the non-linear variational parameters. The explicit shifts employed here enhance greatly the flexibility of the correlated Gaussians specifically for the non-rotationally symmetric system at hand. According to the variational principle the wave function is found by minimizing the expectation value of the Hamiltonian,

E[Ψ] ≡ Ψ|H|Ψ Ψ|Ψ . ( 17
)
The linear parameters c i are determined by solving the generalized eigenvalue problem

Hc = E[Ψ]N c , (18)
where H and N are N basis × N basis matrices,

H ij ≡ G i |H|G j and N ij ≡ G i |G j , (19)
where c is the column-vector of the linear parameters c i . The overlaps and the matrix elements of the kinetic energy operator in equation ( 19) are calculated using the analytical expressions in Appendix B. The matrix elements of the interaction potential are calculated numerically using adaptive Gauss-Kronrod quadratures [40].

The generalized eigenvalue problem equation ( 18) is solved by first performing Cholesky decomposition of the (symmetric and positive definite) matrix N = LL T , and then solving the ordinary symmetric eigenvalue problem,

Hc = E[Ψ]c , (20)
where H = L -1 HL -1 T and c = L T c, using standard linear algebra methods [40].

The non-linear parameters -the elements of the correlation matrices A i and the shift vectors s i -are determined through stochastic sampling. Basically, the minimum of the energy is found by sampling a large number of random sets of N basis Gaussians with randomly generated correlation matrices and shift vectors. The elements of the correlation matrices are generated as ±1/b 2 and the elements of the shift vectors as ±b, where b is a stochastic variable with the dimension of length sampled from an exponential distribution with certain length parameter l. To ensure the positive definiteness of the correlation matrix the following procedure is used: First a diagonal matrix with positive elements 1/b 2 is generated, and then a random orthogonal transformation is performed.

For ordinary quantum systems the length parameter, l, of the exponential distribution should be chosen close to the typical range, d, of the interaction potential between particles. For the dipolar potential considered in this paper, it is in fact the interlayer distance d which explains the abuse of notation. If the binding energy B of the system is of natural size, that is B ∼ weakly bound systems, where B 2 2µd 2 , two different length scales are of importance: One is the interaction range d and the other is the typical size |α| -1 of the tail of the wave function (where α = i 2µB

2 ). We therefore subdivide the basis into N short , shortrange, and N long , long-range, Gaussians, N basis = N short + N long , where the short-range Gaussians are sampled from the exponential distribution with range d and the longrange Gaussians are sampled from the exponential distribution with range |α| -1 . The latter is calculated from the current estimate of the binding energy. Since the long-range Gaussians are introduced specifically to better describe the asymptotics of the wave function, they can be chosen in a much simpler form,

U = 1.8; θ = π/4 N basis EM d 2
G = e -r 2 /b 2 , r ≡ x 2 + y 2 . ( 21
)
The typical convergence plot of the binding energy as function of the basis size is presented in figure 2. It shows that the energy is converged to within four significant digits with the basis size of about 45 Gaussians. The convergence of the tail of the wave function is illustrated in figure 3 where it is compared to the analytic asymptotic form. Clearly, addition of the long-range Gaussians significantly improves the quality of the wave function at very large distances. In this example about 46 Gaussians, 30 short-ranged plus 16 long-ranged, can accurately describe the the wave function up to the radius of about 800d.

Results

The general derivation in Appendix A demonstrates that there always is at least one bound state in 2D for any anisotropic potential with zero net volume as obtained in a different manner in [START_REF] Simon | [END_REF]. When U becomes small we expect universal behavior of 30 + 8 30 + 12 30 + 16 const × K 0 (|α|s) energies and radii [31,34,41]. Using the stochastic variational approach in the small U limit, our results for θ = π/2 approach the universal behavior of the energy which to leading order scales like ln(-E) ∝ -1/U 2 as discussed previously [26,27,31,[START_REF] Simon | [END_REF]. For other values of θ we expect the same scaling for very small U , however, the range of U around zero where this applies has a strong dependence on θ as seen by the energies presented below. We also expect differences for general θ between positive and negative U in the limit U → 0 as for θ = π/2 [26]. The question of how the binding energy approaches universality is investigated in more details in [31].

U = 0.9; θ = θ * c s = ( x 2 + y 2 )/d A 0 Φ 0 (s) √
The energies have been calculated using the correlated gaussian approach. In figure 4 we exhibit the results as a function of U > 0 for a selection of polarization angles. At small U the energy decreases very fast with decreasing U as noted already in [START_REF] Simon | [END_REF], whereas at larger U we find a linear dependence on U as argued in [26] for θ = π/2. The binding energies decrease dramatically as θ approaches zero. We stress that we numerically find a bound state for any value of U also in the particularly unfavorable case of θ = 0 as well as for θ c (sin 2 θ c = 1/3) where the decisive monopole potential vanishes identically. Notice that there are more bound states for larger U but we restricted our discussion to the single bound state regime.

The U < 0 case is also of great interest as that potential can be generated by using microwave-dressed molecules. In [15,16,22] an AC light field directed perpendicular to the layers was used to create the θ = π/2 potential with U < 0. A straightforward calculation shows that if the laser hits the layers at an angle, θ, the potential is the same as for a homogeneous electric field at angle θ but with negative U . For U < 0 we again find numerically that for all values of the strength the two-body system has bound states. The results for the binding energy at different angles are shown in figure 5 as function of |U |. The first thing one notices is that the overall magnitude of the bound state energy is smaller than that for U > 0. At θ = π/2 this can be understood as the potential has a repulsive core at θ = 0, forcing the state to reside in the shallow attractive pocket at intermediate distance. In turn, this gives a much smaller binding energy.

U M d 2 2h 2 E y -2 0 2 -3 -2 -1 0 1 2 3 x y -2 0 2 -2 0 2 x y -2 0 2 -2 0 2 θ = θ c θ = θ c θ = 0 θ = π 2 θ = π 4 θ = π 2 θ = 0
This qualitative behavior of the potential persists until θ decreases below θ c where the monopole changes sign. Then the potential changes overall character to become more attractive with inner attractive pocket and outer repulsive tail. The U < 0 results thus show maximum binding at θ = 0 which is, however, still about a factor of three smaller than the U > 0 case at its most favorable angle of θ = π/2. Numerical and analytical results are compared in figure 6. To get the most accurate comparison we show the exponent in equation (7) multiplied by the square of the strength, U 2 . The approach in the figure of analytical and numerical results is reassuring in the limit of relatively weak potentials. Since we neglected the A 2 -contribution this approach indicates its minor significance.

Wave functions

The structure of the bound state wave functions can be seen from the partial wave decomposition. The results are shown in figure 7 for a strong coupling of U = 10 and a weaker one of U = 4. The probabilities are normalized so that they sum to one. We note that the contribution of m > 2 is only a few percent with a maximum at θ = 0 of 5% in m > 2 terms. As expected we find that m = 0 becomes dominant for θ → π/2 as we approach cylindrical symmetry. Interestingly, close to θ = 0 we also find a very large m = 0 component, no m = 1 content, and a significant m = 2 contribution. The remaining content of the wave function is found in the higher m contributions. The fact that m = 1 has no weight for θ = 0 can be understood from the symmetry of the potential. For x → -x the m = 1 term changes sign, whereas the potential is invariant. Interestingly, as we go away from θ = 0, the m = 1 component raises rapidly and stays on the order of 40% until we reach θ = θ c at which it starts to decline as for m = 2, in line with the restoration of cylindrical symmetry at θ = π/2. For U > 10 the m = 0 component can be even more suppressed in comparison to m > 0 for intermediate θ, whereas for positive U < 4 the m = 0 component will eventually dominate as one approaches the universal limit discussed above. )) for different polarization angles θ. The dashed red curves are calculated numerically and the solid blue curves are from equation (7). The relatively small contribution from the A 2 -term is neglected in this comparison.

We have found similar results for the U < 0 when taking into account that the angle θ for U > 0 correspond to angle π/2 -θ for U < 0 and vice versa. This is in fact an exact symmetry of the dipole part of the potential and an approximate one for the monopole term since θ c is close to π/4. For U = -10 we find that there is a window θ c < θ < 1.1 in which the m = 1 term is around 40%. Interestingly, we find that the partial-wave content for U < 0 is almost exclusively m = 0 and m = 1. This is perhaps surprising as the potential in the m = 2 channel is non-vanishing except at θ = π/2.

The potential has completely different forms for different polarization angles as illustrated on the contour plots in the inset of figure 4. For θ = π/2, the potential is cylindrical while asymmetry appears for decreasing θ-values, and eventually two minima emerge when θ approaches zero. For small θ, the potential looks like a harmonic oscillator along the y-axis for small x and y. The depth of this harmonic well around zero is about twice as large for θ = 0 and U < 0 in comparison to the depth of the two wells in the x-direction for θ = 0 and U > 0 that is shown in the inset of figure 4 and figure 8. A simple gaussian wave function should therefore be a fair approximation to the full problem and in turn the lowest partial-waves should dominate.

The wave functions for strongly bound states mimic the contours of the potential. As the strength, as well as binding energy, decreases the wave function is spreading out to larger distances and approaching cylindrical symmetry, as illustrated in figure 8. The probability decreases in all points of space and approaches zero uniformly outside the potential. However, inside the potential the shape of the potential is maintained even for vanishingly small strengths where the probability also approaches zero. This behavior is necessary to provide binding which arises from the attractive potential at small distances. In turn, the modified Bessel function, K 0 (|α|s), is approached for vanishing strength which corresponds to a wavefunction that is roughly constant in space until the distance, s, is comparable to 1/|α|. V (x, y)

The Many-Body Bilayer System

The bilayer system has an interesting many-body structure with combination of attractive interactions that can induce pairing and repulsive interaction that tend to suppress such effects. This was discussed recently for the θ = π/2 case in [19,20]. Here we consider the strongly-coupled limit (large U ) where the bound two-body dimers are expected to be the relevant degrees of freedom. As the dimers are effectively bosons, they are capable of forming a (quasi)-condensate under the right conditions [20]. However, as is well-known from BCS-BEC crossover studies [42], this is only expected to occur when the density is low. Unfortunately, the Berezinskii-Kosterlitz-Thouless (BKT) transition [43,44] that governs this two-dimensional system has a critical temperature that is proportional to density [19,20]. Therefore, a compromise where the dimer condensate occurs at not too low densities would be optimal to allow experimental access to this unusual many-body state.

A further interesting complication is the fact that each of the layers can accomodate various coherent many-body states when we consider them independently. Proposals for the ground state include p-wave superfluids [14,16] and density waves [17,18]. Of course, when more than one layer is present the long-range inter-and intralayer interactions are competing, and at this point it is not clear what state is favoured for arbitrary directions of the polarization. It is known that the density-wave instability will be enhanced and occurs at smaller coupling strength for bi-and multilayer systems [21]. We note that in the U < 0 bilayer case, a particle-hole coherence reminiscent of ferromagnetism has even been suggested [15]. Here we will be concerned mainly with the BCS-BEC crossover scenario in the strongly-coupled limit, but we will also estimate the appearance of a density-wave state and comment on possible superfluids.

As the criterion for the onset of condensation of dimers, we consider the point at which the chemical potential becomes negative [20], i.e.

µ(U, θ) = 1 2 nV eff (U, θ) + E F - 1 2 E B (U, θ), (22)
where E B is the dimer binding energy and V eff is the long-wavelength (zero momentum) effective momentum-space interaction between two dimers. Here we include both the binding energy and the dimer-dimer interaction, and we also include a term for the Fermi energy, E F , that the constituents of the dimer inherit from their layer. The density of dimers (equal to the single-layer density when the layers have an equal number of molecules) is denoted by n. To obtain the effective interaction, one must in principle integrate out the wave function of the dimer and include all inter-and intralayer twobody terms [20]. However, here we are only interested in the long-wavelength limit (momentum zero) in which the interlayer term vanishes [13]. This gives

V eff (U, θ) = 2 M 4U 3 √ 2π d w 4πP 2 (sin θ), (23)
where P 2 (x) = (3x 2 -1)/2. For the layer width, we take w/d = 0.2 in the following. Notice that V eff is attractive for θ < θ c , vanishes at θ c , and repulsive for θ > θ c . The attraction for θ < θ c results in a negative compressibility in a single layer [14]. We stress again that θ c is much smaller than the angle at which the intralayer repulsion vanishes in a one-dimensional system. In this sense θ c is a special angle for the intralayer repulsion, whereas it has no dramatic effect on the binding energies which vary smoothly around θ = θ c . Combining the formula above, the final expression for µ becomes

Md 2 2 µ = (k F d) 2 2 4U 3 √ 2π d w P 2 (sin θ) + 1 - Md 2 2 2 E B , (24)
where we use the Fermi momentum k 2 F = 4πn for fermions in a single layer in place of n.

The lines of µ = 0 for selected angles are shown in figure 9 in the (U, k F d) plane for 1.5 < U < 5. For U > 2 the dimers have significant binding energy and can we treat them as localized bosonic objects. For θ = π/2 we present results both with and without the intralayer term which is clearly seen to shrink the region of potential dimer condensation. For θ = θ c , the intralayer term vanishes and we find a larger region of µ < 0. For θ < θ c the region would in principle become even larger, however, the intralayer term is then attractive and can lead to instability and collapse [14]. We therefore expect the line for θ = θ c to provide a boundary for how large the BEC region can become when tuning the angle within our approximations. In figure 9 we also show the line above which the density-wave instability appears at angle θ = θ c in the randomphase approximation as discussed in Ref. [21] (above full blue line, denote by DW). To study the crossover outside the density-wave region, we see that low densities are indeed needed.

We now consider the important question of the finite temperature behavior of the system. In the large U limit, the BKT transition temperature is maximal at k B T BKT = E F /8 [15,19,20], or

T BKT = 765 n 10 8 cm -2 amu M nK. ( 25
)
If we consider LiCs molecules (which can have a dipole moment of up to 5.5 Debye) with d = 0.5 µm, then the µ = 0 phase for θ = θ c at U = 10 is at n = 8.1 • 10 7 cm 2 , and thus T BKT ∼ 4.9 nK. While still very low, this is a significant increase over the sub-nano Kelvin temperatures for θ = π/2. The p-wave superfluid for single layers is estimated to appear for angles sin θ < 2/3 [14]. However, density-waves occur for sin 2 θ > 1/3. This implies that there is a rather large regime (2/3 < sin θ ≤ 1) where the single-layer p-wave superfluid should be absent but where BCS-BEC crossover induced by the interlayer interaction is possible. We also note that there is an intermediate region where both p-wave superfluid and densitywave have been proposed as the ground state and one could imagine something like a supersolid phase as recently discussed in a two-dimensional optical lattice with fermionic polar molecules and perpendicular polarization [45]. In the strong-coupling limit we expect that the allowed few-body states will play an important role in determining the system properties, and calculation of the binding energies of states with three or more particles for arbitrary directions of the polarization are being pursued. . Lines of vanishing chemical potential for sin 2 θ c = 1/3 (solid black) and θ = π/2 with (solid red) and without (dashed red) intralayer repulsion. Above the solid lines we expect condensation of dimers to occur (BEC), whereas in the lower right part a many-body paired BCS-like state (BCS) should be the ground-state of the system. Also shown is the line above which a density wave (DW) for θ = θ c is expected in the bilayer system (solid blue).

We expect that the partial-wave analysis presented in the previous section can help indicate what symmetries are possible and relevant for the corresponding manybody problem in the large U limit. The problem is of course still that the intralayer term is attractive in the long-wavelength limit for θ < θ c , and we thus expect that the most stable system require θ > θ c where the decomposition of the wave function is entirely dominated by the m = 0 term. However, close to θ c we still have a substantial m = 1 contribution (immediately to the right of the vertical line in figure 7). We therefore expect a region of interest in which an exotic many-body state with nontrivial symmetry like a p-wave dominated or mixed symmetry superfluid would emerge in the bilayer. These indications are consistent with a partial-wave decomposition of the potential itself [46]. Combining this information strongly suggests that there is a very interesting crossover from weak-to strong-coupling in the corresponding many-body system as recently discussed for the θ = π/2 case [19,20]. Similar considerations hold for the U < 0 case.

Summary and Outlook

We have studied a bilayer system of dipolar molecules for arbitrary orientation of the dipoles with respect to the planes. The two-body bound state structure was calculated, including energies and partial-wave decomposition of the wave function as function of dipolar strength and polarization angle. We proved that there is always a bound twobody state in the system, irrespective of strength and polarization angle of the molecules, and also verified this fact numerically. We argued that this follows from the fact that for small strength, the wave function must reside outside the region where the potential is non-zero. The results apply irrespective of the sign of the interaction strength. Negative strengths invert the dependence of energy on the dipole angle such that perpendicular polarization angle has the smallest binding energy. The structure of the wave function is dominated by the monopole component which decreases with the strength of the interaction. Up to moderate strengths, the monopole component is always larger than 50 percent while the dipole component accounts for most of the remaining probability.

The conclusion that zero net volume potentials always have at least one bound state could perhaps be reached in other ways. First, approaching this limit from small negative net volume with one bound state strongly indicates that the bound state remains. Second, a perturbation argument is tempting, e.g. let us assume that the cylindrical monopole potential always has a bound state (as shown in [START_REF] Simon | [END_REF]), and then treat dipole and quadrupole terms as perturbations. We thus extend from a Hilbert space entirely of s-waves to include also p and d-waves. To second order in perturbation theory this always gives a negative contribution and hence more binding than for the monopole potential alone. However, closer scrutiny of such argements and their practical implementation reveal that in the limit of weak potentials the perturbations are always of the same order as the initial potential. This type of perturbative argumentation always fails. The deeper lying reason is that the energy is a non-analytical function of the strength for zero net volume potentials. This is seen in the expression for the energy where zero'th and first order in the strength contribute to the non-analytical structure at the continuum threshold. Only terms higher than second order in the strength are perturbations around this strong singularity which is a characteristic feature of zero volume potentials in two spatial dimensions.

Implications for the many-body physics of a bilayer were discussed in the limit of strong-coupling where the two-body bound states are expected to be the important degrees of freedom. We conclude that the region where (quasi)-condensation of twomolecule dimers is likely to occur can be enhanced by tuning the angle of the dipoles. In particular, at the critical angle, sin 2 θ c = 1/3, where the long-wavelength part of the intralayer interaction vanishes, we expect the conditions are most favorable for accessing this phase where dimers condense. This critical angle is different from the 'magic' angle at which dipoles moving on a line become non-interacting which has been discussed in a number of previous works [8,28,29]. We also estimated the potential for a density wave instability in the bilayer [21], and demonstrated that BCS-BEC crossover and (quasi)-condensation of dimers discussed here can occur also in an intermediate coupling and low density part of the phase diagram which is below the density wave regime. The possible roton instability in the bilayer system, as discussed for the perpendicular case in [20], is currently under active investigation.

The results presented in this paper indicate that bound complexes of more than two particles must exist in the bilayer. In the case of one-dimensional tubes this was studied in some special cases in [28] and [29], while a more complete investigation of one-dimensional complexes as function of angles and dipole moment can be found in [30]. For the two-dimensional case, the method employed here can be extended to complexes of more particles and we plan such investigations in the near future. We also note that the external trapping potential that is present in each layer in experiments [23] can be easily accomodated in the current approach by introducing one-body harmonic oscillator terms.

In conclusion, we find that bound states of dimers in a bilayer consisting of one particle in each layer are generic for particles interacting through the dipole-dipole force, irrespective of the dipole strength or polarization angle of the dipoles with respect to the layers. In general, the wave function contains several partial-wave components and therefore has interesting spatial structure. This suggests that few-body states with more than two particles will also have rich structure and it also implies that the manybody physics of the system is highly non-trivial. We sketched a phase diagram for the appearance of a dimer condensate as a function of polarization angle and showed that it is enhanced around the so-called magic angle. At this point the dimer contains a large admixture of higher partial waves and we expect the collective behavior of the system to reflect this fact. The many-body problem of a bilayer with polar molecules of arbitrary polarization angle therefore deserve further investigation.

a m Φ m (s) √ s = 1 (1 + δ m0 )π 2π 0 dφ cos(mφ)Ψ(s, φ), (A.2)
where the corresponding contribution of the sin(mφ) terms is zero due to the assumed symmetry and the coefficients a m are real. The normalization of the functions Φ m (s) will be chosen later. The rather artificial separation into coefficients and functions allows a simple procedure to compute higher orders in the strength λ. The 2D radial wave functions, Φ m (s), satisfy the system of coupled radial equations

Φ m + 1 -4m 2 4s 2 Φ m + α 2 Φ m = λ l a l a m Φ l V ml (s) , (A.3)
where the matrix elements, V ml , carrying all information about the potential, are

V ml = 1 (1 + δ m0 )π 2π 0 cos(mϕ) cos(lϕ) V (s, ϕ)dϕ . (A.4)
For cylindrical potentials V ml ∝ δ ml and the different m-values decouple. The regular radial solution to equation (A.3) at the origin provides the usual boundary condition for a centrifugal barrier potential, i.e. we choose the normalization of Φ m (α, s) such that lim s→0 s -1/2-m Φ m (α, s) = 1. We inserted here explicitly the dependence on the energy parameter, α, in Φ m . We assume here the potential, and consequently also the right hand side of equation (A.3), diverges slower than 1/s 2 when s → 0. We use the Green's function formalism which in 2D is described for a cylindrical potential in [47]. For anisotropic potentials we have more generally that the mcomponents, a m Φ m , are given by

a m Φ m (α, s) = a m Φ m0 (α, s) ( A . 5 ) -λ l a l s 0 g m (α, s, s )V ml (s )Φ l (α, s ) ds . (A.6)
The last terms vanish for s = 0 and the boundary condition at s = 0 is obeyed through the first term:

Φ m0 = √ sJ m (αs)(2/α) m m! , (A.7)
where J m is the Bessel function. The Green's function in the last terms of equation (A.6) is expressed as

g m (α, s, s ) = (A.8) iπ 4 √ ss [H (1) m (αs)H (2) m (αs ) -H (1) m (αs )H (2) m (αs)] ,
in terms of Hankel functions, H

m [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. At large distance equation (A.6), with the Green's function from equation (A.8), may be rewritten as

a m Φ m (α, s) = 1 2 √ s 2 α m m! ( A . 9 ) × H (1) m (αs) l (c * ml + δ ml )a l + H (2) m (αs) l (c ml + δ ml )a l ,
where the matrix elements, c ml , are:

c ml ≡ (A.10) λ α 2 m iπ 2m! ∞ 0 √ s H (1)
m (αs )V ml (s )Φ l (α, s )ds .

Equation (A.9) contains incoming and outgoing waves and is applicable for scattering problems. For bound states, α has to be imaginary corresponding to negative energy, and H

m (αs) in equation (A.9) diverges unless the coefficient vanishes, i.e. This is an eigenvalue equation, which only has non-trivial solutions for discrete values of α, and hence for the binding energy. However, the full radial wave functions, Φ m , enter in the definitions of the matrix elements and the equations must be solved selfconsistently.

We look for solutions in the limit of very weak strength, |λ| 1, which implies that |α| (i.e. the energy) must also be very small. The m = 0 component will then dominate in the solution to equation (A.11), because the centrifugal barrier suppresses higher partial waves. This suppression becomes more pronounced with decreasing binding energy.

We now expand both coefficients, a m , and functions, Φ m , in powers of the strength, i.e. a m = λa (1) m + λ 2 a (2) m + ..., (A.12)

Φ m = Φ (0) m + λΦ (1) m + λ 2 Φ (2) m + ..., (A.13)
where we leave the coefficient a 0 (without expansion) for normalization of the total wave function. In total we then find

a (1) m a 0 = - 1 2m ∞ 0 s 1-m V m0 (s)ds , (A.14) a (2) m a 0 = - 1 2m ∞ 0 s 1 2 -m V m0 (s)Φ (1)
0 (0, s)ds (A.15)

- 1 2m i>0 a (1) i a 0 ∞ 0 s 1 2 -m V mi (s)Φ (0) i (0, s)ds , (A.16) Φ (1) 0 (0, s) = - √ s s 0 s V 00 (s ) ln(s /s)ds , (A.17) Φ (0) m (0, s) = s m √ s + s 0 (s ) 2k -(s) 2k (ss ) m s V m0 (s )ds ∞ 0 V m0 (s )(s ) 1-m ds , (A.18) Φ (1) m (0, s) = s 0 g m (0, s, s ) √ s a 0 a (2) m a (1) m 2 V m0 (s ) - a 0 a (1) m
V m0 (s )Φ This behaviour is a consequence of the attractive and repusive centrifugal barriers for m = 0 and m > 0, respectively. The energy can now be found to any order in powers of λ. To second order the results can be written as The potential dependent constants, (A i , B i ), become increasingly more complicated with powers of λ. When the net volume of the potential, λI < 0, is negative the weak binding limit for an arbitrary potential is given by the Landau expression [START_REF] Landau | Quantum Mechanics[END_REF], which also turns out to be valid for anisotropic potentials with the appropriate definition of the volume.

E = -
In the case of λI < 0, it is not necessary to retain A 2 in (A.25) to second order in λ, whereas for I = 0 it must be retained when calculating the corrections to the leading term.

When I = 0 corresponding to zero net volume of the potential the leading order term is given by A 0 , i.e.

A 0 ≡ - ∞ 0 √
sV 00 (s)Φ where only first order from equation (A.14) has to be used in the derivation. This derivation is completely general for two-dimensional, anisotropic and reflection symmetric interactions. The symmetry requirement is only a minor simplification, and omitted in the derivation in [31]. The overall results are that there is always a bound state for very weak potentials with negative or zero net volume, and the weak binding threshold behavior of the energy is given by equation (A. 25) to second order in the potential strength. The leading order term for zero net volume is given by equation (A.27) where only the first term contribute for cylindrical potentials since then V 0m ∝ δ 0m . For non-cylindrical potentials even the leading order expression in equation (A.27) is rather complicated.

Higher orders than those related to I and A 0 are found for small energy by an iteration procedure through eqs.(A.6) and (A.10). The radial solutions are computed from equation (A.6) which in turn are used to determine the C-matrix and the energy. This procedure can be repeated to give higher order corrections of both energy and wave function. Much care is necessary to include consistently all terms up to a given order because the resulting expressions contain many terms. The simplest is B 1 which is found to be The remaining expressions, A 1 , A 2 , B 2 are more complicated and we do not show them here. They can be found by expanding c 00 in equation (A.10) up to fourth order in λ, as well as c 0i up to third order in λ in eqs.(A. 19) and (A.18), and c ml up to second order for m = 0 and l = 0. Finally, we compute the determinant of the matrix c ml + δ ml , and equate to zero, i.e.

Figure 1 .

 1 Figure 1. Illustration of the setup consisting of two dipolar particles of mass M moving in parallel planes separated by a distance d. Their dipole moments, D, are assumed to be aligned by an external field at angle θ with respect to the planes.

Figure 2 .

 2 Figure 2. Illustration of the convergence of the energy with respect to the basis size: the energy E of the ground state as function of the size N basis of the basis.

Figure 3 .

 3 Figure 3. Illustration of the convergence of the radial wave function with respect to the basis size: the radial function Φ 0 (r) from equation (A.17) as function of s = x 2 + y 2 /d calculated with bases of different sizes: N short = 30 and N long = 8, 12, 16. For comparison the analytic form of the tail, K 0 (|α|r), is also shown.

Figure 4 .

 4 Figure 4. Bound state energies as a function of U for different angles. Insets show contour plots of the potentials with valleys in bright (blue) and hills in dark (red) colors.

4 θ = π 2 Figure 5 .

 425 Figure 5. Same as figure 4 but for U < 0. Note that the vertical scale is different from that of figure 4.

Figure 6 .

 6 Figure 6. The function F (U ) = -U 2 ln(|E|Md 2 /(2 2)) for different polarization angles θ. The dashed red curves are calculated numerically and the solid blue curves are from equation(7). The relatively small contribution from the A 2 -term is neglected in this comparison.

Figure 7 .

 7 Figure 7. Partial waves probabilities, P m , for the bound state wave function at U = 10 (solid black) and U = 4 (dashed red) as function of polarization angle θ for m = 0, 1, and 2. The vertical line indicates θ = θ c . Note the different vertical scales.

θ

  = π/4, U = 4

Figure 8 .

 8 Figure 8. Contour diagram for different polarization angles of the dipole-dipole potential divided by U (upper part) followed by the probability distributions below. All lengths are in units of d.

Figure 9

 9 Figure 9. Lines of vanishing chemical potential for sin 2 θ c = 1/3 (solid black) and θ = π/2 with (solid red) and without (dashed red) intralayer repulsion. Above the solid lines we expect condensation of dimers to occur (BEC), whereas in the lower right part a many-body paired BCS-like state (BCS) should be the ground-state of the system. Also shown is the line above which a density wave (DW) for θ = θ c is expected in the bilayer system (solid blue).

∞c

  l=0 ml a l = -a m . (A.11)

V

  mi (s )Φ (0) i (0, s ) ds , s 2m (ss ) m .(A.22)The equations A.17,A.19, and A.20 provide the expansions for the coefficients and the wave functions in eqs. (A.12) and equation (A.13) when s 1/|α|. The behaviour of the wave function at infinity is now given by the non-diverging piece in equation (A.6), i.e. lim s→∞ Φ m (α, s) ∼ √ sH (1) m (αs) , (A.23) or in the particular case of weak binding we have lim s→∞ Φ m (α, s) ∼ √ sH (1) m (αs)δ m0 . (A.24)

2 ( 1 +

 21 λB 1 + λ 2 B 2 ) -λI + λ 2 (A 0 + λA 1 + λ 2 A 2 ),where the leading order constant is given by

  2m V m0 (s)ds ∞ 0 √ s V 0m (s )Φ (0) m (0, s )ds ,

B 1

 1 2mV mm (s)Φ (0) m (0, s)ds .

  1 + λB 1 + λ 2 B 2 (A.29) + λ 2 (A 0 + λA 1 + λ 2 A 2 ) ln( α 2 exp(γ)) = 0 ,which then directly leads to equation (A.25) with the identifiable constants A 1 , A 2 , B 2 .If only A 1 is needed, lower orders are required at each step but then A 2 , B 2 are not obtained correctly.

‡ In SI units we have D 2 = d 2 /4π 0 when d is dipole moment of the molecules. However, in this paper we use d to denote the interlayer distance.

2µd 2 , such sampling proves very effective. However, for

Acknowledgments

We would like to thank J. R. Armstrong, G. M. Bruun, E. Demler, and C. J. M. Mathy for valuable discussions.

Appendix Appendix A. General derivation

We assume reflection symmetry in the x-axis as for the potential in equation (2) but we shall otherwise proceed with the general derivation for any potential with this symmetry. The slightly more general formulation without any symmetry can be found in [31]. We decompose the wave function into partial waves, cos(mϕ), i.e.

Appendix B. Matrix Elements

The different matrix elements with shifted correlated n-dimensional Gaussians can be calculated with the help of the following integrals,

where

For dipole-dipole interaction equation ( 2) the matrix elements of with the long-range Gaussians equation ( 21) are given as

where

is the Tricomi confluent hypergeometric function. Note that only the monopole part, V 0 (r), gives a non-zero contribution to the matrix element.