Quantum magnetism and counterflow supersolidity of up-down bosonic dipoles

We study a gas of dipolar Bosons confined in a two-dimensional optical lattice. Dipoles are considered to point freely in both up and down directions perpendicular to the lattice plane. This results in a nearest neighbor repulsive (attractive) interaction for aligned (anti-aligned) dipoles. We find regions of parameters where the ground state of the system exhibits insulating phases with ferromagnetic or anti-ferromagnetic ordering, as well as with rational values of the average magnetization. Evidence for the existence of a novel counterflow supersolid quantum phase is also presented.


I. INTRODUCTION
Recently, the physics of ultracold dipolar gases has received growing attention. Experimental and theoretical studies are motivated by the long-range and anisotropic dipole-dipole interaction which introduce a rich variety of quantum phases [1]. For instance, supersolid and checkerboard phases are predicted in the phase diagram of polarized dipolar Bosons confined in a two-dimensional (2D) optical lattice [2]. In bilayer samples where dipole interactions can also be attractive more exotic quantum phases are accessible, for example a pair-supersolid [3].
Experiments with ultracold atoms have now highlighted the role of dipole-dipole interactions. Notably, with 52 Cr, control via Feshbach resonances allows one to efficiently reduce contact interactions and enter a regime where magnetic dipole-dipole interactions become dominant [4]. Polar molecules also appear very promising towards implementation of degenerate dipolar gases (see [1] and References therein), in particular since the demonstration of ultra-cold rubidium-potassium, and lithiumcesium molecules prepared in their ground rotovibrational state [5,6].
Note that in most cases, if not all cases considered so far, dipolar gases are polarized, i.e. magnetic or electric dipoles point in the same direction. This, however, does not have to be the case always. A prominent example is that of exciton gases. Excitons, which are bound electron-hole pairs, obviously carry an electric dipole moment that can a priori attain quite arbitrary directions. Nevertheless, in indirect quantum wells, electrons and holes are confined in spatially separated regions such that electric dipoles can be aligned [7] or even anti-aligned in type-II heterostructures [8].
Another extreme case concerns molecules that follow the Hund's rule (a) [9]. For such molecules the electric dipole is either parallel or anti-parallel to the direction of the magnetic moment. If one takes thus a sample of such molecules and polarizes it magnetically, say in the up direction, one will obtain in this way a dipolar gas with certain fraction of electric dipoles pointing up, and the remaining fraction pointing down. This is the situation that we study in the present paper. In fact there was recently a spectacular progress in cooling and trapping of magnetically confined neutral OH molecules. This progress opens important perspectives to study novel quantum phases with ultracold dipoles [10,11]. In its ro-vibronic ground state OH is indeed a Hund's case (a) molecule, and it features both an electric d = 1.67D and a magnetic dipole moment µ = 1.2µ B (with µ B being the Bohr magneton), which can be aligned or anti-aligned independently. Using these two degrees of freedom, OH molecules may be confined in a 2D optical lattice, where the dipolar interactions lead to rich quantum phases, as we show in this work.
We study here a sample of bosonic dipoles confined in a 2D lattice. We consider that dipoles are free to point in both directions normal to the lattice plane. This constitutes a novel ingredient compared to our previous works [12,13], and results in a nearest neighbor interaction either repulsive for aligned dipoles, or attractive for anti-aligned ones. We consider the case of dipolar interactions to be relatively weak compared to the contact interactions, and remarkably we find that the dipole-dipole interactions dominate the physics of the system. Using a mean-field Gutzwiller approach, we show that the system presents Mott-insulating phases with ferromagnetic or anti-ferromagnetic ordering, as well as with fractional values of the average magnetization, depending on the imbalance between the population of dipoles pointing in opposite directions. We also find evidences of a novel counterflow supersolid quantum phase. The latter exhibits broken translational symmetry, namely a modulation of the order parameter on a length scale larger than the one of the lattice spacing, analogously to supersolid phases.

