Many body population trapping in ultracold dipolar gases

First we consider the simple non-trivial situation capturing the essential physics: two neighbouring sites i and j share a single pair localized in either of the sites. Due to the action of ${{H}_{WB}}$, states in the neighbouring sites will be coupled via transitions to the common continuum. The state of a pair distributed among sites i and j may be written as,

$$\begin{eqnarray} \left| \Phi \right\rangle & = & \underset{l=1}{\overset{3}{\mathop{\mathop{\sum }^{}}}}\,{{C}_{l}}{{\left| l \right\rangle}_{i}}{{\left| 0 \right\rangle}_{j}}\left| \mathbf{0} \right\rangle+\underset{l=4}{\overset{6}{\mathop{\mathop{\sum }^{}}}}\,{{C}_{l}}{{\left| 0 \right\rangle}_{i}}{{\left| l-3 \right\rangle}_{j}}\left| \mathbf{0} \right\rangle \\ {} & + & \underset{{{n}_{1}}{{n}_{2}};{{q}_{1}}{{q}_{2}}}{\mathop{\mathop{\sum }^{}}}\,\alpha _{{{q}_{1}}{{q}_{2}}}^{{{n}_{1}}{{n}_{2}}}{{\left| 0 \right\rangle}_{i}}{{\left| 0 \right\rangle}_{j}}\left| {{\mathbf{1}}_{{{q}_{1}}}}{{\mathbf{1}}_{{{q}_{2}}}} \right\rangle, \\ \end{eqnarray} \tag{ 7 }$$

where ${{\left| l \right\rangle}_{i}}$ denotes the state of the system at site i [following the notation of (6)]; $\left| \mathbf{0} \right\rangle$ denotes the vacuum for the continuum and $\left| {{\mathbf{1}}_{{{q}_{1}}}}{{\mathbf{1}}_{{{q}_{2}}}} \right\rangle$ denotes the state with both particles in the continuum corresponding to the quantum numbers ${{n}_{1}}{{q}_{1}}$ and ${{n}_{2}}{{q}_{2}}$. The time-dependent Schrödinger equation for $\left| \Phi \right\rangle$ leads to a set of coupled equations for probability amplitudes ${{C}_{l}}$, grouped in a 6-component vector $\mathbf{C}$, corresponding to discrete states, as well as for continuum amplitudes $\alpha _{{{q}_{1}}{{q}_{2}}}^{{{n}_{1}}{{n}_{2}}}$.

$$\begin{eqnarray} i\mathbf{\dot{C}} & = & {{\mathbf{U}}_{1}}\mathbf{C}+\underset{n,{{q}_{1}}{{q}_{2}}}{\mathop{\mathop{\sum }^{}}}\,\left[ \mathbf{P}_{ij,{{q}_{1}}{{q}_{2}}}^{n} \right]\alpha _{{{q}_{1}}{{q}_{2}}}^{nn}i\dot{\alpha }_{{{q}_{1}}{{q}_{2}}}^{nn} \\ {} & = & \left[ {{E}_{n}}\left( {{q}_{1}} \right)+{{E}_{n}}\left( {{q}_{2}} \right) \right]\alpha _{{{q}_{1}}{{q}_{2}}}^{nn}+{{\left[ \mathbf{P}_{ij,{{q}_{1}}{{q}_{2}}}^{n} \right]}^{\dagger }}\mathbf{C} \\ {} & - & \frac{\pi D{{\Omega }_{\text{eff}}}}{12}\underset{{{q}_{3}},{{q}_{4}}}{\mathop{\mathop{\sum }^{}}}\,\alpha _{{{q}_{3}}{{q}_{4}}}^{nn}-\frac{\pi D{{\Omega }_{\text{eff}}}}{6}\underset{n\ne n\prime ,{{q}_{3}},{{q}_{4}}}{\mathop{\mathop{\sum }^{}}}\,\alpha _{{{q}_{3}}{{q}_{4}}}^{n\prime n\prime }, \\ \end{eqnarray} \tag{ 8 }$$

where $|i-j|=1$ and ${{\mathbf{U}}_{1}}$ is the interaction matrix between the discrete states originating from the Hamiltonian equation (4):

$$\begin{eqnarray}{{\mathbf{U}}_{1}}=\left( \begin{array}{cc} \mathbf{U} & \mathbf{0} \\ \mathbf{0} & \mathbf{U} \\ \end{array} \right)\end{eqnarray} \tag{ 9 }$$

with

$$\begin{eqnarray}\mathbf{U}=\left( \begin{array}{cc} {{U}_{\text{ss}}} & 0 & {{T}_{\text{ps}}} \\ 0 & {{E}_{1}}+{{U}_{\text{ps}}} & 0 \\ {{T}_{\text{ps}}} & 0 & 2{{E}_{1}}+{{U}_{\text{pp}}} \\ \end{array} \right).\end{eqnarray} \tag{ 10 }$$

Due to a lack of the direct coupling between the Wannier states at different sites, ${{\mathbf{U}}_{1}}$ is block diagonal. The Bloch−Wannier Hamiltonian in equation (5) will give rise to the discrete-continuum coupling array $\mathbf{P}_{ij,{{q}_{1}}{{q}_{2}}}^{n}={{\left[ \mathbf{P}_{i,{{q}_{1}}{{q}_{2}}}^{n},\mathbf{P}_{j,{{q}_{1}}{{q}_{2}}}^{n} \right]}^{T}}$.

To find the time evolution of the pair in the continuum, we make the ansatz that $\alpha _{{{q}_{1}}{{q}_{2}}}^{nn}\approx {{\alpha }^{n}}$. This is justified as the attractive interaction is momentum independent and much larger than the bandwidth of the each Bloch band n, so that the population amplitudes have weak momentum independence. Moreover, the last term in equation (8) denotes coupling of population amplitude of a Bloch band n to that of another Bloch band $n\prime$. The corresponding coupling strength $\sim D$ is of the same order of magnitude as the energy difference of the nearest Bloch bands which will be strongly coupled. Accordingly, we have assumed that for the last term in equation (8), $n-n\prime =\pm 1$ and ${{\alpha }^{n}}\approx {{\alpha }^{n-1}}\approx {{\alpha }^{n+1}}$. Within these approximations, one can rewrite equation (8) as,

$$\begin{eqnarray} i{{{\dot{\alpha }}}^{n}} & \approx & \left[ {{E}_{n}}\left( {{q}_{1}} \right)+{{E}_{n}}\left( {{q}_{2}} \right)-\pi D{{\Omega }_{\text{eff}}} \right]{{\alpha }^{n}} \\ {} & + & {{\left[ \mathbf{P}_{ij,{{q}_{1}}{{q}_{2}}}^{n} \right]}^{\dagger }}\mathbf{C}, \\ \end{eqnarray} \tag{ 11 }$$

where strong dipolar interaction effectively shifts the dispersion of each Bloch band. As initially the pairs were prepared in the discrete states in the limit of weak polarizing field, by performing Laplace transform of equation (11) we get,

$$\begin{eqnarray}&&{{\alpha }_{n}}\left( s \right)=-i\frac{{{\left[ \mathbf{P}_{ij,{{q}_{1}}{{q}_{2}}}^{n} \right]}^{\dagger }}\mathbf{C}}{s-\left[ {{E}_{n}}\left( {{q}_{1}} \right)+{{E}_{n}}\left( {{q}_{2}} \right)-\pi D{{\Omega }_{\text{eff}}} \right]}.\end{eqnarray} \tag{ 12 }$$

Then one can do similar Laplace transform for the discrete state amplitudes in equation (8) and eliminate the continuum amplitudes by equation (12). Subsequently, in the time evolution of the discrete state amplitudes, one gets expressions like,

$$\begin{eqnarray}&&\underset{n \gt {{n}_{0}},{{q}_{1}}}{\overset{{{n}_{\text{cut}}}}{{q}_{2}}}{\mathop{\mathop{\sum }^{}}}\,\left[ \mathbf{P}_{ij,{{q}_{1}}{{q}_{2}}}^{n} \right]\frac{{{\left[ \mathbf{P}_{ij,{{q}_{1}}{{q}_{2}}}^{n} \right]}^{\dagger }}\mathbf{C}}{s-\mathbf{i}\left[ {{E}_{n}}\left( {{q}_{1}} \right)+{{E}_{n}}\left( {{q}_{2}} \right)-\pi D{{\Omega }_{\text{eff}}} \right]},\end{eqnarray} \tag{ 13 }$$

where we have introduced a cut off ${{n}_{\text{cut}}}\sim 20$ in the band index and $\mathbf{i}=\sqrt{-1}$. Any excitations to higher bands than ${{n}_{\text{cut}}}$, will be lost due to formation of strongly bound molecular pairs (The origin of this cut off—the abundance of sticking collisions [37]—is discussed in detail in appendix B). Due to the shift of the energy of the continuum, the minimum of the continuum energy, ${{E}_{{{n}_{0}}}}-\pi D{{\Omega }_{\text{eff}}}/20$. Then by transforming the summation over energy level n to integration, one integrates over the range $-\infty \to \infty$.

