Quantum halo states in two-dimensional dipolar clusters

A halo is an intrinsically quantum object defined as a bound state of a spatial size which extends deeply into the classically forbidden region. Previously, halos have been observed in bound states of two and less frequently of three atoms. Here, we propose a realization of halo states containing as many as six atoms. We report the binding energies, pair correlation functions, spatial distributions, and sizes of few-body clusters composed by bosonic dipolar atoms in a bilayer geometry. We find two very distinct halo structures, for large interlayer separation the halo structure is roughly symmetric and we discover an unusual highly anisotropic shape of halo states close to the unbinding threshold. Our results open avenues of using ultracold gases for the experimental realization of halos composed by atoms with dipolar interactions and containing as many as six atoms.

One of the most remarkable aspects of ultracold quantum gases is their versatility, which permits to bring ideas from other areas of physics and implement them in a clean and highly controllable manner. Some of the examples of fruitful interdisciplinary borrowings include Efimov states, originally introduced in nuclear physics and observed in alkali atoms [1][2][3] , lattices created with counter-propagating laser beams 4,5 as models for crystals in condensed matter physics, and Bardeen-Cooper Schrieffer (BCS) pairing theory first introduced to explain superconductivity and later used to describe two-component Fermi gases 6,7 . In the present work, we exploit the tunability of ultracold gases to demonstrate the existence of halo states composed by a large number of atoms with dipolar interactions. Originating in nuclear physics [8][9][10] , halo dimer states have been studied and experimentally observed in ultracold gases 11 .
A halo is an intrinsically quantum object and it is defined as a bound state with a wave function that extends deeply into the classically forbidden region 12,13 . These states are characterized by two simultaneous features: a large spatial size and a binding energy which is much smaller than the typical energy of the interaction. One of the most dramatic examples of a halo system, experimentally known, is the Helium dimer ( 4 He 2 ) 14 , which is about ten times more extended than the size of a typical diatomic molecule 11 .
While most of the theoretical and experimental studies of halos have been carried in three dimensions (3D) [15][16][17][18][19][20][21][22][23][24] , there is an increasing interest in halos in two dimensions (2D) 12,25,26 . In fact, two dimensions are especially interesting as halos in 2D have different properties 12 of the 3D ones. A crucial difference between 3D and 2D geometries is that lower dimensionality dramatically enhances the possibility of forming bound states. If the integral of the interaction potential V (ρ) over all the space is finite and negative, V k=0 = V (ρ)dρ < 0 , this is always sufficient to create a two-body bound state in 2D but not necessarily in 3D. Furthermore, the energy of the bound-state is exponentially small in 2D and it can be expressed as E = − 2 /(2ma 2 ) exp[− 2 |V k=0 |/(2πm)] 27 , where a is the typical size of the bound state. An intriguing possibility arises when such integral is exactly equal to zero, V k=0 = 0 , as a priori it is not clear if a bound state exists. This situation exactly happens in a dipolar bilayer in which atoms or molecules are confined to two layers separated by a distance h and the dipolar moments d are aligned perpendicular to the plane of motion by an external field. The interaction between atoms of different layers is given by V (ρ) = d 2 (ρ 2 − 2h 2 )/(ρ 2 + h 2 ) 5/2 , where ρ is the in-plane distance. The vanishing Born integral has first lead to conclusions that the two-body bound state disappears when the distance between the layers is large 28 although later it was concluded that the bound state exists for any separation [29][30][31][32][33][34] , consistently with Ref. 35 . A peculiarity of this system is that the bound state is extremely weakly bound in h → ∞ limit. That is, a potential with depth V (ρ = 0) = −d 2 /h 3 and width h would be expected to have binding energy equal to E = − 2 /(2ma 2 ) exp(−const · r 0 /h) where r 0 = md 2 / 2 is the characteristic distance associated with the dipolar interaction and m is the particle mass. Instead, the correct binding energy, E = −4 2 /mh 2 exp(−8r 2 0 /h 2 + O(r 0 /h)) 31,32 is much smaller as it has h −2 in the exponent and not the usual h −1 . This suggests that the bilayer configuration is very promising for the creation of a two-body halo state. Moreover, the peculiarity of the bilayer problem has resulted in the controversial claim that the three-and www.nature.com/scientificreports/ four-body 36 bound states never exist in this system, and only very recently it has been predicted that actually, they are formed 37 .
In the present work, we analyze the ground-state properties of few-body bound states of dipolar bosons within a two-dimensional bilayer setup, as candidates for halo states. In particular, we study the ground-state of up to six particles occupying the A and B layers, with A and B denoting particles in different planes. In order to find the exact system properties, we rely on the diffusion Monte Carlo (DMC) method 38 with pure estimators 39 , which has been used previously to get an accurate description of quantum halo states in Helium dimers 17 , trimers and tetramers 18,19 . In addition, we report relevant structure properties of the clusters, such as the spatial density distributions and the pair distribution functions for characteristic interlayer separations.

