Many body population trapping in ultracold dipolar gases

A system of interacting dipoles is of paramount importance for understanding of many-body physics. The interaction between dipoles is {\it anisotropic} and {\it long-range}. While the former allows to observe rich effects due to different geometries of the system, long-range ($1/r^3$) interactions lead to strong correlations between dipoles and frustration. In effect, interacting dipoles in a lattice form a paradigmatic system with strong correlations and exotic properties with possible applications in quantum information technologies, and as quantum simulators of condensed matter physics, material science, etc. Notably, such a system is extremely difficult to model due to a proliferation of interaction induced multi-band excitations for sufficiently strong dipole-dipole interactions. In this article we develop a consistent theoretical model of interacting polar molecules in a lattice by applying the concepts and ideas of ionization theory which allows us to include highly excited Bloch bands. Additionally, by involving concepts from quantum optics (population trapping), we show that one can induce frustration and engineer exotic states, such as Majumdar-Ghosh state, or vector-chiral states in such a system.


Introduction
In recent years the ultracold gases have been used as a tool to quantum engineer various novel states of matter with an unprecedented precision and control. In this regard, particularly challenging is the engineering of frustrated systems for ultracold gases trapped in optical lattices. Frustration can either be induced by the lattice geometry, which can lead to kinetic frustration, or by higher order exchange processes due to strong interactions [1,2]. Polar molecules are particularly interesting in this context, as they can interact via long-range dipolar forces, which can induce yet another kind of frustration. In particular, dipolar lattice gases have been proposed to simulate various quantum phases and exotic phenomena, such as supersolidity [4,5], quantum magnetism [6], topological states [7,8], exotic pair-superfluidity [9], etc. Experimental progress towards creation of quantum degenerate gas of ground state polar molecules has been spectacular over the last years [10,11,12,13], leading, for instance, to realization of quantum spin models using fermionic molecules [14] or dipolar Chromium atoms [15].
One of the important properties of the polar molecules is that their dipole moment can be tuned by applying an electric field. The more polarized these molecules get, the stronger becomes the dipolar interaction between them. Theoretically it is a challenge to investigate the properties of these strongly interacting molecules trapped in an optical lattice. The standard approach based on Bose-Hubbard models limited to the lowest Bloch band [4,5,6,7,8,9] becomes inapplicable due to strong interaction induced coupling between the bands. In this paper we provide a novel route to describe such strongly interacting systems. Specifically, we consider bosonic polar molecules trapped in a one dimensional optical lattice. We find that the system can be modeled with effective couplings between the localized states at lattices sites and the continuum of highly excited states. This connects our approach to the extensive studies of strong laser field induced ionization of atoms and molecules. In particular we find analogies to auto-ionization processes, in which multi-configuration interactions couple discrete states with continua, as in the celebrated Fano model [16]. Usually, due to the coupling to the continuum, the electrons in atoms or molecules are transferred from the bound states to the continuum, which leads in the long-time limit to the irreversible decay of bound state population. Strong laser field, however, enables efficient couplings between different ionization paths leading to various interference phenomena. For strong field auto-ionization it may lead to the so called confluence of coherence [17,18], which slows down very efficiently the ionization process.
Similarly, if several (at least two) bound states are coupled to a common continuum, a phenomenon of coherent population trapping may occur -the ionization is incomplete and a significant part of the system population is trapped in the bound subspace [19]. The resulting stable bound configuration is a superposition of original bound states with properties depending on the details of the coupling to the continua. The coherent population trapping phenomenon appears also for multi-level discrete systems when coherent driving may create non-absorbing states (often called "dark states") -for a review of coherent population trapping see [20]. Most importantly, in our system of polar molecules, we find that similar phenomenon can give rise to frustration in lattice systems, as the population trapping can involve particles trapped in different sites of the optical lattice. Specifically, we find that for a half-filling, the many-body population trapped state is a dimer state known in the condensed-matter physics as Majumdar-Ghosh state [21]. Majumdar-Ghosh state is a paradigmatic example in the study of frustrated models, since it retains basic properties of spin-liquid phases, such as fractional excitations [22]. For lower filling we find that the effective model can be written as a J 1 − J 2 Hamiltonian with nearest and next-nearest neighbour tunneling, along with the long range dipolar interactions. Similar models, restricted only to nearest and next-nearest neighbor tunneling, have been investigated for long in connection with various magnetic materials [22]. But, in solid-state materials [23] as well as in optical lattices [2], such next-nearest neighbour tunneling can only come from higher order exchange processes, which makes it considerably weaker than the nearest neighbour tunneling. The corresponding temperature is thus very low. Amazingly, the temperature scale associated with population trapped frustration remains comparable to the characteristic temperature scale of the system. An alternative way to achieve long range "tunneling" in spin models is offered in ultracold ions setting [24,25,26], but such systems are not easily scalable to macro-or even meso-scopic sizes.