The procedure described above takes into account the continuum-continuum transitions in a mean-field way. In effect, we obtain the effective coupled equations for the time evolution of discrete Wannier states amplitudes $\mathbf{\dot{C}}=\mathcal{M}\mathbf{C}$. The coupling matrix $\mathcal{M}$ is expressed as

$$\begin{eqnarray} \mathcal{M} & = & -\left[ i{{\mathbf{U}}_{1}}+\frac{\pi }{2}{{\left[ \frac{D{{\Omega }_{\text{eff}}}}{3\pi } \right]}^{2}}\Gamma \right], \\ \Gamma & = & \underset{n}{\mathop{\mathop{\sum }^{}}}\,\mathop{\int }^{}\mathop{\int }^{}d{{q}_{1}}d{{q}_{2}}\left[ \mathbf{P}_{{{q}_{1}}{{q}_{2}}}^{n} \right]{{\left[ \mathbf{P}_{{{q}_{1}}{{q}_{2}}}^{n} \right]}^{\dagger }}, \\ \end{eqnarray} \tag{ 14 }$$

where we have introduced the decay matrix $\Gamma$ and the effective trapping strength ${{\Omega }_{\text{eff}}}=\hbar \Omega /2{{E}_{R}}$. In the expression above $\mathbf{P}_{{{q}_{1}}{{q}_{2}}}^{n}$ is a vector of couplings of 6 Wannier discrete states (3 per site—compare equation (6)) with the continuum. The non-zero elements linking different sites of the discrete-continuum coupling array will induce an additional effective hopping terms for the pairs from site i to site j. One immediately notices that in the absence of interference effects, the decay rate of each channel will be proportional to ${{D}^{2}}$ (compare equation (14)). Thus deviation from this behaviour may serve as an indicator of important interference terms affecting the dynamics.

The full time dependent solution of the problem now reads $\mathbf{C}\left( t \right)=\mathop{\sum }^{}_{l=1}^{6}{{c}_{l}}\exp \left[ -{{\Gamma }_{l}}t-i{{\epsilon }_{l}} \right]{{\mathbf{u}}_{l}}$, where ${{\mathbf{u}}_{l}}$ is the eigenvector of the matrix $\mathcal{M}$ with ${{\Gamma }_{l}}$ and ${{\epsilon }_{l}}$ being the decay rate and the energy of the l-th eigenstate for the neighbouring sites. In figure 2 (left panel) we plot the decays rates for two neighbouring sites $|i-j|=1$. Let us concentrate on the states with the low decay rates (black line and black-circled line). All the other channels (denoted by red and blue curves) have decay rates proportional to ${{D}^{2}}$, which points towards absence of interference effects. The states with low decay rates show a much different and slower scaling as a function of D. The corresponding eigenstates can be approximately expressed as symmetric and anti-symmetric combinations of the single-site eigenstates ${{\left| \pm \right\rangle}_{ij}}=\left( {{\left| \phi \right\rangle}_{i}}\pm \left| \phi \right\rangle \right)/\sqrt{2}$ with energies ${{\epsilon }_{\pm }}$ with ${{\epsilon }_{+}}<{{\epsilon }_{-}}$:

$$\begin{eqnarray}&&\left| \phi \right\rangle_{i}=\left[ {\beta }_{1}\left( \hat{s}\,_{i}^{\dagger } \right)^{2}+{\beta }_{2}\left( \hat{p}\,_{i}^{\dagger } \right)^{2} \right]\left| 0 \right\rangle\end{eqnarray} \tag{ 15 }$$

expressed in terms of s and p orbitals. The overlap of these approximate combinations with the exact eigenstates: $\left| \left\langle \sigma \right|\sigma \prime \right\rangle_{exact}\left| _{i}\right.\;j\approx \mathcal{F}{\delta }_{\sigma {{\sigma }^{\prime }}}$ is large with $\mathcal{F}\sim 0.95$, where $\sigma ,\sigma \prime =\pm$. The deviation from the perfect overlap is due to the fact that there is an additional continuum induced off-site transition between states with opposite parity, $\left| \phi \right\rangle_{i}\leftrightarrow sign\left( i-j \right)\hat{s}_{j}^{\dagger }\hat{p}\,_{j}^{\dagger }\left| 0 \right\rangle$.

Zoom In Zoom Out Reset image size

The state with the lowest decay rate (black line in figure 2 (left panel)) corresponds to the state ${{\left| - \right\rangle}_{ij}}$ with highest energy. For this state we find that the ratio between the decay rate and the energy lies in the range, ${{\Gamma }_{-}}/{{\epsilon }_{-}}=0.01\to .05$ as the dipolar strength changes from $10\to 50$. On the other hand, for the state ${{\left| + \right\rangle}_{ij}}$, for the same dipolar range, ${{\Gamma }_{+}}/{{\epsilon }_{+}}=0.05\to 0.1$.

It follows that on the timescale of $\sim 1/{{\Gamma }_{+}}\approx 10/{{E}_{R}}$, only the ${{\left| - \right\rangle}_{ij}}$ survives and will be populated. What is the origin of this surprizing stabilization? What slows down the decay in such a spectacular way? A clue lies in the fact that the analogous analysis of the fate of a pair localized in a single site only indicates a much faster decay (see figure 2 (right panel)). Therefore, we find a surprising situation in which a state is stabilized by delocalizing between two neighbouring sites in the presence of continuum-induced tunneling—a coupling between sites. Such a situation is well known from single bound electron quantum optics studies—it is the phenomenon of population trapping [20]. While the physics seems to be quite similar to a strong laser field induced trapping [20] let us stress that the 'dark state' in our situation entangles two distinct lattice sites. We like to point out that in our scenario both the decay and delocalization is induced by strong coupling to the continuum. Similar analysis may be carried out for separated sites with $|i-j|>1$. It shows that in that case the effect of the continuum-assisted coupling is much smaller within the regime of dipolar strengths studied.

Next we discuss the creation of dimer states due to population trapping for the half-filling of the pairs. It is known that the strong dipolar interaction induces a density-wave phase where the pairs arrange in a checkerboard pattern [40]. As such pairs are pinned to the sites, the checkerboard configuration will not be stable as each pair occupied site will decay rapidly to the continuum. The stable configuration can only have states containing the delocalized state ${{\left| - \right\rangle}_{ij}}$. Then, in the limit of strong interaction and for half-filling of pairs, ${{\left| - \right\rangle}_{ij}}$ will cover the whole region of lattice sites. The resulting many-body state is a checkerboard state of nearest-neighbour dimers, ${{\left| \Psi \right\rangle}_{A}}={{\Pi }_{i}}{{\left| - \right\rangle}_{2i,2i+1}}$ or ${{\left| \Psi \right\rangle}_{B}}={{\Pi }_{i}}{{\left| - \right\rangle}_{2i-1,2i}}$. These dimer states are the ground states of the celebrated Majumder−Ghosh (MG) model [21]. This paradigmatic model consists of a frustrated one-dimensional spin chain consisting of nearest and next-nearest neighbour hopping with a particular ratio. The dimer state is characterized by an absence of long-range correlations, ${\left\langle\hat{b}\,_{i}^{\dagger }{{\hat{b}\,}_{j}}\right\rangle}_{\Psi }=\left\langle{{\hat{n}\,}_{i}}{{\hat{n}\,}_{j}}\right\rangle=0$ for $|i-j|>1$. This dimerized state can be thought of as the simplest form of the valance-bond solid with short-range correlations and with double the period of the original lattice. As a further support for our claim, in appendix B, we have presented many-body calculation for small systems which shows that the state with lowest decay has almost unit overlap with the MG state.

To prepare the MG state, initially one prepares half-filling molecular repulsively bound pairs in the regime of low dipolar interaction. Then one can switch on the strong electric field to create a strong dipolar interaction. This couples the Wannier states to the continuum. Such a coupling usually reduces the population of molecules in the Wannier states. But in our case, due to the coherent population trapping, the initial particle density in the Wannier states will be maintained within the decay time of the population-trapped state. Any small deviation of the initial density from half-filling will manifest themselves as excitations to the final MG state.