II. HAMILTONIAN OF THE SYSTEM
We consider a sample of dipoles in the presence of a 2D optical lattice, and an extra confinement in the perpendicular direction. The dipoles are free to point in both directions normal to the lattice plane, as shown in Fig. 1. The particles feel repulsive intra-species Uaa, U bb , and inter-species U ab repulsive on-site energies. The first nearest-neighbor interaction is repulsive UNN > 0 for aligned dipoles, while it is attractive −UNN for anti-aligned particles while the hopping term J is equal for both the species.

The system is described by the Hamiltonian
where σ = a, b indicates the type of species, i.e. dipoles pointing in the up and down direction perpendicular to the 2D plane of the lattice, respectively. U aa and U bb are the on-site energies for particles of the same species, while U ab is the on-site energy for different species. The long-range dipolar interaction potential decays as the inverse cubic power of the relative distance r ij , which we express it in units of the lattice spacing. For computational reasons, in most theoretical approaches the range is cutoff at a certain neighbors. In the present work, we consider a range of interaction up to the 4th nearest neighbor. The first nearest neighbor dipolar interaction is attractive (repulsive) for particles of the same (different) species, with strength U N N > 0 (or respectively −U N N ), and we consider an equal tunneling coefficient (J) for both the species, while the densities of the dipoles pointing upwards and downwards are fixed by the corresponding chemical potentials µ σ .
The on-site interactions have two contributions: one is arising from the s-wave scattering length U s = 4πh 2 as m d 3 rρ(r) 2 , and the second one is due to the on-site dipole-dipole interaction where U dd (k) and ρ 2 (k) are the Fourier transform of the dipole potential and the density respectively. We consider that the swave scattering length is independent on the orientation of the dipoles. Instead, the on-site dipolar contribution U dd depends both on the orientation of the dipoles and on the geometry of the trapping potential, and it can be varied by changing the ratio between the vertical to the axial confinement. For simplicity we will focus on the specific case of a spherically symmetric confinement, where the on-site dipolar interactions average out to zero U dd = 0, and the resulting on-site interactions are all equal to U . We consider the case of dipole-dipole interactions to be 600 times weaker with respect to the on-site interaction, i.e. U NN = U/600 [14].

A. Filling factor and imbalance
The properties of the system are conveniently extracted using the operators given by the sum (filling factor) and by the difference (imbalance) of the two species number operators at each site of the lattice, namely bŷ which are simultaneously diagonal on a given Fock state |ν, m i . Notice that the eigenvalues of these two operators are not independent. In fact, by fixing ν the eigenvalues ofm i can only assume 2ν + 1 values given by m = {−ν, −ν + 1, ..., +ν}, in complete analogy with the angular momentum operatorŜ 2 i and its projection along the z axisŜ z i , as we will discuss in Sec. II B. It is also useful to introduce the average magnetization of the system, defined as where N S is the total number of lattice sites. Substituting Eqs. (2) into equation (1) allows us to express the system Hamiltonian asĤ =Ĥ ν 0 +Ĥ m 0 +Ĥ νm 1 , whereĤ H νm In Eqs. (4,5), we have introduced the chemical potentials which respectively fix the eigenvalues of the filling factor (+), and the imbalance operators (−) in Eq. (2). In the following we considerĤ νm 1 to be a small perturbation on the interaction terms (4) and (5). In the limit where U ≫ U N N ,J , the ground state of the system is found to be a uniform distribution of constant filling factor ν i =ν at each site of the lattice. The value ofν is fixed by µ + , and can be integer as well as semi-integer. This is better understood at J = 0, where we can calculate the expectation value ofĤ ν 0 on a given classical distribution of atoms in the lattice |Φ = i |ν i , m i i , as follows where ν i = Φ|ν i |Φ . In the right hand-side of Eq. (8) all sites i are equal, and like in the homogeneous case of a Bose-Hubbard Hamiltonian at J = 0 [15], the minimum of Eq. (8) is provided by a uniform distribution in the lattice, where ν i =ν at each site. Instead, for a givenν, finding the distribution of m i = Φ|m i |Φ , which minimize the expectation value is non-trivial due to the presence of the long-range interaction. However, we can qualitatively argue that for |µ − | ≫ U N N , the minimum of the energy (9) is obtained for m i = ν × sign µ − ∀i, which corresponds to a ferromagnetic phase (FM) of average magnetization M = ν × sign µ − , where only particles of one species are present. Instead, for µ − = 0, a succession of nearest neighbors with m i = ν and m j = −ν provides the minimum of Eq. (9), and the phase is anti-ferromagnetic (AM), i.e. M = 0. In other words, the spatial distribution of the particles is given by sites occupied from the species a alternated with sites occupied by the species b in a checkerboard-like structure. In Fig. 2 we plot the ground state phase diagram at J = 0, in the µ − vs. µ + plane, where the text in parentheses indicates filling factor and average magnetization, i.e. ν and M respectively. Note that between the AM, and the FM phases, we find magnetic phases which present non trivial rational values of the average magnetization (RM). The precise choice of the cutoff range in the dipolar interaction, and the lattice size, determine the fractional character of the allowed ground state magnetization.
In the next section, we include the presence of tunneling. The physical properties of the system are described by an effective Hamiltonian, supported by the existence of a low-energy subspace. Hence, the theoretical description of the system cannot be based on standard mean-field theory, which is not suitable to correctly describe the ground state of the system, as we discuss in the following section.