The model
We consider bosonic polar molecules trapped in an optical potential inducing a onedimensional lattice geometry, where V 0 denotes the lattice depth and a is the lattice constant. Ω denotes the harmonic (strong) trapping frequency along the y and z direction. The molecules are polarized by an electric field along the z axis. To describe this system we make two assumptions: i) along the trapping directions, only the lowest harmonic oscillator eigenstate is occupied, and ii) at time t < 0, repulsively bound pairs of molecules in the limit of weak dipolar strength are prepared by tuning the lattice depth [27], or by applying a weak electric field. The molecules of the pair repel each other and cannot separate due to the energy conservation -a separation would imply populating single particle states in the band gap. Then at t = 0, we switch on a strong polarizing electric field to induce a strong dipolar interaction between the molecules. The strength of the dipolar interaction in dimensionless units is denoted by D = m b µ 2 ind /2 0 2 a, where µ ind is the effective dipole moment controlled by the external electric field, 0 is the vacuum permittivity and m b is the mass of the molecules. The dipolar interaction is given by To describe the strongly interacting regime for our system, one needs to go beyond the simple single band tight-binding approximations [28]. The single-particle motion in a periodic potential results in energy bands, known as Bloch bands which can be expressed in terms of quasi-momentum q. For each such a band, one constructs localized basis states or orbitals [the so called Wannier functions(WF)] from the Bloch states [29]. By taking into account only the lowest energy Wannier states, one arrives at a Hamiltonian containing density-density interactions terms (both on-site and long-range) and nearest neighbour tunneling processes along the x direction. In the presence of strong interactions such an approximation breaks down due to two primary reasons: i) the interaction mixes different bands or orbitals and different sites (specially for the higher orbitals), and ii) for higher orbitals, one has to take into account long-range tunneling matrix elements. These problems have been partially addressed taking into account higher excited bands in the tight-binding approximation and considering only the onsite interactions [30,31,32,33,34,35,36]. For strong interactions though, serious complications appear due to the lack of convergence of results as a function of number of bands taken into consideration. Subsequently, standard approaches become questionable and impractical. The effective description that interaction effectively increase or decrease the width of the Wannier functions [33] does not hold. For such strong interactions an entirely new approach is needed. The method initiated in this work paves the way for efficient description of such systems.
The essential observation, forming the core of our approach, is that for typical optical lattice depths, only few lowest bands are separated from each other energetically with forbidden gaps in between. The higher bands, in reality, form a continuum of energies. The simplest situation occurs for relatively small lattice depths say of few recoil energies, as shown in Fig. 1(a) for V 0 = 7E R . Here two lowest, s and p, bands are separated from the continuum formed by other bands. The Wannier states (called often orbitals) of the first two bands are relatively well localized, and the mean energy calculated for them is lower than the optical lattice depth V 0 . In this situation it is natural to express the motion of the particles in a mixed basis, where only the lowenergy motion is expressed in terms of localized Wannier orbitals. Instead of using Wannier basis for the other bands too, as in the standard approaches, the remaining higher energy states will be treated by continuous Bloch functions. For much deeper lattices a natural generalization of our approach will be to take more than two discrete bands into account; in this work we limit ourselves to the simplest situation. Therefore, we may write down the field operator in the chosen mixed Wannier-Bloch basis as: where ω α i (x) is the localized Wannier function at site i corresponding to α = s or α = p orbital while φ 0 is the lowest harmonic oscillator eigenfunction for trapping directions.ŝ † i ,ŝ i ,p † i ,p i are the creation and annihilation operators for the bosons in the s-and p-orbitals. U n (q) denotes the Bloch functions for band n with quasi-momentum q (n > n 0 = 2 the latter counts the number of bands treated using Wannier basis). Consequently,â † nq ,â nq denote the boson creation and annihilation operators in the bands considered in the Bloch basis with a quasi-momentum q. To define dimensionless quantities, we first rescale the distance πx/a → x, that defines the scale for the energy E R = π 2 2 /2m b a 2 . So in the limit of E n 0 V 0 ,the simplified Hamiltonian in the Wannier-Bloch basis is given by, with and where E p gives the single particle energy of the p-orbital (E s = 0 is assumed). Note that, from the very begining we omit single particle tunneling terms between sites despite the lattice depth being low. That assumption is due to the fact that we shall consider a specific preparation of the system (see below) in form of pairs. The tunneling of pairs can be possible due to second order processes only. The single particle tunnelings, on the other hand, are reduced for dipoles by interaction mediated density-dependent (bond-charge) tunneling terms as discussed in [9]. The full Hamiltonian is given in the Appendix A while the estimates of the effects due to single-particle and correlated tunneling terms are discussed in Appendix B.
The interaction between particles in localized orbitals H Int contains (with σ, σ = s, p denoting the orbitals) the onsite intra-orbital interactions U σσ , the possible transitions of a pair between orbitals with the strength T ps and the long range dipolar interaction (again, additional terms in the Hamiltonian have negligible effect as discussed in the Appendix B). H Bloch in Eq.(3) contains the kinetic energy of the molecules in the continuous band, q,n>n 0 E n (q)â † nqâ nq as well as interaction between particles in the continuum (see Methods section). The Wannier-Bloch Hamiltonian part H WB describes the coupling between Wannier-described sites with two particles and the Bloch continuum. P n i,ss (q 1 q 2 ), P n i,sp (q 1 q 2 ) and P n i,pp (q 1 q 2 ) are the corresponding coupling constants of two particles at site i and the continuum for the ss-, sp-and pp-orbitals respectively. A cartoon of these various transition processes is shown in Fig.1(b). We would like to stress that in Eq. (4) we have taken into account the contribution from the dipolar interactions only. There are additional Van-der Waals terms arising from the mixing of rovibrational levels of molecules. Such contributions can potentially lead to a formation of long-lived molecular complexes as described in Ref. [37] for RbCs molecules resulting in additional loss processes which will limit the density of molecules in a lattice. Though for molecules with low density of bound molecule-molecule states such loss rate can be considerably lower.
To simplify the notation we denote the basis states for zero or two particles on a site as where the state |n 1 n 2 denotes n 1 particles in the s-orbital and n 2 particles at p-orbital.
We refer to these states as Wannier states in the following sections. Before considering the physics generated by the postulated Hamiltonian let us mention also that we treat the molecules rather brutally, considering them as simple dipoles. In particular we neglect the rotational structure of molecular energy levels and the induced rotational level mixing (with the effective van der Waals potential) [50]. In unfavorable situations that may lead to creation of deeply bound molecular pairs [37] whose large kinetic energies allows them to leave the optical lattice potential resulting in a strong loss. These effects are discussed in more detail in Appendix B, we believe that in the parameters regime discussed below these effects can be neglected.