The doubling of periodicity in a MG state can form an experimental signature in the time of flight image due to the reduction of the Brillouin zone. The required temperature to reach this phase depends on the delocalization energy which is given by the energy difference $\delta E$ between the single-site state ${{\left| \phi \right\rangle}_{i}}$ and the dimer state ${{\left| - \right\rangle}_{ij}}$. For a dipolar strength of $D\sim 20$ (near the lowest decay rate in figure 2) this energy difference is of the order of $0.4{{E}_{R}}$. For RbCs molecules, these parameters correspond to a dipole moment of $\sim 0.7$ Debye with a lattice constant $\sim 500$ nm. Then the relevant temperature scale to observe this phase is $\sim 50$ nK. Such a temperature is much larger than the one needed to reach the super-exchange regime for the ultracold atoms, and thus it is much easier to access experimentally. The price to pay in our present case is the meta-stability of the dimerized state with the lifetime $\sim 10$ ms. One way to increase the stability is by decreasing the electric field strength within the decay time, which makes all the interaction terms small. At the same time, by increasing the lattice depth one can decrease the tunneling amplitudes. This will make the dimer state frozen in time felicitating the characterization of it.

Let us extend our calculation of a single pair distributed in two sites to a larger system size. We have performed an exact diagonalization for half-filled pairs distributed over 8 sites. Following the same procedure, we have found an effective equation of motion for the many-body discrete state probability amplitudes denoted by ${{\mathbf{C}}_{\text{mb}}}$ with modified continuum induced transition matrix ${{\mathbf{P}}_{\text{mb}}}$ where we have taken into account continuum induced long-range coupling. The resulting equation of motion has the form, ${{\mathbf{\dot{C}}}_{\text{mb}}}={{\mathcal{M}}_{\text{mb}}}{{\mathbf{C}}_{\text{mb}}}$, and the many body coupling matrix ${{\mathcal{M}}_{\text{mb}}}$ is given by

$$\begin{eqnarray}&&{{\mathcal{M}}_{\text{mb}}}=-\left[ i{{\mathbf{U}}_{\text{mb}}}+\frac{\pi }{2}{{\left[ \frac{D{{\Omega }_{\text{eff}}}}{3\pi } \right]}^{2}}{{\mathbf{P}}_{\text{mb}}} \right],\end{eqnarray} \tag{ 16 }$$

where the discrete states interaction matrix ${{\mathbf{U}}_{\text{mb}}}$ now also includes the long-range dipolar interaction. We then find the eigenvalues and eigenstates of the matrix ${{\mathcal{M}}_{\text{mb}}}$. The real part of the eigenvalues describe the decay rate of the respective eigenstates. We then concentrate on the state with the lowest decay rate which shows similar decrease in decay strength as the PT state discussed in the manuscript. Next, we find the overlap of this state with the Majumdar–Ghosh (MG) states (${{\left| \Psi \right\rangle}_{A}},{{\left| \Psi \right\rangle}_{B}}$) defined in the manuscript. We find that for larger dipolar strength D, there is large overlap of the lowest decay state with the antisymmetric MG state ${{\left| \Psi \right\rangle}_{-}}=\left[ {{\left| \Psi \right\rangle}_{A}}-{{\left| \Psi \right\rangle}_{B}} \right]/\sqrt{2}$. We denote this overlap by the function ${{\mathcal{F}}_{MG}}$ and plot it against the dipolar strength in figure 3. Manifestly, in spite of a strong repulsive long-range interaction between the off-site molecular pairs, the phenomenon of PT can result in a creation of the frustrated MG state. Such a result shows, furthermore, a possibility of a creation of resonating valence-bond MG state. A detailed discussion of such a possibility is beyond the scope of the present paper.

Zoom In Zoom Out Reset image size

In this section we discuss a possible way to construct an effective Hamiltonian in terms of local operators for low density of the pairs. To do this, we consider a simple system where one pair of atoms is moving in three sites coupled to the continuum. Following the same procedure as before we derive the full coupling matrix $\mathcal{M}\left( i,i+1,i+2 \right)$ for three sites. Studying eigenstates related to the lowest decay rates we find, as before, that the coherent population trapping occurs due to the coupling of neighbouring sites via ${{\left| \pm \right\rangle}_{ij}}$ states. Subsequently, a tunneling Hamiltonian in terms of the states ${{\left| \pm \right\rangle}_{ij}}$ is given by, ${{H}_{i,i+1,i+2}}=-{{J}_{\text{eff}}}\left( {{\left| - \right\rangle}_{i,i+1}}{{\left\langle - \right|}_{i+1,i+2}}+\alpha {{\left| + \right\rangle}_{i,i+1}}{{\left\langle + \right|}_{i+1,i+2}}+h.c. \right)$, where $\alpha$ can be extracted from the eigenvalues of the effective coupling matrix $\mathcal{M}\left( i,i+1,i+2 \right)$. As the states ${{\left| - \right\rangle}_{i,i+1}},{{\left| \pm \right\rangle}_{i+1,i+2}}$ are not orthogonal, it is convenient to rewrite ${{H}_{i,i+1,i+2}}$ in terms of local orthogonal operators. To do that we define a local pair operator, ${{\left| \phi \right\rangle}_{i}}=b_{i}^{\dagger }\left| 0 \right\rangle$ which creates a pair at site i in the lowest decay state. The pair operators satisfy bosonic commutation relations $\left[ {{b}_{i}},b_{j}^{\dagger } \right]={{\delta }_{ij}}$. In terms of these pair operators we can rewrite the states as ${{\left| \pm \right\rangle}_{ij}}=\frac{1}{\sqrt{2}}\left[ b_{i}^{\dagger }\pm b_{j}^{\dagger } \right]\left| 0 \right\rangle$. Subsequently, the Hamiltonian ${{H}_{i,i+1,i+2}}$ is re-expressed as,

$$\begin{eqnarray} \frac{{{H}_{i,i+1,i+2}}}{{{J}_{\text{eff}}}} & = & -\frac{1+\alpha }{2}\left[ b_{i}^{\dagger }{{b}_{i+1}}+b_{i+1}^{\dagger }{{b}_{i+2}}+h.c \right]+\left( 1-\alpha \right)b_{i+1}^{\dagger }{{b}_{i+1}} \\ {} & {} & +\frac{1-\alpha }{2}\left[ b_{i}^{\dagger }{{b}_{i+2}}+b_{i+2}^{\dagger }{{b}_{i}} \right]. \\ \end{eqnarray} \tag{ 17 }$$

The values of ${{J}_{\text{eff}}}$ and α are derived by comparing the energies of Hamiltonian equation ( 17) and the energies of the states with three lowest decay rates derived from the full coupling matrix $\mathcal{M}\left( i,i+1,i+2 \right)$. For small values of the dipolar strength D we find that $\alpha \approx 1$, thus the long-range tunneling is small and one recovers the usual picture with nearest-neighbour tunneling only. But for higher dipolar strengths $\alpha \ne 1$, due to the different decay rates of $\left| \pm \right\rangle$ states. For such values of α one obtains, therefore, (compare equation (17)) an effective model with next-nearest neighbour tunneling leading to frustration. We would like to point out that, in the present situation, the origin of such a frustration is entirely different from the usual origin of such terms due to the higher order processes in solid-state systems [23].

At this point, we write down the effective many-body Hamiltonian including long-range dipolar interaction and involving all sites as,

$$\begin{eqnarray} {{H}_{\text{eff}}} & = & \underset{i}{\mathop{\mathop{\sum }^{}}}\,{{H}_{i,i+1,i+2}}+\frac{D}{{{\pi }^{3}}}\underset{ij}{\mathop{\mathop{\sum }^{}}}\,\frac{{{n}_{i}}{{n}_{j}}}{|i-j{{|}^{3}}}-\mu \underset{i}{\mathop{\mathop{\sum }^{}}}\,{{n}_{i}} \\ {} & = & -{{J}_{\text{eff}}}\left( 1+\alpha \right)\underset{\langle ij\rangle }{\mathop{\mathop{\sum }^{}}}\,b_{i}^{\dagger }{{b}_{j}}+{{J}_{\text{eff}}}\frac{1-\alpha }{2}\underset{\langle \langle ij\rangle \rangle }{\mathop{\mathop{\sum }^{}}}\,b_{i}^{\dagger }{{b}_{j}} \\ {} & + & \frac{2D}{{{\pi }^{3}}}\underset{i\ne j}{\mathop{\mathop{\sum }^{}}}\,\frac{{{n}_{i}}{{n}_{j}}}{|i-j{{|}^{3}}}-\mu \underset{i}{\mathop{\mathop{\sum }^{}}}\,{{n}_{i}}, \\ \end{eqnarray} \tag{ 18 }$$