B. Low-energy subspace and effective Hamiltonian
The ground state of the system at J = 0 is described by a product over single-site Fock states of the type with uniform total on-site occupation 2ν. As single particle hopping changes the total on-site population, it breaks the translational invariance of the ground state with respect to the total on-site occupation 2ν. The energy cost of these excitations is of the order of the on-site interaction energy U , and is therefore very costly in the limit where U ≫ U N N , J . On the contrary, exchanging two particles from nearest neighboring sites does not require such a large amount of energy. This defines a low-energy subspace spanned by the |α configurations at constant filling factor ν of Eq. (10), which is energetically well separated from the rest of the Hilbert space in the limit of parameters we consider. Thus, a successful description of such a system is obtained through an effective HamiltonianĤ eff restricted to the low-energy subspace, where single-particle hopping is suppressed and tunneling is included at second order in perturbation theory. The validity of the effective Hamiltonian relies on the existence of this low-energy subspace well separated in energy from the subspace of virtual excitations, to which it is coupled via single-particle hopping. The relevant virtual subspace is then obtained from the states |α via single particle hopping, and it is spanned by the states as schematically represented in Fig. 3.
FIG. 3: Schematic representation of a two-particle hopping between the states |α and |β , belonging to the low-energy subspace at ν = 1/2. These states are coupled through virtual excitations to the states |γ by single-particle hopping.
This situation is in fact similar to the one discussed in [3] for a bilayer optical lattice, and therefore we apply the same technique to computeĤ eff . In the basis of constant total on-site population 2ν, the matrix elements of such a Hamiltonian at second order in perturbation theory are given by whereĤ 0 =Ĥ ν 0 +Ĥ m 0 , given by the sum of the interaction terms (4) and (5), is diagonal on the states |α , and the single-particle tunneling termĤ νm 1 of Eq. (6) is treated at second order. For a given state |α , whereĉ i =â ib † i andĉ † i =â † ib i are composite operators, corresponding to the creation of a particle of one species and a hole of the other species, such that while their commutation relation obeys For a given ν, the second line of the Hamiltonian (14) can be equivalently written in terms of the spin operators at site i [17,18], given bŷ where σ = (σ x , σ y , σ z ) are the Pauli matrices, and u = a, b indicates each species. Thus, the creation and annihilation operators (15) becomeĉ † i =Ŝ x i + iŜ y i and c i =Ŝ x i − iŜ y i respectively, while the imbalance operator is given bym i =Ŝ z i as already anticipated in Sec. II A. Therefore, in a spin representation as in Eq. (16), the second line of the Hamiltonian (14) acts as a Heisenberg spin Hamiltonian (see e.g. [16]). The chemical potential µ − plays the role of an external magnetic field along the z axis. The interplay between µ − , and the long-range interaction determines the magnetic ordering of the system, as we will discuss in the next section.