Hamiltonian
We consider two-dimensional systems consisting from two to six dipolar bosons of mass m and dipole moment d confined to a bilayer setup. All the dipole moments are oriented perpendicularly to the layers making the system always stable. In this configuration the angular dependence of the dipolar interaction vanishes. In our model, we suppose that the confinement to each plane is so tight that there is no interlayer tunneling and excitations into the excited levels of the tight confinement are suppressed. The Hamiltonian of this system is where h is the distance between the layers. The terms in the first row of the Hamiltonian (1) are the kinetic energy of N A dipoles in the bottom layer and N B dipoles in the top layer; the first two terms in the second row correspond to the intralayer dipolar interactions of N A and N B bosons; and the last one accounts for the interlayer interactions. The in-plane distance between two bosons belonging to the same layer is denoted by ρ ij(αβ) = |ρ i(α) − ρ j(β) | , and belonging to different layers by ρ iα = |ρ i − ρ α | , where ρ i is the in-plane position. We use the characteristic dipolar length r 0 = md 2 / 2 and energy E 0 = 2 /(mr 2 0 ) as units of length and energy, respectively. We use ρ for 2D in-plane distances and r for 3D distances.
Dipoles in the same layer are repulsive, with an interaction decaying as 1/ρ 3 . However, for dipoles in different layers the interaction is attractive for small in-plane distance ρ and repulsive for larger ρ . In other terms, a dipole in the bottom layer induces attractive and repulsive zones for a dipole in the top layer. Importantly, the area of the attractive cone increases with the distance between layers h, making the formation of few-body bound states more efficient.

Structure of the bound states
We first analyze the structure of few-body clusters, composed of two, three or four particles, as a function of the interlayer separation h. To this end, we calculate the pair distribution function g σ σ ′ (r) , which is proportional to the probability of finding two particles at a relative distance r. In the case of the ABB trimers and AABB tetramers, we also determine the ground-state density distributions for different values of the interlayer separation. AB dimer. The AB dimer is strongly bound for h r 0 and its energy decays exponentially in the limit of large interlayer separation 32 . In order to understand how the cluster size changes with h/r 0 , we show in Fig. 1 the interlayer pair distributions g A B (r) (left-axis, blue curves) for three values of h/r 0 . The strong-correlation peak of g A B at h/r 0 is due to the interlayer attraction V A B (r) at short distances, also shown in the right-axis of the figure (red curves). For the cases shown in Fig. 1 we notice that g A B are very wide in comparison to the interlayer distance h/r 0 reflecting the exponential decay of the bound state. The tail at large distances becomes longer as the interlayer distance increases. ABB trimer. The ABB trimer is bound for large enough separation between the layers h/r 0 > 0.8 while, for smaller separations, it breaks into a dimer and an isolated atom 37 . The trimer binding energy is vanishingly small for h ≈ h c , with h c ≃ 0.8r 0 , and it becomes larger as h is increased, reaching its maximum absolute value at h/r 0 ≈ 1.05 . Then, it vanishes again in the limit of h → ∞ 37 . We report the intralayer and interlayer pair distributions, g B B (r) and g A B (r) respectively, in Fig. 2 for strongly-(a, b) and weakly-bound (c, d) trimers. We observe that the g A B distributions are very wide in comparison to h, similarly to what has been observed in Fig. 1 for dimers. The same-layer distribution g B B vanishes when r/r 0 → 0 as a consequence of the strongly repulsive dipolar intralayer potential at short distances. As r increases, g B B exhibits a maximum, next it monotonically decreases with r/r 0 . For a weakly-bound trimer ( h/r 0 = 1.6 ), both g A B and g B B produce long tails at large distances.
The trimer is weakly bound close to the threshold, h → h c , and for large interlayer separation, h → ∞ , but its internal structure in those two limits is significantly different. This can be seen in Fig. 3a, b, where we plot the trimer ground-state spatial distribution for h/r 0 = 1.05 and 1.6. The spatial distribution is shown as a function of the distance between two dipoles in the same layer |r B 1 − r B 2 | (horizontal axis) and the minimal distance between dipoles in different layers min{|r A 1 − r B 1 |, |r A 1 − r B 2 |} (vertical axis). For large separation between layers, shown in Fig. 3b for h/r 0 = 1.6 , the distances between AB and BB atoms are all of the same order, revealing an approximately symmetric structure. However, by decreasing the distance between layers the particle distribution becomes significantly asymmetric. For h/r 0 = 1.05 (Fig. 3a), we observe that the trimer spatial distribution is elongated: two dipoles in different layers are close to each other while the third one is far away. Regardless of the   www.nature.com/scientificreports/ the trimer becomes more extended and breaks into a dimer and a single atom at h/r 0 ≈ 0.8.