where we have introduced the chemical potential μ for the pairs and $\langle \langle ij\rangle \rangle$ is a shorthand for next nearest neighbour summation index. The Hamiltonian in equation (18) contains two sources of frustration: (i) the effective next-nearest neighbour tunneling, and (ii) long-range dipolar interaction. For our present system, the deviation of dipolar interaction from the cubic power law is negligible [39]. The Hamiltonian, equation (18), is a generalization of the ${{J}_{1}}-{{J}_{2}}$ model where the interaction is present to the next-nearest neighbours only. The ${{J}_{1}}-{{J}_{2}}$ model is a prototype for studying the effect of frustration and emergence of various proposed exotic phases in magnetic materials [1]. The single particle dispersion relation for this Hamiltonian is given by ${{\epsilon }_{q}}={{J}_{\text{eff}}}\left( 1+\alpha \right)\cos qa+{{J}_{\text{eff}}}\frac{1-\alpha }{2}\cos 2qa$. For $\frac{1-\alpha }{1+\alpha }>1/2$ it shows two minima at wavevectors $\pm Qa={{\cos }^{-1}}\left[ -\frac{1+\alpha }{2\left( 1-\alpha \right)} \right]$. In our case, the two-minima limit corresponds to $D>18$. In the low-density limit, one way to treat the problem is by going to the two-component homogenous Bose gas limit [39] with the effective Hamiltonian,

$$\begin{eqnarray}&&{{H}_{\text{eff}}}=\mathop{\int }^{}\left[ \frac{1}{2}{{T}_{1}}\left( \rho _{1}^{2}+\rho _{2}^{2} \right)+{{T}_{12}}{{\rho }_{1}}{{\rho }_{2}}-\mu \left( {{\rho }_{1}}+{{\rho }_{2}} \right) \right]dx,\end{eqnarray} \tag{ 19 }$$

where ${{\rho }_{1,2}}$ are the densities of the two component Bose gas centered around the the minima $\pm Q$ and ${{T}_{1}},{{T}_{12}}$ are the renormalized intra-component and inter-component interaction. A detailed discussion of the Hamiltonian equation (19) is presented in the methods section. For a short-range ${{J}_{1}}-{{J}_{2}}$ model, the phase diagram from such a procedure shows qualitative agreement with more involved density-matrix renormalization group simulations [39]. When ${{T}_{1}}<{{T}_{12}}$, the mean-field ground state solution is given by the phase-separated state ${{\rho }_{1}}\ne 0,{{\rho }_{2}}=0$ or ${{\rho }_{1}}=0,{{\rho }_{2}}\ne 0$. Choosing one of the ground state will break the discrete symmetry which will result in true long-range order (LRO) even in one-dimension. The nature of this phase can readily be observed by writing the wavefunction in phase space, ${{\psi }_{s}}=\sqrt{{{\rho }_{s}}}\exp \left[ -i{{\theta }_{s}} \right]$, with s = 1, 2. When ${{\rho }_{1}}=0$, we see that $\left\langle \hat{b}\,_{i}^{\dagger }\right\rangle=\sqrt{{{\rho }_{1}}}\exp \left[ -iQx+\theta \right]$. Such a 'cone' phase is identified as a the vector-chiral (VC) phase which breaks the ${{\mathcal{Z}}_{2}}$ symmetry. In contrast when ${{T}_{1}}>{{T}_{12}}>0$, we have a mixed state with equal density from both components. This homogeneous solution with ${{\rho }_{1}}={{\rho }_{2}}$ is known as the two-component Tomonaga−Luttinger (TLL$_{2}$) liquid. There can be another possibility when the effective inter-species interaction is attractive ${{T}_{1}}>0,{{T}_{12}}<0$. In this situation, intra-component bound states with emerge with center of mass momentum $\sim 2Q$. Such bound states with finite momenta are usually not present in the anti-ferromagnetic model [39]. In the present case, these bound states are a direct consequence of the long-range nature of the dipolar interaction which can induce resonances [41]. The quasi-condensate of such bound pairs can give rise to a spin-nematic phase [42, 43], or spin-density wave phase [44], a detailed discussion of which is beyond the scope of current article. In figure 4, we have plotted the resulting phase diagram in the $D-\mu$ parameter space for vanishingly small μ. We find that the vector-chiral phase is stable for smaller and larger values of the dipolar strength D. In between the homogeneous TLL$_{2}$ phase is the ground state. For larger values of chemical potential μ, one finds that there is a bound state phase due to ${{T}_{12}}<0$.

Zoom In Zoom Out Reset image size

Here we write down the Hamiltonian originating from the single-particle kinetic energy and dipolar interaction between the discrete states,

$$\begin{eqnarray}&&{{H}_{\text{d}isc}}={{H}_{\text{T}}}+{{H}_{\text{pair}}}+{{H}_{\operatorname{int}}}\end{eqnarray} \tag{ A.3 }$$

with ${{H}_{T}}$ describing standard and interaction induced (density-dependent) single particle tunneling terms

$$\begin{eqnarray} H_{\rm T}&=& \sum_{\langle ij\rangle}\left [ -J_0\hat{s}^{\dagger}_i\hat{s}_j+J_1 \hat{p}^{\dagger}_i\hat{p}_j \right ] + \sum_{\langle ij\rangle} \left [ T_0\hat{s}^{\dagger}_i \left(\hat{n}_{s i}+\hat{n}_{s j}\right)\hat{s}_j+T_1 \hat{p}^{\dagger}_i\left(\hat{n}_{p i}+\hat{n}_{p j}\right)\hat{p}_j \right ] \notag\\ &&+ \sum_{\langle ij\rangle} \left [ T_{00}\hat{s}^{\dagger}_i \left(\hat{n}_{p i}+\hat{n}_{p j}\right)\hat{s}_j \right\rangle + \left\langle T_{11} \hat{p}^{\dagger}_i\left(\hat{n}_{s i}+\hat{n}_{s j}\right)\hat{p}_j \right ] + T_{01}\sum_{\langle ij\rangle}{\it f}_{ij} \left[\hat{p}^{\dagger}_i\hat{n}_{s i}\hat{s}_j + h.c. \right ]\notag\\ &&+ T_{10}\sum_{\langle ij\rangle}{\it f}_{ij} \left[\hat{p}^{\dagger}_i\hat{n}_{p j}\hat{s}_j + h.c. \right ] + {T^{\prime}}_{01}\sum_{\langle ij\rangle}{\it f}_{ij} \left[\hat{p}^{\dagger}_j\hat{n}_{s i}\hat{s}_i + h.c. \right ]\notag\\ &&+{T^{\prime}}_{10}\sum_{\langle ij\rangle}{\it f}_{ij} \left[\hat{p}^{\dagger}_j\hat{n}_{p j}\hat{s}_i + h.c. \right ] \end{eqnarray} \tag{ A.4 }$$

while the correlated pair hopping part of the Hamiltonian reads

$$\begin{eqnarray} &&{H}_{\text{pair}}=\sum_{ij}\,\left[ \frac{1}{2}{T}_{\text{p},0}\hat{s}\,_{i}^{\dagger }\hat{s}\,_{i}^{\dagger }\hat{s}\,_{j}\hat{s}\,_{j}+\frac{1}{2}{T}_{\text{p},1}\hat{p}\,_{i}^{\dagger }\hat{p}\,_{i}^{\dagger }\hat{p}\,_{j}\hat{p}\,_{j}+\frac{1}{2}{T}_{\text{p},01}\left( \hat{s}\,_{i}^{\dagger }\hat{s}\,_{i}^{\dagger }\hat{p}\,_{j}\hat{p}\,_{j}+h.c. \right)+{T}_{\text{p},10}\hat{s}\,_{i}^{\dagger }\hat{p}\,_{i}^{\dagger }\hat{p}\,_{j}\hat{s}\,_{j} \right]\end{eqnarray} \tag{ A.5 }$$

where ${{J}_{0}},{{J}_{1}}>0$ denote the single particle nearest neighbor tunneling amplitudes in the $s,p$-orbital respectively. Intra-orbital interaction-induced tunneling amplitudes are denoted by ${{T}_{0}},{{T}_{1}},{{T}_{00}},{{T}_{11}}$. The interaction-induced inter-orbital tunneling amplitudes are given by ${{T}_{01}},T{{\prime }_{01}},{{T}_{10}},T{{\prime }_{10}}$. The staggered nature of the inter-orbital tunneling is denoted by ${{f}_{ij}}=\pm 1$ when $i-j=\mp 1$. The pair tunneling Hamiltonian is denoted by ${{H}_{\text{pair}}}$ and the corresponding pair-tunneling amplitudes are given by ${{T}_{\text{p},0}},{{T}_{\text{p},1}},{{T}_{\text{p},01}},{{T}_{\text{p},10}}$. All these terms, partially canceling each other, are neglected in our simplified Hamiltonian equation (3). The extended analysis taking into account single particle tunneling is discussed in appendix B, below.

Next, we rewrite the dipolar interaction between the various Wannier orbitals from equation(3) in the main text,