III. MEAN-FIELD
In this section, we provide a mean-field solution to the effective Hamiltonian (14) to investigate quantum phases of the system. We identify the different phases through the composite order parameters ĉ i , as well as both the single-particle ones â i , and b i .
For every subspace at constant filling, we find that the system presents three different kinds of phases. The Mott-insulating phase (MI), with a well defined number of particles at each site of the lattice, and absence of any low-energy transport [15]. The MI is characterized by vanishing ĉ i = â i = b i = 0, and depending on the value of µ − presents FM, AM or RM ordering. The second quantum phase, a super-counter-fluid phase (SCF), exhibits on-site density fluctuations where a net transport of atoms is suppressed but a counterflow is present, for which the currents of the two species have equal amplitudes but opposite directions [17]. In the SCF phase, while the single-particle order parameters still vanish â i = b i = 0, the composite order parameters are nonzero ĉ i = 0, indicating the presence of counterflow. We also find evidence of a third and novel quantum phase, namely a counterflow supersolid phase (CSS). The CSS is characterized by vanishing single-particle order parameters â i = b i = 0, and non-vanishing composite order parameters ĉ i = 0, coexisting with broken translational symmetry, namely, a modulation of both m i , and ĉ i on a scale larger than the one of the lattice spacing, analogously to the supersolid phase.
To determine the insulating phases we perform a perturbative treatment at first order in the composite order parameters ψ i = ĉ i , which allows us to compute the boundaries of the insulating lobes. Furthermore, we solve the time dependent Gutzwiller equations in imaginary time to determine the nature of the SCF-CSS phases outside the lobes.

A. Insulating lobes
The low-energy subspace is spanned by the classical distribution of atoms in the lattice |α of Eq. (10). Similarly to the two layer system discussed in [3], in the limit where U ≫ U NN , asymptotically all classical states |α become stable with respect to single-particle-hole excitations and develop an insulating lobe at finite J. The energy of single particle-hole excitations is of the order of U at J = 0 and is given by the width of the lobes at finite J (see, e.g., the thin blue/black lobes in Fig. 4).
Instead, the low-lying excitations remain within the subspace and are obtained by adding (PH) or removing (HP) one composite, made of a particle of the upperpolarized dipoles (species a), and a hole of the lowerpolarized dipoles (species b), at the i-th site of the lattice. This corresponds to flipping the direction of a dipole at the site i, respectively from down to up (PH) or from up to down (HP). For any given configuration |α , one can calculate the corresponding energy costs using the diagonal terms of the effective Hamiltonian (14), which are respectively given by Note that in the last expressions, there is no explicit dependence on the chemical potential µ + . This is because by adding or removing one composite, we remain within the subspace at filling factor ν, and therefore the contribution of µ + vanishes in the calculation of the excitations (17). By using a perturbative mean-field method, we can calculate the order parameters ψ i = ĉ i for |α , which satisfy the equations whereψ i = j i ψ j . With Eqs. (18) one can calculate the mean-field lobes of any distribution of atoms in the lattice |α , provided that the elementary excitations (17), are positive in some parameters range. For every site i of the lattice, one can adequately flip the direction of a dipole, and depending on m i either one of the following conditions apply For example, suppose the site i is occupied only by one particle of the species a, i.e m i = 1/2, then for this site the conditions (19) since the only possible excitation at this site corresponds to remove a composite. Conditions (19) are necessary for the existence of an insulating lobe and provide its boundaries at J = 0. For each site of the lattice one has such a condition (19), and an Eq. (18). The latter constitutes a set of coupled equations, which can be written in the matrix form M(µ − , U, U NN , J) · ψ = 0, with ψ ≡ (· · · ψ i · · · ) T being the vector of the order parameters at each site of the lattice. For every µ − , a nontrivial solution is provided by the smallest J for which det[M(µ − , U, U NN , J)] = 0, that is the insulating lobe of the |α configuration in the J vs. µ − plane. In Fig. 5 we plot the ground state insulating lobes calculated in this way for ν = 1/2 (left) and ν = 1 (right). For all filling factors ν, we find an AM ground state (ν, M = 0), which presents a spatial distribution of alternating sites occupied by particles of species a and b resembling a checkerboard structure. Remarkably, the larger the ν, and the more stable is the AM ordering with respect to flipping the direction of a dipole. Indeed, using the inequalities (19), it is not difficult to calculate the boundaries at J = 0 of such a checkerboard structure, which, for a fixed ν are given by −2zνℓ . Instead the tip of the AM lobes is found to be ν-independent, and given by J/U = ℓ U NN /2U at µ − = 0. By increasing the absolute value of µ − we find RM ground states with rational values of the average magnetization, corresponding to M = (±2ν, ±4ν, ±6ν)/8. The exact fractional values of M in the ground state, depend on the cutoff range in the dipolar interactions, and on the size of the lattice. We have used a 4 × 4 elementary cell with periodic boundary conditions, and dipolar interaction range cut at the 4th nearest neighbor. By considering more neighbors in the interactions, and larger lattices, we expect to find RM states appearing at all rational M , asymptotically approaching a Devil's staircase as recently shown in [21,22]. Finally we find a FM ground state (ν, M = ±ν), in which only particles of one type are present.
It is worth noticing that the insulating lobes calculated in this way, do not contain any dependence on µ + , which does not enter into Eqs. (18) as previously discussed. Therefore, for any given value of µ + , in order to obtain the ground state phase diagram one has to compare the energies of the ground state configurations at different ν. Using the effective Hamiltonian (14), for any value of µ + , J, and µ − , we compare the energies of the ground state configurations for different ν, and select the state with the smaller energy. In this way we have obtained the phase diagram at J = 0 of Fig. 2.