AABB tetramer.
As we have shown in the previous Section, an ABB trimer dissolves into a dimer and an atom when h ≈ h c . Here, we address the structure properties of the balanced case for a tetramer, in which the number of A and B atoms is the same. The AABB tetramer is weakly bound for large values of h/r 0 . When the distance between layers decreases, the tetramer becomes unbound at h/r 0 ≈ 1.1 and splits into two AB dimers 37 . The pair distributions g A B and g B B for AABB are shown in Fig. 2e-h for two characteristic values of the interlayer distance h/r 0 . We observe a behavior that is similar to that previously reported for ABB. That is, both g A B and g B B are compact for the deepest bound state ( h/r 0 = 1.3 for AABB) and become diffuse, showing long tails at large distances, when it turns to a weakly-bound state ( h/r 0 = 1.6).
The ground-state spatial distributions for the symmetric tetramer are shown in Fig. 3c, d. We observe that for large separation h, i.e., when the tetramer is weakly bound, it has large spatial extension and the distances between AA and AB pairs are of the same order. As the interlayer separation is progressively decreased, the tetramer size decreases and its structure becomes anisotropic. In this case, the distance between dipoles in the same layer is several times larger than the distance between dipoles in different layers. When the tetramer approaches the threshold for unbinding the cluster becomes even more elongated and it breaks into two AB dimers at h/r 0 ≈ 1.1.