$$\begin{eqnarray} H_{\operatorname{int}} & = & \underset{i,\sigma =s,p}{\mathop{\mathop{\sum }^{}}}\,\frac{U_{\sigma \sigma }}{2}\hat{n}\,_{\sigma i}\left( \hat{n}\,_{\sigma i}-1 \right)+U_{\text{ps}}\underset{i}{\mathop{\mathop{\sum }^{}}}\,\hat{n}\,_{si}\hat{n}\,_{pi} \\ && + \frac{T_{\text{ps}}}{2}\underset{i}{\mathop{\mathop{\sum }^{}}}\,\left[ \hat{p}\,_{i}^{\dagger }\hat{p}\,_{i}^{\dagger }\hat{s}\,_{i}\hat{s}\,_{i}+H.c \right]+\frac{D}{2{{\pi }^{3}}}\underset{\sigma ,\sigma \prime ,i\ne j}{\mathop{\mathop{\sum }^{}}}\,\frac{\hat{n}\,_{\sigma i}\hat{n}\,_{\sigma \prime j}}{|i-j|^{3}}. \end{eqnarray} \tag{ A.6 }$$

The different amplitudes in the discrete subspace Hamiltonian are obtained by appropriate integrals of the dipole−dipole interaction potential and the mode functions (compare equation (2)) that contain Wannier functions for orbitals along x with product of ground state Gaussians in perpendicular direction. For completness we list these integrals explicitly, assuming a shorthand notation

$$\begin{eqnarray*}&&\mathcal{W}\left( \mathbf{r},\mathbf{r}\prime \right)={{V}_{\text{dd}}}\left( \mathbf{r}-\mathbf{r}\prime \right)\phi _{0}^{2}\left( y \right)\phi _{0}^{2}\left( y\prime \right)\phi _{0}^{2}\left( z \right)\phi _{0}^{2}\left( z\prime \right)\end{eqnarray*}$$

and assuming the Wannier functions to be real:

$$\begin{eqnarray}&&{{U}_{\text{ss}}}=\mathop{\int }^{}d\mathbf{r}d\mathbf{r}\prime |\omega _{i}^{s}\left( x \right)\omega _{i}^{s}\left( x\prime \right){{|}^{2}}\mathcal{W}\left( \mathbf{r},\mathbf{r}\prime \right),\end{eqnarray} \tag{ A.7 }$$

$$\begin{eqnarray} U_{\rm pp}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&|\omega^{p}_i\left(x\right)\omega^{p}_i\left(x^{\prime}\right)|^2 \mathcal{W}\left(\mathbf{r,r^{\prime}}\right),\notag\\ T_{0}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\left[\omega^{s}_j\left(x\right)\right]^3\omega^{s}_i\left(x^{\prime}\right) \mathcal{W}\left(\mathbf{r,r^{\prime}}\right),\notag\\ T_{1}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\left[\omega^{p}_j\left(x\right)\right]^3\omega^{p}_i\left(x^{\prime}\right)\mathcal{W}\left(\mathbf{r,r^{\prime}}\right), \notag\\ T^{\prime}_{01}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\omega^{p}_j\left(x\right)\omega^{s}_i\left(x\right)\left[\omega^{s}_i\left(x^{\prime}\right)\right]^2\mathcal{W}\left(\mathbf{r,r^{\prime}}\right), \notag\\ T^{\prime}_{10}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\omega^{p}_j\left(x\right)\omega^{s}_i\left(x\right)\left[\omega^{p}_i\left(x^{\prime}\right)\right]^2 \mathcal{W}\left(\mathbf{r,r^{\prime}}\right), \notag\\ T_{\rm ps}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\omega^{p}_i\left(x\right)\omega^{p}_i\left(x^{\prime}\right)\omega^{s}_i\left(x\right)\omega^{s}_i\left(x^{\prime}\right) \mathcal{W}\left(\mathbf{r,r^{\prime}}\right),\notag\\ T_{\rm p,0}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\omega^{s}_i\left(x\right)\omega^{s}_j\left(x\right)\omega^{s}_i\left(x^{\prime}\right)\omega^{s}_i\left(x^{\prime}\right) \mathcal{W}\left(\mathbf{r,r^{\prime}}\right),\notag\\ T_{\rm p,1}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\omega^{p}_i\left(x\right)\omega^{p}_j\left(x\right)\omega^{p}_i\left(x^{\prime}\right)\omega^{p}_i\left(x^{\prime}\right)\mathcal{W}\left(\mathbf{r,r^{\prime}}\right),\notag\\ T_{\rm p,01}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\omega^{s}_i\left(x\right)\omega^{p}_j\left(x\right)\omega^{s}_i\left(x^{\prime}\right)\omega^{p}_j\left(x^{\prime}\right) \mathcal{W}\left(\mathbf{r,r^{\prime}}\right),\notag\\ T_{\rm p,10}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\left [ \omega^{s}_i\left(x\right)\omega^{s}_j\left(x\right)\omega^{p}_i\left(x^{\prime}\right)\omega^{p}_j\left(x^{\prime}\right) + \omega^{s}_i\left(x\right)\omega^{p}_j\left(x\right) \omega^{p}_i\left(x^{\prime}\right)\omega^{s}_j\left(x^{\prime}\right)\right]\mathcal{W}\left(\mathbf{r,r^{\prime}}\right), \notag\\ U_{\rm ps}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}&\left [ |\omega^{p}_i\left(x\right)\omega^{s}_i\left(x^{\prime}\right)|^2 + \omega^{p}_i\left(x\right)\omega^{s}_i\left(x\right)\omega^{p}_i\left(x^{\prime}\right)\omega^{s}_i\left(x^{\prime}\right)\right] \mathcal{W}\left(\mathbf{r,r^{\prime}}\right),\notag\\ T_{00}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}& \left[\omega^{s}_i\left(x\right)\omega^{s}_j\left(x\right)\left[\omega^{p}_i\left(x^{\prime}\right)\right]^2 + \omega^{s}_i\left(x\right)\omega^{p}_i\left(x\right)\omega^{p}_i\left(x^{\prime}\right)\omega^{s}_j\left(x^{\prime}\right)\right] \mathcal{W}\left(\mathbf{r,r^{\prime}}\right), \notag\\ T_{11}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}& \left[\omega^{p}_i\left(x\right)\omega^{p}_j\left(x\right)\left[\omega^{s}_i\left(x^{\prime}\right)\right]^2 + \omega^{p}_i\left(x\right)\omega^{s}_i\left(x\right)\omega^{s}_i\left(x^{\prime}\right)\omega^{p}_j\left(x^{\prime}\right)\right] \mathcal{W}\left(\mathbf{r,r^{\prime}}\right), \notag\\ T_{01}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}& \left[\omega^{p}_i\left(x\right)\omega^{s}_j\left(x\right)\left[\omega^{s}_i\left(x^{\prime}\right)\right]^2 + \omega^{p}_i\left(x\right)\omega^{s}_i\left(x\right)\omega^{s}_i\left(x^{\prime}\right)\omega^{s}_j\left(x^{\prime}\right)\right]\mathcal{W}\left(\mathbf{r,r^{\prime}}\right), \notag\\ T_{10}&=&\int d\mathbf{r}d\mathbf{r^{\prime}}& \left[\omega^{p}_i\left(x\right)\omega^{s}_j\left(x\right)\left[\omega^{s}_j\left(x^{\prime}\right)\right]^2 + \omega^{p}_i\left(x\right)\omega^{p}_j\left(x\right)\omega^{p}_j\left(x^{\prime}\right)\omega^{s}_j\left(x^{\prime}\right)\right]\mathcal{W}\left(\mathbf{r,r^{\prime}}\right) , \notag\\ \end{eqnarray} \tag{ A.8 }$$

Next we consider relevant properties of Bloch states. Let us denote the Bloch band ${n_{0}}$ as the start of the continuous bands with ${{E}_{{{n}_{0}}}}$ as the minimum of energy and ${{E}_{{{n}_{0}}}}\gg 1$. Then the energy of the Bloch band $n={{n}_{0}}+m$ can be written as, ${{E}_{n={{n}_{0}}+m}}\left( q \right)={{E}_{{{n}_{0}}}}+2{{n}_{0}}m+{{m}^{2}}+2\left( {{n}_{0}}+m \right)|q|+{{q}^{2}}$ for even m and one can get similar results for odd m. Moreover, we write the Bloch wavefunctions in the ${{E}_{{{n}_{0}}}}\gg 1$ limit as [48], ${{\mathcal{U}}_{nq}}\left( x \right)\approx \sqrt{\frac{2}{L}}\exp \left[ iqx \right]\cos \left[ \sqrt{{{E}_{n}}-q{{.}^{2}}-V/2}x \right],$ for even n and ${{\mathcal{U}}_{nq}}\left( x \right)\approx \sqrt{\frac{2}{L}}\exp \left[ iqx \right]\sin \left[ \sqrt{{{E}_{n}}-q{{.}^{2}}-V/2}x \right]$ for n odd. As these functions are eigenstates, they are also orthogonal, $\mathop{\int }^{}U_{nq}^{*}\left( x \right)U_{mq\prime }^{*}\left( x \right)dx={{\delta }_{n,m}}{{\delta }_{q,{{q}^{\prime }}}}.$ Now we write down the discrete-continuum coupling matrix elements as,