B. Counterflow superfluid/supersolid
In the low-energy subspace at constant ν, the Gutzwiller Ansatz on the wave function of the system

reads
where we allow the Gutzwiller amplitudes to depend on time, i.e. f (i) ν,m (t). We obtain the equations of motion for the amplitudes by minimizing the action of the system, given by S = dtL, with respect to the variational parameters f is the Lagrangian of the system in the quantum state |Φ [23]. Therefore setting to zero the variation of the action with respect to f * (i) ν,m leads to the equations j =i m j /|r ij | 3 have to be calculated in a self consistent way, and the order parameter ψ i = Φ|ĉ i |Φ is given by We solve Eqs. (22) in imaginary time τ = it, which due to dissipation is supposed to converge to the ground state. In Figure 5 we show the ground state phase diagram of the system for ν = 1/2 (left) and ν = 1 (right), computed in this way for U NN = U/600. The phase diagram is symmetric with respect to the µ − = 0 axis. This is because the description of the system is identical under the interchange of the species a with the species b, and vice versa. For ν = 1/2, in the region immediately outside the insulating AM lobe, depending on the values of J and µ − we find either SCF or CSS. Note that we don't find any evidence of CSS at µ − = 0, which indicates that a finite imbalance between the two species is a necessary condition for the system in order to sustain CSS. Let us underline that for ν = 1/2, our effective Hamiltonian (14) can be mapped onto a hard-core dipolar Bose-Hubbard Hamiltonian, provided that we neglect a constant factor proportional to ν, and rescale the first nearest neighbor interaction by the small quantity −2J 2 /U . Contrary to mean-field predictions, Monte Carlo studies of the hard-core Bose-Hubbard Hamiltonian on square lattices, with interactions extended to first and second nearest neighbors, have shown that no supersolid phase is obtained by doping the checkerboard solid (AM ordering) [19,20]. However, it was recently demonstrated by Monte Carlo simulations that considering an infinite range in the dipolar interactions stabilizes the supersolid phase, obtained by doping the checkerboard solid (AM) either with particles or holes [21], in agreement with mean-field predictions. In our treatment we are approaching this limit, and we therefore expect the CSS phase to be stabilized by the long-range interactions. Furthermore, because J > 0, the rescaling of the first nearest neighbor interaction by the small quantity −2J 2 /U in the mapping process, should enforce further this effect. Although our treatment is approaching infinity in the interaction range, it would be important to verify our predictions with first-principle quantum Monte Carlo simulations.
To get reliable results, one should combine the Gutzwiller predictions with an estimate of the limits of validity of H (0) eff , beyond which the subspace of constant ν looses its meaning. Before starting the discussion on the validity of the subspace, let us explain how we define the dominant classical configurations of a given state |Φ . It is not difficult to see that Eq. (20) can be equivalently written as where m = (m 1 , ..., m i , ..., m N S ) is a collection of the indices m at each site, and we have introduced the notation ν,mi . The advantage of writing the Gutzwiller state |Φ in the form (24), lays in the product over singlesite Fock states |α = i |ν, m i i , which is nothing but a classical distribution of atoms in the lattice. Therefore we can rewrite Eq. (24) as For each point of the phase diagram, from the ground state Gutzwiller wavefunction, we define the dominant classical configurations with the criteria |g m | = | i f (i) ν,mi | > (0.02) 4 , and we require |f (i) ν,mi | 2 > 0.05, implying that each of the contributing f (i) ν,mi should also be sufficiently large. For each of these configurations, we calculate the lobe with respect to single particle-hole excitations [27]. If the system at this given point of the phase diagram turns out to be stable against all dominant single particle-hole excitations (in other words, if this point is inside all selected single particle-hole lobes), H (0) eff is considered valid. This procedure is shown for µ + = 0.4U , J = 0.035U and µ − = 0 in Fig. 4, and gives the thick black lines of Fig. (5). On the right hand side of these black lines, the low energy subspace is not well defined and therefore the effective Hamiltonian looses its meaning, leaving the description of the system to the domain of single-particle single-hole excitation theory that predicts SF and SS phases for each component separately.
We have already mentioned, that the boundaries of the lobes calculated with the effective Hamiltonian do not show any dependence on the chemical potential µ + , which does not give any contribution in the expression of the  (17). This is not true in the case of the single-particle-hole insulating lobes, since adding or removing a single particle results in a change of both µ + , and µ − . This makes the process of estimating the limits of validity ofĤ (0) eff more complicated and leads to a 3D phase diagram in J, µ + , and µ − highly non-trivial. In Fig. 6 we present slices of the 3D phase diagram, calculated at constant values of µ + , for ν = 1/2 only. From Fig. 6, it is evident that the insulating magnetic phases persist for a wide range of µ + values. Instead, the SCF and the CSS phases are more affected by the limits of validity ofĤ (0) eff , which may imply high degrees of control in order to experimentally observe this quantum phases.