Interesting configurations
A large variety of different situations may be considered for the model studied. Let us imagine the situation when the system is prepared (for typical weak interactions) in an insulating state, for example the Mott state. We assume that at t = 0 we suddenly switch on the electric field which strongly polarizes the molecules inducing large dipoles along the static field direction (assumed perpendicular to the lattice). The interaction between dipoles becomes strong making the analysis of the system difficult.
Whether strong interactions will destabilize the system if the interaction energy becomes comparable to binding in the lattice? May be some metastable states still survive leading to interesting effects? These are the basic questions we want to address.

A single pair of molecules in neighboring sites
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 W B , 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, where |l i denotes the state of the system at site i [following the notation of (6)]; |0 denotes the vacuum for the continuum and |1 q 1 1 q 2 denotes the state with both particles in the continuum corresponding to the quantum numbers n 1 q 1 and n 2 q 2 . The timedependent Schrödinger equation for |Φ leads to a set of coupled equations for probability amplitudes C l , grouped in a 6-component vector C, corresponding to discrete states, as well as for continuum amplitudes α n 1 n 2 q 1 q 2 .
where |i−j| = 1 and U 1 is the interaction matrix between the discrete states originating from the Hamiltonian (4): Due to a lack of the direct coupling between the Wannier states at different sites, U 1 is block diagonal. The Bloch-Wannier Hamiltonian in Eq.(5) will give rise to the discretecontinuum coupling array P n ij,q 1 q 2 = [P n i,q 1 q 2 , P n j,q 1 q 2 ] T . To find the time evolution of the pair in the continuum, we make the ansatz that α nn q 1 q 2 ≈ α 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 Eq.(8) denotes coupling of population amplitude of a Bloch band n to that of another Bloch band n . The corresponding coupling strength ∼ 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 Eq. (8), n − n = ±1 and α n ≈ α n−1 ≈ α n+1 . Within these approximations, one can rewrite Eq. (8) as, 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 Eq.(11) we get, Then one can do similar Laplace transform for the discrete state amplitudes in Eq. (8) and eliminate the continuum amplitudes by Eq. (12). Subsequently, in the time evolution of the discrete state amplitudes, one gets expressions like, where we have introduced a cut off n cut ∼ 20 in the band index and i = √ −1. Any excitations to higher bands than n 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 − πDΩ eff /2 0. Then by transforming the summation over energy level n to integration, one integrates over the range −∞ → ∞.
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Ċ = MC. The coupling matrix M is expressed as where we have introduced the decay matrix Γ and the effective trapping strength Ω eff = Ω/2E R . In the expression above P n q 1 q 2 is a vector of couplings of 6 Wannier discrete states [3 per site -compare (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 Eq. (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 where u l is the eigenvector of the matrix M with Γ l and l being the decay rate and the energy of the l-th eigenstate for the neighbouring sites. In Fig. 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 (the black line and the blackcircled 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 |± ij = (|φ i ± |φ ) / √ 2 with energies ± with + < − : expressed in terms of s and p orbitals. The overlap of these approximate combinations with the exact eigenstates: | σ| σ exact | i j ≈ Fδ σσ is large with F ∼ 0.95, where σ, σ = ±. 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, The state with the lowest decay rate [the black line in Fig.2 (left panel)] corresponds to the state |− ij with highest energy. For this state we find that the ratio between the decay rate and the energy lies in the range, Γ − / − = 0.01 → .05 as the dipolar strength changes from 10 → 50. On the other hand, for the state |+ ij , for the same dipolar range, Γ + / + = 0.05 → 0.1.
It follows that on the timescale of ∼ 1/Γ + ≈ 10/E R , only the |− 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 Fig.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.

Continuum-assisted creation of dimer states
Next we discuss the creation of dimer states due to population trapping for the halffilling of the pairs. It is known that the strong dipolar interaction induces a density-wave phase where the pairs arrange in a checkerboard pattern [41]. 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 |− ij . Then, in the limit of strong interaction and for half-filling of pairs, |− ij will cover the whole region of lattice sites. The resulting manybody state is a checkerboard state of nearest-neighbour dimers, . 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, b † ib j Ψ = n inj = 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 δE between the single-site state |φ i and the dimer state |− ij . For a dipolar strength of D ∼ 20 (near the lowest decay rate in Fig.2) this energy difference is of the order of 0.4E R . For RbCs molecules, these parameters correspond to a dipole moment of ∼ 0.7Debye with a lattice constant ∼ 500nm. Then the relevant temperature scale to observe this phase is ∼ 50nK. 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 ∼ 10ms. 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.

Many-body effects due to long-range dipolar interaction
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 C mb with modified continuum induced transition matrix P mb where we have taken into account continuum induced long-range coupling. The resulting equation of motion has the form, C mb = M mb C mb , and the many body coupling matrix M mb is given by where the discrete states interaction matrix U mb now also includes the long-range dipolar interaction. We then find the eigenvalues and eigenstates of the matrix M 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 (|Ψ A , |Ψ 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 |Ψ − = [|Ψ A − |Ψ B ] / √ 2. We denote this overlap by the function F M G and plot it against the dipolar strength in Fig. 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.

Constructing effective Hamiltonian for low filling
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 M(i, i + 1, i + 2) 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 |± ij states. Subsequently, a tunneling Hamiltonian in terms of the states |± ij is given by, , where α can be extracted from the eigenvalues of the effective coupling matrix M(i, i + 1, i + 2). As the states |− i,i+1 , |± 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, |φ i = b † i |0 which creates a pair at site i in the lowest decay state. The pair operators satisfy bosonic commutation relations In terms of these pair operators we can rewrite the states as |± ij = 1 . Subsequently, the Hamiltonian H i,i+1,i+2 is re-expressed as, The values of J eff and α are derived by comparing the energies of Hamiltonian (17) and the energies of the states with three lowest decay rates derived from the full coupling matrix M(i, i + 1, i + 2). For small values of the dipolar strength D we find that α ≈ 1, thus the long-range tunneling is small and one recovers the usual picture with nearestneighbour tunneling only. But for higher dipolar strengths α = 1, due to the different decay rates of |± states. For such values of α one obtains, therefore, [compare (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 longrange dipolar interaction and involving all sites as, where we have introduced the chemical potential µ for the pairs and << ij >> is a shorthand for next nearest neighbour summation index. The Hamiltonian in Eq.(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, (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 q = J eff (1 + α) cos qa + J eff 1−α 2 cos 2qa. For 1−α 1+α > 1/2 it shows two minima at wavevectors ±Qa = cos −1 − 1+α 2(1−α) . 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 [40] with the effective Hamiltonian, where ρ 1,2 are the densities of the two component Bose gas centered around the the minima ±Q and T 1 , T 12 are the renormalized intra-component and inter-component interaction. A detailed discussion of the Hamiltonian (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 [40]. When T 1 < T 12 , the mean-field ground state solution is given by the phase-separated state ρ 1 = 0, ρ 2 = 0 or ρ 1 = 0, ρ 2 = 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, ψ s = √ ρ s exp[−iθ s ], with s = 1, 2. When ρ 1 = 0, we Such a "cone" phase is identified as a the vectorchiral (VC) phase which breaks the 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 ρ 1 = ρ 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 ∼ 2Q. Such bound states with finite momenta are usually not present in the anti-ferromagnetic model [40]. In the present case, these bound states are a direct consequence of the long-range nature of the dipolar interaction which can induce resonances [42]. The quasi-condensate of such bound pairs can give rise to a spin-nematic phase [43,44], or spin-density wave phase [45], a detailed discussion of which is beyond the scope of current article. In Fig. 4, we have plotted the resulting phase diagram in the D − µ 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.

Discussion
Summarizing, in the present article we have demonstrated a novel approach to the problem of strongly interacting molecules in optical lattices. We have explored a mathematical analogy between the system studied and strong bound-continuum couplings present in the theory of strong field ionization. We have found that the phenomenon of coherent population trapping, a well known interference effect in quantum optics, is responsible for frustration in our system in a form of dimerization and next-nearest neighbour tunneling. One strong point of our proposal is that the required temperature scale is much higher than the one corresponding to the usual super-exchange regime. Our results can be generalized to higher dimensions, where one can look for simulation of spin liquids, and valance bond crystals [22]. Our method can also be extended to other strongly interacting systems, such as atoms in optical lattices, strongly-coupled cavity-QED systems [46], recently proposed nanoplasmonic lattices [47], and possible lattice geometries for the indirect excitons with strong dipolar interactions [48]. We hope that further progress can be obtained in studies of strongly interacting systems by exploring analogies with strongly coupled quantum optics problems in general, and strong field ionization theory in particular.
The different amplitudes in the discrete subspace Hamiltonian are obtained by appropriate integrals of the dipole-dipole interaction potential and the mode functions [compare (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 E n=n 0 +m (q) = E n 0 + 2n 0 m + m 2 + 2(n 0 + m)|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 1 limit as [49], As these functions are eigenstates, they are also orthogonal, U * nq (x)U * mq (x)dx = δ n,m δ q,q . Now we write down the discrete-continuum coupling matrix elements as, is the Fourier transform of the Wannier function ω σ i (x). In deriving the above form, we have used the orthogonality condition between the Bloch functions and assumed that E n 1. Additionally, in the Hamiltonian (3), we have neglected terms corresponding to processes likeâ † nq 1ŝ † iŝ iŝi where one particle is coupled to the continuum. The transition amplitudes for such processes contains convolution sums of the form As E n 1, such terms are negligibly small. Thus we ignored them in comparison to the leading two-particle transition amplitudes.

Hamiltonian for the continuum states
The Hamiltonian for the continuum Bloch states reads the continuous band index n, n > n 0 and the momentum index q = [q 1 , q 2 , q 3 , q 4 ]. The second term in the Hamiltonian (28) 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 ∼ D. Furthermore, from Hamiltonian (28), we notice that the interaction is strongly attractive in the higher Bloch bands and for strong interaction (D 1), the dipolar strength can exceed the width of the first few continuous Bloch bands.

Many-body effects in the continuum
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, 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 (28) in terms of the composite operators as, where the average dispersion energy of a pair in Bloch band n is given by n,avg = 2 E n (q)dq. From the Hamiltonian (29), by taking a mean-field type approximation for the composite operator will again result is the effective shift in the dispersion.
Appendix B: Testing the approximations

Small system analysis of a single pair
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 Eq. (23) and Eq. (25). 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 Fig.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 singleparticle states, the decay rates of various channel remains unchanged with respect to the case discussed in the paper.
Effect of Van der Waals (VdW) potential due to rotational level mixing 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 [50]. 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 E along the z direction. In the limit of (µE/ B e ) 1, one can write the rotational Hamiltonian in the M = 0 sector (M is the projection of angular momentum along the molecular axis) as, 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 rot is denoted by the index m = 0, 1, 2, ..... with energy E rot,m = (2m + 1) B e /d 2 θ and wavefunction is the Hermite polynomial of order m, N m is the normalization constant and the width d θ = [2 B e /µE] 1/4 . The rotational state of the polar molecule is denoted by the lowest energy rotational wavefuntion Φ 0 (θ) which induced a dipole moment of µ ind µ cos θΦ 2 0 (θ)dθ. 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, where the distance are in the units of a/π and the effective dimensionless VdW length VdW is given by, = where the maximum dipolar strength is given by D max = m b µ 2 /2 0 2 a. For RbCs molecule, the rotational constant is given by B e = 0.014cm −1 and the permanent dipole moment is given by µ = 1.27 Debye. Then for an electric field strength of E = 10kV/cm, the angular width reads d θ = 0.7. Correspondingly, the VdW length is given by VdW ≈ 0.17 when the lattice constant is a = 500nm. From this we can also define a short distance cutoff scale sr where rotational mixing effect of the dipoles becomes similar magnitude to the rotational splitting [50]. In our units, this cut off is given by sr ≈ .03 for dipolar strength D = 20. For length scales r > sr , the perturbative form of the VdW interaction in Eq.(30) remain valid and for r < 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 VdW,int ), continuous (H VdW,Bloch ) and discretecontinuous (H VdW,WB ) sector. The interaction in discrete sector is weak compare to the dipolar interaction. This can be easily seen by Fourier transforming Eq. (30), V VdW (k) ≈ 4 VdW k 3 F(k VdW ), with the function F ∼ 1. The widths of the Wannier functions in the momentum space are of the order of k ∼ 1. Then as 4 VdW ∼ 10 −3 1, we can neglect the VdW interaction in the discrete states compare to the dipolar strength in Eq. (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(−1/ VdW )ρ 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 Ref. [37], for RbCs molecue in the rotational ground state, the density of states is large, ρ ∼ 40. Accordingly, the decay rate will be proportional to the overlap which is of the order of 0.1E R which gives a timescale of ∼ 1.0ms. 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 ∼ √ E n and the corresponding strength of the VdW interaction in the continuum band n is in the order of − 4 VdW E 3/2 n . Whereas from Eq. (28), we find that the strength of the dipolar attraction in this band is of the order of ∼ D. For dipolar strength of D ∼ 20, the VdW interaction gets prominent only for very high Bloch bands with n cutoff 1/ sr ≈ 30. As such bands probes distance shorter than 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 Γ VdW The situation remains similar also for the discrete to continuous transitions. There the transitions happens between states with momentum k ∼ √ E n continuous states and discrete states with momentum k ∼ 1. As for continuous states, E n 1, the corresponding VdW discrete-continuous transition strength is of the order of − 4 VdW E 3/2 n . Subsequently, for intermediate momentum, the discrete-continuum transition is again dominated by the dipolar terms in Hamiltonian (5)and molecular sticking due to VdW interaction involves very high energy Bloch bands (or shorter distance) with band index n 30.
Accordingly, while integrating out the high Bloch bands, one get additional terms (equivalent to the term in Eq.(31)), n>ncut P n ij,q 1 q 2 as n 2 cut D, and the decay rate Γ VdW ∼ D for a lattice constant of 500nm and dipolar strength D ∼ 20. We have calculated Γ VdW from Ref. [37] but assuming temperature in the nano-Kelvin regime which suppresses the d-wave resonances.

Two-component Bose gas limit of Eq.(18)
We rewrite our Hamiltonian (18) in the conventional J 1 − J 2 form as, 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 [40]. Here we extend this treatment to include long-range dipolar interaction. To do that we transform the Hamiltonian to the momentum space, where the dispersion relation is given by q = 2J 1 cos qa + 2J 2 cos 2qa and the interaction energy in momentum space is given by, V (q) = U +2V ∞ n=1 cos nqa/n 3 , where the hardcore constraint is given by U → ∞. We only consider the dilute limit, µ → 0. When J 2 > J 1 /4, the dispersion relation has two minima at wavevectors, Qa = cos −1 [−J 1 /4J 2 ]. Around these minima, we can write the dispersion relation as, Q+k = Q + 2 k 2 /2m * , where m * is the effective mass. Then we expand the boson operator near the two minima, b k = φ 1,Q+k + φ 2,−Q+k + φ k , where φ 1 and φ 2 are the two-component Bose gas centered around momentum ±Q respectively, while φ k denotes the high momentum contribution, which is integrated out. Then one can re-express the Hamiltonian (33) in terms of the φ 1,2 which in position space reads, where T 1 and T 12 are renormalized interactions. To find these renormalized interactions, we first write down the full Bethe-Salpeter equation, In the dilute limit we can substitute Ω = 2µ. Then the respective renormalized interaction is given by, T 1 = T (Q, Q, 0) and T 12 = T (Q, −Q; 0) + T (Q, −Q; 2Q). Imposing the hard-core constraint with U → ∞, we get an additional equation, T (k, k ; p) k+p + k −p + Ω dp 2π = 1.
Due to the form of the interaction V (q), we expand the full renormalized interaction as, T (k, k ; q) = A 0 + n A n cos nqa, where the coefficients A 0 , A n depends on k, k . Putting this ansatz in Eq.(35), we get a set of coupled equations for m > 0, A m cos mpa cos m pa 2 p 2 /m * + Ω dp 2π .
and from the constraint condition, A m cos mpa 2 p 2 /m * + Ω dp 2π .
We found that in the limit of Ω → 0, the magnitude of the integral like cos mpa 2 p 2 /m * +Ω dp 2π falls off when m > 1. Then we get the following relation, A m cos(m − m )pa 2 p 2 /m * + Ω dp 2π .
We numerically find convergent solution for the A m by taking m max = 100.