$$\begin{eqnarray*} P_{i,{{p}_{1}}{{p}_{2}}}^{n}\left( {{q}_{1}}{{q}_{2}} \right) & = & \mathop{\int }^{}U_{n{{q}_{1}}}^{*}\left( x \right)U_{n{{q}_{2}}}^{*}\left( x\prime \right){{V}_{\text{dd}}}\left( \mathbf{r}-\mathbf{r}\prime \right)\omega _{i}^{{{p}_{1}}}\left( x\prime \right)\omega _{i}^{{{p}_{2}}}\left( x \right) \\ {} & \times & |{{\phi }_{0}}\left( z \right){{\phi }_{0}}\left( y \right){{|}^{2}}|{{\phi }_{0}}\left( z\prime \right){{\phi }_{0}}\left( y\prime \right){{|}^{2}}d\mathbf{r}d\mathbf{r}\prime , \\ {} & \approx & \frac{1}{4}\mathop{\int }^{}dk{{V}_{\text{dd}}}\left( k \right)\left[ \mathcal{W}_{i}^{{{\sigma }_{1}}}\left( k-{{q}_{1}}+\sqrt{{{E}_{{{n}_{1}}}}\left( {{q}_{1}} \right)} \right) \right. \\ {} & \times & \mathcal{W}_{i}^{{{\sigma }_{2}}}\left( -k-{{q}_{2}}-\sqrt{{{E}_{{{n}_{2}}}}\left( {{q}_{1}} \right)} \right) \\ {} & + & \mathcal{W}_{i}^{{{\sigma }_{1}}}\left( k-{{q}_{1}}-\sqrt{{{E}_{{{n}_{1}}}}\left( {{q}_{1}} \right)} \right) \\ {} & \times & \left. \mathcal{W}_{i}^{{{\sigma }_{2}}}\left( -k-{{q}_{2}}+\sqrt{{{E}_{{{n}_{2}}}}\left( {{q}_{1}} \right)} \right) \right], \\ \end{eqnarray*}$$

where $\mathcal{W}_{i}^{\sigma }\left( k \right)$ is the Fourier transform of the Wannier function $\omega _{i}^{\sigma }\left( x \right)$. In deriving the above form, we have used the orthogonality condition between the Bloch functions and assumed that ${{E}_{n}}\gg 1$. Additionally, in the Hamiltonian equation (3), we have neglected terms corresponding to processes like $\hat{a}\,_{n{{q}_{1}}}^{\dagger }\hat{s}\,_{i}^{\dagger }{{\hat{s}\,}_{i}}{{\hat{s}\,}_{i}}$ where one particle is coupled to the continuum. The transition amplitudes for such processes contains convolution sums of the form

$$\begin{eqnarray*}&&{{S}_{n}}\sim \mathop{\int }^{}{{\mathcal{W}}^{{{\sigma }_{1}}}}\left( k+q+\sqrt{{{E}_{n}}} \right){{\mathcal{W}}^{{{\sigma }_{2}}}}\left( -k \right){{V}_{\text{dd}}}\left( k \right)dk.\end{eqnarray*}$$

As ${{E}_{n}}\gg 1$, such terms are negligibly small. Thus we ignored them in comparison to the leading two-particle transition amplitudes.

The Hamiltonian for the continuum Bloch states reads

$$\begin{eqnarray} {{H}_{\text{Bloch}}} & \approx & \underset{n,q}{\mathop{\mathop{\sum }^{}}}\,{{E}_{n}}\left( q \right)\hat{a}\,_{nq}^{\dagger }{{\hat{a}\,}_{nq}}-\frac{\pi D{{\Omega }_{\text{eff}}}}{12}\underset{n,\mathbf{q}}{\mathop{\mathop{\sum }^{}}}\,\hat{a}\,_{n{{q}_{1}}}^{\dagger }\hat{a}\,_{n{{q}_{2}}}^{\dagger }{{\hat{a}\,}_{n{{q}_{3}}}}{{\hat{a}\,}_{n{{q}_{4}}}} \\ {} & - & \frac{\pi D{{\Omega }_{\text{eff}}}}{6}\underset{n\ne n\prime ,\mathbf{q}}{\mathop{\mathop{\sum }^{}}}\,\hat{a}\,_{n{{q}_{1}}}^{\dagger }\hat{a}\,_{n{{q}_{2}}}^{\dagger }{{\hat{a}\,}_{{{n}^{\prime }}{{q}_{3}}}}{{\hat{a}\,}_{{{n}^{\prime }}{{q}_{4}}}} \\ {} & - & \frac{\pi D{{\Omega }_{\text{eff}}}}{6}\underset{n\ne n\prime ,\mathbf{q}}{\mathop{\mathop{\sum }^{}}}\,\hat{a}\,_{n{{q}_{1}}}^{\dagger }\hat{a}\,_{n\prime {{q}_{2}}}^{\dagger }{{\hat{a}\,}_{n{{q}_{3}}}}{{\hat{a}\,}_{{{n}^{\prime }}{{q}_{4}}}}, \\ \end{eqnarray} \tag{ A.9 }$$

the continuous band index $n,n\prime >{{n}_{0}}$ and the momentum index $\mathbf{q}=\left[ {{q}_{1}},{{q}_{2}},{{q}_{3}},{{q}_{4}} \right]$. The second term in the Hamiltonian equation (A.9) denotes the dipolar interaction between the molecules in the same Bloch band n whereas the next term denotes the transition of pairs between two Bloch bands and the last term denotes interaction between molecules from different Bloch bands. We only include the leading terms whose strength is of the order of $\sim D$. Furthermore, from Hamiltonian equation (A.9), we notice that the interaction is strongly attractive in the higher Bloch bands and for strong interaction ($D\gg 1$), the dipolar strength can exceed the width of the first few continuous Bloch bands.

Consider the effect of dipolar interaction when many pairs decay into the continuum. Again, within each Bloch band, the dipolar attraction is larger than the respective bandwidth of the Bloch band. This suggests strong binding of the molecular pairs. To denote this we introduce a composite operator for the pairs,

$$\begin{eqnarray*}&&\hat{b}\,_{n}^{\dagger }=\frac{\mathop{\int }^{}\mathop{\int }^{}\hat{a}\,_{n{{q}_{1}}}^{\dagger }\hat{a}\,_{n{{q}_{2}}}^{\dagger }d{{q}_{1}}d{{q}_{2}}}{\mathop{\int }^{}\mathop{\int }^{}d{{q}_{1}}d{{q}_{2}}}.\end{eqnarray*}$$

As the molecules can scatter to any quasi-momentum state with equal strong probability, one can assume that each quasi-momentum level in the band n is at most occupied by one molecule. Then, in terms sof the pairing operator, one can find an momentum average representation Hamiltonian equation (A.9) in terms of the composite operators as,

$$\begin{eqnarray} {{H}_{\text{Bloch}}} & \approx & \underset{n}{\mathop{\mathop{\sum }^{}}}\,\left[ {{\epsilon }_{n,\text{avg}}}-\frac{\pi D{{\Omega }_{\text{eff}}}}{3} \right]\hat{b}\,_{n}^{\dagger }{{\hat{b}\,}_{n}}-\frac{2\pi D{{\Omega }_{\text{eff}}}}{3}\underset{n\ne {{n}^{\prime }}}{\mathop{\mathop{\sum }^{}}}\,\hat{b}\,_{n}^{\dagger }{{\hat{b}\,}_{{{n}^{\prime }}}} \\ {} & - & \frac{2\pi D{{\Omega }_{\text{eff}}}}{3}\underset{n\ne {{n}^{\prime }}}{\mathop{\mathop{\sum }^{}}}\,\hat{b}\,_{n}^{\dagger }{{\hat{b}\,}_{n}}\hat{b}\,_{{{n}^{\prime }}}^{\dagger }{{\hat{b}\,}_{{{n}^{\prime }}}}, \\ \end{eqnarray} \tag{ A.10 }$$

where the average dispersion energy of a pair in Bloch band n is given by ${{\epsilon }_{n,\text{avg}}}=2\mathop{\int }^{}{{E}_{n}}\left( q \right)dq$. From the Hamiltonian equation (A.10), by taking a mean-field type approximation for the composite operator will again result is the effective shift in the dispersion.

Let us reconsider the model of a pair distributed over neighbouring sites. This time we include the effect of pair breaking due to the single particle tunneling matrix in Hamiltonian equations (A.4) and (A.6). To do that, within the two-site model, we have reevaluated the dynamics of the pairs by taking into account states with single molecule per site. Our initial state consists of the situation where only one of the site contains a pair.