Quantum halo characteristics
A halo is a quantum bound state in which particles have a high probability to be found in the classically forbidden region, outside the range of the interaction potential. The key characteristics of a halo are its extended size and binding energies much smaller than the typical energy of the interaction. In order to classify a system as a halo state, one typically introduces two scaling parameters with which the size and the energy are compared. The first parameter is the scaling length R. For two-body systems one commonly chooses R as the outer classical turning point. The second parameter is the scaling energy µBR 2 / 2 , where µ is the reduced mass and B is the absolute value of the ground-state energy of the cluster. The size of a cluster is usually quantified through its mean-square radius ρ 2 , where ρ is the interparticle distance. A two-body quantum halo is then defined by the condition ρ 2 /R 2 > 2 , which means that the system has a probability to be in the classically forbidden region larger than 50% 12 .
The dipolar interaction in the bilayer geometry has vanishing Born integral and thus, the AB dimer can show an enhancement of its halo properties. In Fig. 4e, we show the scaling plot for the dipolar dimers, corresponding to interlayer distance from h/r 0 = 0.14 to 1.6, as indicated on the upper axis. All dimers which lie above the halo limit �ρ 2 �/R 2 = 2 (horizontal line in Fig. 4e) are halo states and follow a universal scaling law �ρ 2 �/R 2 = 2 /(3µBR 2 ) , shown with a dashed line in the figure 12 . This is exactly the case for all dimers with interlayer separations h/r 0 > 0. 45 . This threshold value is close to the characteristic value, h/r 0 = 0.5 , for which the dimer binding energy is approximately equal to the typical energy of the dipolar interaction E A B ≈ 2 /(mr 2 0 ). While AB dimers exist for any interlayer separation, ABB trimers and AABB tetramers are self-bound for large h values, where AB dimers are in fact halo states. Thus, it can be anticipated that these few-body bound states are also halos. To verify that, the sizes of three-and four-body systems should be compared to the meansquare hyperradius 12 , www.nature.com/scientificreports/ where m * is an arbitrary mass unit, M is the total mass of the system, and m i is the mass of particle i. In analogy to the hyperradius Eq. (2), Jensen et al. 12 defined the scaling size parameter ρ 0 as with R ik the two-body scaling length of the i − k system, which is calculated as the outer classical turning point for the i − k potential. Considering the problem of a dipolar dimer in the bilayer, the two-body scaling length R E,2 is related to the interlayer interaction potential V AB (ρ, h) and the dimer energy E 2 according to For a given value of the interlayer distance h the classical turning points are given by the roots of Eq. (4) which is a polynomial of tenth degree in R E,2 . The exact value of the dimer energy E 2 depends on h and can be obtained by solving numerically the Schrödinger equation with the potential V AB (ρ, h) or using the DMC method. In Table 1, we report the classical turning points R E,2 corresponding to interlayer distance from h/r 0 = 0.4 to 1.6. On the other hand, for repulsive potentials we choose R ik equal to zero. The condition for three-and four-body quantum halos is then ρ 2 HR /ρ 2 0 > 2 12 . While there is a general agreement on the definition of the classical turning point R E,2 in the two-body problem, Eq. (4), in the few-particle case there is not a single way to determine the scaling length R ik . It has been suggested 24 to consider R ik as the distance at which the attractive interaction Eq. (4) equals a characteristic energy of the system. One possibility is to obtain the classical turning point from the relation, where the brackets �· · · � denote averaging of the two-body interaction potential over the many-body bound state. Another possibility 24 , denoted by R E,N , is based on comparing the two-body potential to the bound-state energy B N divided by the number of pairs N(N − 1)/2,  Table 2, we report the values of R E,2 , R E,N , and R pot,N for h/r 0 = 1.2 and for N = 2, . . . , 6 . One notices that in the bilayer system the three different definitions, Eqs. (4-6), result in very similar values. Therefore, choosing one of these radii over the other two to calculate ρ 0 is a minor effect. In the following we will use R E,2 to calculate the scaling size parameter ρ 0 . The dependence of the scaled size on the scaled energy for ABB and AABB are shown in Fig. 4f. We find a non-monotonic behavior, in clear contrast with the dependence observed in the dimer case (see Fig. 4e). That is, the cluster size decreases with increasing energy and reaches a minimum and then it starts to grow again. The minima correspond to the deepest bound states 37 . This resurgence appears as the clusters approach the thresholds, where trimers eventually break into a dimer and an atom, and tetramers into two dimers. We want to emphasize that all the trimers and tetramers analyzed in Fig. 4 are halo states, although they are organized in significantly different spatial structures. On the left side of the minima, the clusters are almost radially symmetric and all the interparticle distances are of the same order. However, at the minima and on the right side of the minima the cluster structures are elongated and highly asymmetric. The AABBB pentamer and AAABBB hexamer are selfbound and are manifestly halo states. Their mean square size has a similar behavior to the one observed before for the trimer and tetramer, that is a minimum corresponding to the larger binding energy which separates a regime of nearly symmetric particle distribution from another one, more elongated, and thus asymmetric. It is important to notice that the existence of halo tetramers, pentamers and hexamers is very unusual as it contradicts the usual tendency of self-bound clusters to shrink and lose the halo character as the number of particles is increased. As we demonstrate in the present work, the bilayer geometry is very promising for creation of halo states with up to six particles.
Experimental signatures of halo and Efimov states come from the measurement of three-and four-body recombination loss rate in ultracold gases 2,42,43 or by using techniques such as matter-wave diffraction combined with laser-induced Coulomb explosion imaging 14,44 . Using the latter techniques the halo states of helium dimer 14  A possible issue of concern may be the validity of our findings in an experimental realization, where the bilayer has a quasi-two-dimensional geometry and not strictly two-dimensional one as in our model. To this aim, we have compared the dimer energy of our model with the dimer energy of a quasi 2D model. For typical experimental parameters, we have found that the change in the dimer energy is at most 20% . Therefore, it can be concluded that the effects of considering a quasi-two dimensional model do not change our main findings.