low-lying excitations
Finally, it is useful to estimate the critical temperatures at which we expect these quantum phases to be experimentally observable. Bose-Einstein condensation in alkali atoms occurs at temperatures of the order ∼ 100 nK. Assuming that condensation of OH molecules occurs at similar temperatures we expect the SCF to be observed at ∼ 100 nK. For the AM insulating phases, the gap at J = 0 is of the order ∼ 2ν × 10 −2 U , and the critical temperature is expected to be of the order of the standard insulator-superfluid transition temperature, i.e. 10 to 50 nK [24]. Therefore, the CSS phase is expected to be observed somewhere in between 50 to 100 nK. The lowest measured temperature of ultracold atoms in optical lattices is ∼ 1 nK [25], where superexchange interactions, of the order ∼ J 2 /U , are dominant. Recently, a two-component hard-core Bose-Hubbard Hamiltonian was investigated with Monte Carlo techniques [26], and insulating AM as well as SCF phases where observed. In [26], the largest critical temperature observed for the AM phase is T c ∼ 0.12J, at J/U = 0.0125, and T c ∼ 0.104J at U/J = 13 for the SCF phase.

IV. CONCLUSIONS
In conclusion, in this work we have studied a gas of dipolar Bosons confined in a two-dimensional optical lattice, where the dipoles are free to point in both directions perpendicular to the lattice plane. We have found regions of parameters where the ground state presents insulating phases with antiferromagnetic or ferromagnetic ordering, as well as with rational values of the average magnetization. Our mean-field calculations predict the existence of a novel counterflow super solid quantum phase, which presents a spatial modulation of the density coexisting with the presence of counterflow. Our work in a sense is the first step in the studies of dipolar gases with non-polarized dipoles, and is relevant for experimental studies of the ultracold OH molecules. We expect that the methods and ideas developed here will turn out to be useful also for the study of more complex systems, such as excitons in indirect quantum well structures, or completely unpolarized dipolar gases.