With this initial condition, we have carried out the full dynamics within the two-site case and the result is presented in figure 5. There we have plotted the total population of the single-particle states. We see that the maximum population of the single particle states are less than $<0.1$. The main reason for such anobservation is that within the Wannier orbitals, the effective single-particle tunneling terms are much smaller (due to the aspect ratio of a site in the lattice) than the continuum induced pair tunnelings that are independent of any local aspects of the Wannier function. This justifies our assumption of neglecting the pair-breaking effect of the single-particle tunneling Hamiltonian. Moreover, due to such a negligible population of the single-particle states, the decay rates of various channel remains unchanged with respect to the case discussed in the paper.

Zoom In Zoom Out Reset image size

We discuss here the effect of rotational level mixing due to quantum nature of the dipolar interaction, the effect neglected in the main text. Such a mixing gives rise to an effective vdW like potential which decays with distance r as $-1/{{r}^{6}}$ [49]. To look into its effect, we first consider a polar molecule with dipole moment μ, rotational constant ${{B}_{e}}$ is polarized by a strong electric field $\mathbf{E}$ along the z direction. In the limit of $\left( \mu \mathbf{E}/\hbar {{B}_{e}} \right)1$, one can write the rotational Hamiltonian in the M = 0 sector (M is the projection of angular momentum along the molecular axis) as,

$$\begin{eqnarray*}&&{{H}_{\text{rot}}}=\hbar {{B}_{e}}{{\hat{J}\,}^{2}}-\mu \mathbf{E}\cos \theta \approx -\hbar {{B}_{e}}\partial _{\theta }^{2}+\mu \mathbf{E}{{\theta }^{2}}/2,\end{eqnarray*}$$

where θ is the angle between the molecular axis and the electric field direction and μ is the permamnent dipole moment. The energy levels of the Hamiltonain ${{H}_{\text{rot}}}$ is denoted by the index $m=0,1,2,.....$ with energy ${{E}_{\text{rot},\text{m}}}=\left( 2\;m+1 \right)\hbar {{B}_{e}}/d_{\theta }^{2}$ and wavefunction ${{\Phi }_{m}}\left( \theta \right)={{N}_{m}}\exp \left( -{{\theta }^{2}}/2d_{\theta }^{2} \right){{H}_{m}}\left( \theta /{{d}_{\theta }} \right)$ where ${{H}_{m}}\left( . \right)$ is the Hermite polynomial of order m, ${{N}_{m}}$ is the normalization constant and the width ${{d}_{\theta }}={{\left[ 2\hbar {{B}_{e}}/\mu \mathbf{E} \right]}^{1/4}}$. The rotational state of the polar molecule is denoted by the lowest energy rotational wavefuntion ${{\Phi }_{0}}\left( \theta \right)$ which induced a dipole moment of ${{\mu }_{\text{ind}}}\mu \mathop{\int }^{}\cos \theta \Phi _{0}^{2}\left( \theta \right)d\theta$. This results in dipolar interaction between the ground state molecules which falls of as $1/{{r}^{3}}$. Additionally, dipolar interaction also induces excitations to higher energy rotational states. Within second order perturbation theory, the resulting effective interaction between the ground state molecules, in the units of recoil energy, is given by,

$$\begin{eqnarray}&&{{V}_{\text{vdW}}}\left( \mathbf{r} \right)\approx -\frac{\ell _{\text{vdW}}^{4}}{{{r}^{6}}}{{\left[ 1-3{{z}^{2}}/{{r}^{2}} \right]}^{2}},\end{eqnarray} \tag{ B.1 }$$

where the distance are in the units of $a/\pi$ and the effective dimensionless vdW length ${{\ell }_{\text{vdW}}}$ is given by,

$$\begin{eqnarray*}&&{{\ell }_{\text{vdW}}}={{\left[ d_{\theta }^{2}\exp \left( -d_{\theta }^{2}/2 \right)\frac{D_{\max }^{2}{{E}_{R}}}{4\sqrt{2}\hbar {{B}_{e}}} \right]}^{1/4}}/\pi ,\end{eqnarray*}$$

where the maximum dipolar strength is given by ${{D}_{\max }}={{m}_{b}}{{\mu }^{2}}/2{{\epsilon }_{0}}{{\hbar }^{2}}a$. For RbCs molecule, the rotational constant is given by ${{B}_{e}}=0.014$ cm$^{-1}$and the permanent dipole moment is given by $\mu =1.27$ Debye. Then for an electric field strength of E = 10 kV/cm, the angular width reads ${{d}_{\theta }}=0.7$. Correspondingly, the vdW length is given by ${{\ell }_{\text{vdW}}}\approx 0.17$ when the lattice constant is a = 500 nm. From this we can also define a short distance cutoff scale ${{\ell }_{\text{sr}}}$ where rotational mixing effect of the dipoles becomes similar magnitude to the rotational splitting [49]. In our units, this cut off is given by ${{\ell }_{\text{sr}}}\approx .03$ for dipolar strength D = 20. For length scales $r>{{\ell }_{\text{sr}}}$, the perturbative form of the vdW interaction in equation (B.1) remain valid and for $r<{{\ell }_{\text{sr}}}$, the rotational level of the molecules becomes strongly mixed and the deeply bound molecular pairs appears [37].

Following the discussion in the main text and the above sections, we write the vdW Hamiltonian in the discrete (${{H}_{\text{vdW},\operatorname{int}}}$), continuous (${{H}_{\text{vdW},\text{Bloch}}}$) and discrete-continuous (${{H}_{\text{vdW},\text{WB}}}$) sector. The interaction in discrete sector is weak compare to the dipolar interaction. This can be easily seen by Fourier transforming equation (B.1), ${{V}_{\text{vdW}}}\left( k \right)\approx \ell _{\text{vdW}}^{4}{{k}^{3}}\mathcal{F}\left( k{{\ell }_{\text{vdW}}} \right)$, with the function $\mathcal{F}\sim 1$. The widths of the Wannier functions in the momentum space are of the order of $k\sim 1$. Then as $\ell _{\text{vdW}}^{4}\sim {{10}^{-3}}1$, we can neglect the vdW interaction in the discrete states compare to the dipolar strength in equation (4).

Moreover, one can estimate the loss rate due to the coupling of the bound molecular complex by evaluating the overlap between the Wannier orbitals and the bound state wave function which is of the order of $\exp \left( -1/{{\ell }_{\text{vdW}}} \right)\rho$ where ρ is the density of bound states in the units of recoil energy. Here we have assumed that the bound state decays exponentially for a large distance. From [37], for RbCs molecue in the rotational ground state, the density of states is large, $\rho \sim 40$. Accordingly, the decay rate will be proportional to the overlap which is of the order of $0.1{{E}_{R}}$ which gives a timescale of $\sim 1.0$ ms. For other species of molecules it is possible that the density of bound states is lower which can result in an increased stability.

In the continuous Bloch band, the corresponding momentum scale is given by $k\sim \sqrt{{{E}_{n}}}$ and the corresponding strength of the vdW interaction in the continuum band n is in the order of $-\ell _{\text{vdW}}^{4}E_{n}^{3/2}$. Whereas from equation (A.9), we find that the strength of the dipolar attraction in this band is of the order of $\sim D$. For dipolar strength of $D\sim 20$, the vdW interaction gets prominent only for very high Bloch bands with ${{n}_{\text{cutoff}}}\gtrsim 1/{{\ell }_{\text{sr}}}\approx 30$. As such bands probes distance shorter than ${{\ell }_{\text{sr}}}$, this will result in strong overlap (or in other words, strong coupling) with the molecular bound states which can give rise to phenomenon of molecular sticking [37]. Subsequently, any population in those Bloch bands will result in loss due to formation of strongly bound molecular pairs and we denote this by loss rate ${{\Gamma }_{\text{vdW}}}$.

The situation remains similar also for the discrete to continuous transitions. There the transitions happens between states with momentum $k\sim {{\sqrt{E}}_{n}}$ continuous states and discrete states with momentum $k\sim 1$. As for continuous states, ${{E}_{n}}\gg 1$, the corresponding vdW discrete-continuous transition strength is of the order of $-\ell _{\text{vdW}}^{4}E_{n}^{3/2}$. Subsequently, for intermediate momentum, the discrete-continuum transition is again dominated by the dipolar terms in Hamiltonian equation (5) and molecular sticking due to vdW interaction involves very high energy Bloch bands (or shorter distance) with band index $n\gtrsim 30$.

Accordingly, while integrating out the high Bloch bands, one get additional terms (equivalent to the term in equation (B.2)),