Conclusions
We used the diffusion Monte Carlo method to study the ground-state properties of few-body bound states of dipolar bosons in a two-dimensional bilayer setup. We have studied clusters composed by up to six particles, for different values of the interlayer distance, as candidates for quantum halo states.
In the case of dimers, we find that for values of the interlayer separation larger than h/r 0 = 0.45 the clusters are halo states and they follow a universal scaling law. In the cases of trimers up to hexamers, we find two very different halo structures. For large values of the interlayer separation the halo structures are almost radially symmetric and the distances between dipoles are all of the same scale. In contrast, in the vicinity of the threshold for unbinding, the clusters are elongated and highly anisotropic. To our knowledge halo states have been experimentally observed for up to four particles 2,14,42-44 . The addition of particles to a two or three body halo www.nature.com/scientificreports/ states typically makes them shrink towards a more compact liquid structure. Importantly, our results prove the existence of stable halo states composed by atoms with dipolar interactions and containing up to six particles. As commented before, there are reasons to believe that our system is experimentally viable with current technology as compared to other theoretical predictions, where zero-and finite-range interactions are used 23,24 . We conclude that the bilayer geometry is advantageous for the observation of halo states in future experiments. We hope that these results will stimulate experimental activity in this setup, composed by atoms with dominant dipolar interaction, to bring evidence of these quantum halo states. In outlook, our results can stimulate further theoretical and experimental research of halo states in ultracold gases. It could be interesting to understand how the possibility of tunneling between the layers affects the stability of the halo states. An interesting new path of research could be to study halo states in a bilayer system of fermionic dipoles 51,52 .

Method
To investigate the structural properties of the dipolar clusters, we use a second-order DMC method 38 with pure estimators 39 . This method allows for an exact estimation of the ground-state energy, as well as other properties, within controllable statistical errors. The DMC stochastically solves the imaginary-time Schrödinger equation, where τ = it/ is a imaginary time, E s is an energy shift and the walker ρ = (ρ 1 , . . . , ρ N ) is a vector containing positions of N particles. Importance sampling is used to reduce the statistical noise of the calculation, which consists on rewritten the Schrödinger equation, Eq. (7), for the mixed distribution �(ρ, τ ) = �(ρ, τ )ψ(ρ) . We use a trial wave function ψ of the form which is suitable for describing systems with short-range correlations and as well as long-range asymptotic behavior.
The trial wave function for intraspecies correlations is built from the zero-energy two-body scattering solution K 0 (ρ) being the modified Bessel function and C 0 a constant. The interspecies interactions are described by the dimer wave function f A B (ρ) up to R 0 . From the variational distance R 0 on we took the free scattering solution f A B (ρ) = CK 0 ( √ −mE A B ρ/ ) . We impose continuity of the logarithmic derivative at the matching point R 0 , yielding the following equality In a DMC calculation, the expectation value of an observable Ô is obtained for long enough imaginary time propagation with 0 the ground-state wave function. The last equation is known as the mixed estimator. Equation (11) gives the exact expectation value for the Hamiltonian and for observables that commute with it. In the case of operators that do not commute with Ĥ , the result obtained from Eq. (11) can be biased by ψ . In this case, it is possible to obtain exact expectation values using the pure estimators technique 39 . In the present study, pure estimators are used for the calculation of the pair distribution functions, the spatial distribution functions, and the size of the clusters.