$$\begin{eqnarray} &&\underset{n>{{n}_{\text{cut}}}}{\mathop{\mathop{\sum }^{}}}\,\quad \left[ \mathbf{P}_{ij,{{q}_{1}}{{q}_{2}}}^{n} \right]\frac{{{\left[ \mathbf{P}_{ij,{{q}_{1}}{{q}_{2}}}^{n} \right]}^{\dagger }}\mathbf{C}}{s-\mathbf{i}\left[ {{E}_{n}}\left( {{q}_{1}} \right)+{{E}_{n}}\left( {{q}_{2}} \right)-\pi D{{\Omega }_{\text{eff}}} \right]+{{\Gamma }_{\text{vdW}}}} \\ &&\approx {{\tan }^{-1}}\left[ \frac{{{\Gamma }_{\text{vdW}}}}{n_{\text{cut}}^{2}-\pi D{{\Omega }_{\text{eff}}}} \right]\approx 0, \\ \end{eqnarray} \tag{ B.2 }$$

as $n_{\text{cut}}^{2}\gg D$, and the decay rate ${{\Gamma }_{\text{vdW}}}\sim D$ for a lattice constant of 500 nm and dipolar strength $D\sim 20$. We have calculated ${{\Gamma }_{\text{vdW}}}$ from [37] but assuming temperature in the nano-Kelvin regime which suppresses the d-wave resonances.

We rewrite our Hamiltonian equation (18) in the conventional ${{J}_{1}}-{{J}_{2}}$ form as,

$$\begin{eqnarray} {{H}_{\text{eff}}} & = & {{J}_{1}}\underset{\langle ij\rangle }{\mathop{\mathop{\sum }^{}}}\,b_{i}^{\dagger }{{b}_{j}}+{{J}_{2}}\underset{\langle \langle ij\rangle \rangle }{\mathop{\mathop{\sum }^{}}}\,b_{i}^{\dagger }{{b}_{j}} \\ {} & + & V\underset{ij}{\mathop{\mathop{\sum }^{}}}\,\frac{{{n}_{i}}{{n}_{j}}}{|i-j{{|}^{3}}}-\mu \underset{i}{\mathop{\mathop{\sum }^{}}}\,{{n}_{i}}, \\ \end{eqnarray} \tag{ B.3 }$$

where ${{J}_{1}},{{J}_{2}}$ are the nearest and next-nearest neighbour tunneling and V is the strength of the long-range interaction. In the dilute limit, such a system, with nearest and next-nearest neighbour interaction only, has been solved qualitatively by mapping the problem to a two-component Bose gas model [39]. Here we extend this treatment to include long-range dipolar interaction. To do that we transform the Hamiltonian to the momentum space,

$$\begin{eqnarray}&&{{H}_{\text{eff}}}=\underset{q}{\mathop{\mathop{\sum }^{}}}\,{{\epsilon }_{q}}b_{q}^{\dagger }{{b}_{q}}+\underset{k,{{k}^{\prime }},q}{\mathop{\mathop{\sum }^{}}}\,V\left( q \right)b_{k+q}^{\dagger }b_{{{k}^{\prime }}-q}^{\dagger }{{b}_{{{k}^{\prime }}}}{{b}_{k}}-\mu \underset{i}{\mathop{\mathop{\sum }^{}}}\,b_{q}^{\dagger }{{b}_{q}},\end{eqnarray} \tag{ B.4 }$$

where the dispersion relation is given by ${{\epsilon }_{q}}=2{{J}_{1}}\cos qa+2{{J}_{2}}\cos 2qa$ and the interaction energy in momentum space is given by, $V\left( q \right)=U+2\;V\mathop{\sum }^{}_{n=1}^{\infty }\cos nqa/{{n}^{3}}$, where the hard-core constraint is given by $U\to \infty$. We only consider the dilute limit, $\mu \to 0$. When ${{J}_{2}}>{{J}_{1}}/4$, the dispersion relation has two minima at wavevectors, $Qa={{\cos }^{-1}}\left[ -{{J}_{1}}/4{{J}_{2}} \right]$. Around these minima, we can write the dispersion relation as, ${{\epsilon }_{Q+k}}={{\epsilon }_{Q}}+{{\hbar }^{2}}{{k}^{2}}/2\;{{m}^{*}}$, where ${{m}^{*}}$ is the effective mass. Then we expand the boson operator near the two minima, ${{b}_{k}}={{\phi }_{1,Q+k}}+{{\phi }_{2,-Q+k}}+{{\phi }_{k}}$, where ${{\phi }_{1}}$ and ${{\phi }_{2}}$ are the two-component Bose gas centered around momentum $\pm Q$ respectively, while ${{\phi }_{k}}$ denotes the high momentum contribution, which is integrated out. Then one can re-express the Hamiltonian equation (B.4) in terms of the ${{\phi }_{1,2}}$ which in position space reads,

$$\begin{eqnarray*} {{H}_{\text{eff}}} & = & \mathop{\int }^{}dx[\underset{\sigma =1,2}{\mathop{\mathop{\sum }^{}}}\,\left[ -\phi _{\sigma }^{\dagger }\frac{{{\hbar }^{2}}}{2{{m}^{*}}}\nabla _{x}^{2} \right]{{\phi }_{\sigma }} \\ {} & {} & +\frac{1}{2}{{T}_{1}}(\rho _{1}^{2}+\rho _{2}^{2})+{{T}_{12}}{{\rho }_{1}}{{\rho }_{2}}-\mu ({{\rho }_{1}}+{{\rho }_{2}})], \\ \end{eqnarray*}$$

where ${{T}_{1}}$ and ${{T}_{12}}$ are renormalized interactions. To find these renormalized interactions, we first write down the full Bethe−Salpeter equation,

$$\begin{eqnarray}&&T\left( k,k\prime ;q \right)=V\left( q \right)-\mathop{\int }^{}\frac{V\left( p-q \right)T\left( k,k\prime ;p \right)}{{{\epsilon }_{k+p}}+{{\epsilon }_{{{k}^{\prime }}-p}}+\Omega }\frac{dp}{2\pi }.\end{eqnarray} \tag{ B.5 }$$

In the dilute limit we can substitute $\Omega =2\mu$. Then the respective renormalized interaction is given by, ${{T}_{1}}=T\left( Q,Q,0 \right)$ and ${{T}_{12}}=T\left( Q,-Q;0 \right)+T\left( Q,-Q;2Q \right)$. Imposing the hard-core constraint with $U\to \infty$, we get an additional equation,

$$\begin{eqnarray*}&&\mathop{\int }^{}\frac{T\left( k,k\prime ;p \right)}{{{\epsilon }_{k+p}}+{{\epsilon }_{{{k}^{\prime }}-p}}+\Omega }\frac{dp}{2\pi }=1.\end{eqnarray*}$$

Due to the form of the interaction V(q), we expand the full renormalized interaction as, $T\left( k,k\prime ;q \right)={{A}_{0}}+{{\mathop{\sum }^{}}_{n}}{{A}_{n}}\cos nqa$, where the coefficients ${{A}_{0}},{{A}_{n}}$ depends on $k,k\prime$. Putting this ansatz in equation (B.5), we get a set of coupled equations for $m>0$,

$$\begin{eqnarray*}&&{{A}_{m}}=\frac{2\;V}{{{m}^{3}}}-\frac{2\;V}{{{m}^{3}}}\underset{m\prime =0}{\overset{\infty}{\mathop{\mathop{\sum }^{}}}}\,{{A}_{{{m}^{\prime }}}}\mathop{\int }^{}\frac{\cos mpa\cos m\prime pa}{{{\hbar }^{2}}{{p}^{2}}/{{m}^{*}}+\Omega }\frac{dp}{2\pi },\end{eqnarray*}$$

and from the constraint condition,

$$\begin{eqnarray*}&&{{A}_{0}}=\underset{m=0}{\overset{\infty}{\mathop{\mathop{\sum }^{}}}}\,\mathop{\int }^{}\frac{{{A}_{m}}\cos mpa}{{{\hbar }^{2}}{{p}^{2}}/{{m}^{*}}+\Omega }\frac{dp}{2\pi }.\end{eqnarray*}$$

We found that in the limit of $\Omega \to 0$, the magnitude of the integral like $\mathop{\int }^{}\frac{\cos mpa}{{{\hbar }^{2}}{{p}^{2}}/{{m}^{*}}+\Omega }\frac{dp}{2\pi }$ falls off when $m>1$. Then we get the following relation,

$$\begin{eqnarray*}&&{{A}_{m}}=\frac{2\;V}{{{m}^{3}}}-\frac{2\;V}{{{m}^{3}}}\underset{{{m}^{\prime }}=m-1}{\overset{m+1}{\mathop{\mathop{\sum }^{}}}}\,{{A}_{{{m}^{\prime }}}}\mathop{\int }^{}\frac{\cos \left( m-m\prime \right)pa}{{{\hbar }^{2}}{{p}^{2}}/{{m}^{*}}+\Omega }\frac{dp}{2\pi }.\end{eqnarray*}$$

We numerically find convergent solution for the ${{A}_{m}}$ by taking ${{m}_{\max }}=100$.