{'fraction_non_alphanumeric': 0.048672340859283424, 'fraction_numerical': 0.01641098822689975, 'mean_word_length': 4.1393399685699315, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Particle-hole bound states of dipolar molecules in optical lattice

We investigate the particle-hole pair excitations of dipolar molecules in optical lattice, which can be described with an extended Bose-Hubbard model. For strong enough dipole-dipole interaction, the particle-hole pair excitations can form bound states in one and two dimensions. With decreasing dipole-dipole interaction, the energies of the bound states increase and merge into the particle-hole continuous spectrum gradually. The existence regions, the energy spectra and the wave functions of the bound states are carefully studied and the symmetries of the bound states are analyzed with group theory. For a given dipole-dipole interaction, the number of bound states varies in momentum space and a number distribution of the bound states is illustrated. We also discuss how to observe these bound states in future experiments.


I. INTRODUCTION
Cold and ultracold molecules are attracting more and more attention due to their broad applications in the fields of high-precision measurement, quantum chemistry, quantum information and many-body physics [1]. Recent years the experimental techniques for cold and ultracold molecules have been developed greatly. The homonuclear molecular Bose-Einstein condensation (BEC) was realized experimentally [2,3] and the quantum state with exactly one molecule at each site of an optical lattice was created with a Feshbach resonance and STIRAP (stimulated Raman adiabatic passage) techniques [4,5]. Meanwhile the ultracold heteronuclear molecules were also produced [6][7][8][9][10].
The heteronuclear molecules, such as SrO, RbCs, or NaCs, prepared in their electronic and vibrational ground states, have considerable permanent electric dipole moment. The dipole-dipole interaction between molecules can be generated by the external applied electric field [1]. Moreover, the magnitude of the dipole-dipole interaction can be controlled by the strength of the electric field [11,12]. The tunable long range dipole-dipole interaction may significantly modify the ground state and collective excitations of trapped condensates. For example, in trapped dipolar gases, it was shown theoretically [13] that the mean-field inter-particle interaction and, hence, the stability diagram are governed by the trapping geometry. With increasing dipolar interaction, the ground state of rotating atomic Bose gases undergoes a series of transitions between vortex lattices of different symmetries: triangular, square, "stripe", and "bubble" phases [14]. In rapidly rotating Fermion gas, the dipole-dipole interaction may even result in the fractional quantum Hall-like states [15].
The particles with long-range interactions in optical lattice can be described with extended Hubbard model [16][17][18]. Compared with the regular Bose-Hubbard model with on-site interaction, the extended Bose-Hubbard model has richer ground state phases, such as Mott insulator, particle density wave, superfluidity or supersolid phase [19][20][21][22][23][24][25][26]. On the other hand, the gapful particle and hole excitations in the insulating phase [27][28][29] may bind together and form bound states due to the long-range interactions. Although the excitons (holon-doublon pairs) in one-dimensional fermionic Hubbard model have been extensively studied [30], the excitons in higher dimensions, and especially, in bosonic Hubbard models are less studied.
In this paper, we study the particle-hole pair excitations of dipolar bosonic molecules in optical lattice, especially, the possible bound states due to the dipole-dipole interaction. The paper is organized as follows. In Sec. II we introduce the extended Bose-Hubbard model to describe the polarized bosonic molecules in optical lattice, and derive the eigen equations to describe the single particle-hole pair excitation. In Sec. III, the existence regions, the energies and the wave functions of the particle-hole bound states are calculated in one and two dimensions. The symmetries of the bound states are analyzed and possible experimental observation of these bound states is also discussed. A summary is presented in Sec. IV.

II. THE PARTICLE-HOLE BOUND STATES IN THE EXTENDED BOSE-HUBBARD MODEL
Considering only the nearest-neighbor interaction V to simulate the effect of dipole-dipole interaction, we write the extended Bose-Hubbard model as, FIG. 1. Sketches of optical lattice and interactions of ultracold heteronuclear molecules. (a) The landscape of optical lattice potential in the x-y plane. The square optical lattice potential Vo(x, y) = Vopt[sin 2 (kxx) + sin 2 (kyy)] can be created by counter propagating far detuned laser beams, where Vopt is the depth of lattice and k x(y) is the laser wave vector along x(y) direction. (b) The interactions between dipolar molecules in the optical lattice (see equation (1)). The arrows on the molecules indicate the polarization of the electric dipole moments.
where n r = b † r b r is the number operator at site r with b r (b † r ) the annihilation (creation) operator of particle. σ denotes the nearest neighbor vectors of site r. The onsite interaction U , the nearest neighbor interaction V and the hopping t are expressed as where W ( r) is Wannier function corresponding to the lowest energy band and V dd ( is the dipole-dipole interaction of two particles separated with a distance | x|. For the electric dipole-dipole interaction, C dd = d 2 ǫ0 with d the electric dipole moment, ǫ 0 is vacuum permittivity. We assume that the polarization of dipole moment is along the z axial direction. In the following studies, we consider particles confined in one-dimensional chain or two-dimensional square lattice and treat hopping term H 1 as a perturbation. In the atomic limit (t = 0), n r is a good quantum number and the eigenstates of H 0 can be written as direct product of number states. When the filling factor is one, V ≥ 0 and zV < U , the Mott insulating ground state can be written as |11111111... with the ground state energy E 0 = 1 2 N zV , where z is the coordination number and N the number of lattice sites. The excited states contain one or more particle-hole pairs. A single particle-hole pair state with a hole at r and a particle at r ′ can be represented as |r, r ′ = |110 r 111112 r ′ 11... . Taking into account the hopping term t, the particle and hole will move in the lattice. Similar to the two-magnon states in ferromagnetic system [31], the single particle-hole pair state can be written as a linear combination of |r, r ′ : where φ r, r ′ will be determined by solving the approximate eigen equation in single particle-hole pair subspace Calculating H| r, r ′ and considering the boundary condition φ r, r = 0 (a particle and a hole do not share the same lattice site), we get where δ r ′ − r,σ is the Kronecker delta function and ε = E − E 0 − U . The particle-hole excitation energy is ω = Owing to the translational invariance, it is convenient to apply a transformation φ r, and ρ = r ′ − r denote the center-of-mass and the relative coordinates respectively and K the total momentum of the particle and hole. At fixed K, the eigen equation (5) becomes Considering the finite range character of the nearest neighbor interaction V , we utilize the Green's function Energy spectra of the bound states of particlehole pair in one dimension for different nearest-neighbor interactions V = 2t (black circles), V = 3t (red triangles), V = 5t (black diamonds) and V = 8t (red squares) respectively (K = Kx). As comparison, the continuous spectrum of a particle-hole pair excitation is also shown (green shaded region).
approach to solve the eigen equation. To this end, we introduce an effective Hamiltonian where to describe the motion of a particle around a hole, where | ρ is basis set in relative coordinate spaces. The Kdependent H ph,0 describes the kinetic energy of the particle and V ph denotes the interaction between particle and hole. With |φ = ρ φ( ρ)| ρ , it is easy to reproduce Eq. (6) from the eigen equation H ef f |φ = ε|φ . After a Fourier transformation, the kinetic energy is obtained as where q is the relative momentum of particle and hole. When a particle is adjacent to a hole, the energy offset −V in V ph may result in the formation of particle-hole bound states.

Introducing retarded Green functions
In the real space, the Green's function G has a general form where the contributions of the continuous spectra are neglected and E m 's are the eigenvalues of the bound states with φ m 's the corresponding eigenfunctions [32]. We can determine the eigenvalues and eigenfunctions by analyzing the poles and residues of the obtained Green's function. From Eq. (7), we have For specific ρ′ = 0 or σ, Eq. (9) can be reduced to a set of simultaneous linear equations. Green's functions ρ ′ |G| ρ with ρ′ = 0 or σ can thus be obtained exactly. The residues of 0|G| ρ is always vanishing for non-vanishing bound states energies ε in our calculations, which is consistent with the boundary condition φ(0) = 0. We present our results in the next section.

A. Bound states in one dimension
In the Mott insulating phase, the excitation energy of a single particle/hole was calculated as with the dynamical Gutzwiller approach [29]. Up to the first order of t U , the single particle-hole pair excitation energy is with the total momentum K = k p + k h and relative momentum q = ( k h − k p )/2. For the extended Bose-Hubbard model, however, a particle-hole pair may form a bound state due to the interaction (See Eq. (7)). In the following calculations, we take the lattice constant as 1 and use t > 0 as the unit of energy. Conditions of t/U ≪ 1 and zV < U are assumed to make sure that the ground state of the system is the deep Mott insulating phase. In one dimensional case, ε 0 ( K, q) = −2tx cos(q + θ x ) with tx = t 9cos 2 ( Kx 2 ) + sin 2 ( Kx 2 ) and θx = (see Fig.2). Compared with particle-hole pairs in fermion Hubbard model with half filling [34], the minimum of band width of the continuous spectrum is non-zero at K = ±π. This is because that there are no particle-hole symmetries in Bose-Hubbard model.
For a specific nearest-neighbor interaction V and total momentum K, we search for bound state solutions outside the continuum. In the case of ε < −2tx, the free Green function is calculated as [33] , γ = e iθx and β = e −κ with κ = arccosh( |ε| 2tx ). With ρ ′ = 0, ±1, Eq. (9) is reduced to 3 linear equations. The Green functions of −1|G|ρ , 0|G|ρ and 1|G|ρ can be exactly obtained as where θ(ρ) is Heaviside step function. From these Green functions, we obtain two degenerate bound states corresponding to the particle on the left (right) of the hole respectively. The bound state energy is Accordingly, β is simplified as β = tx V < 1 and the condition for the existence of the bound state is obtained as When 0 ≤ V ≤ t, the interaction between particle and hole is too weak to bind them together and no bound state is found. When t < V < 3t, two degenerate bound states exist in the regions of 2arccos The region get larger with the increase of V . When V ≥ 3t, bound states may be found in the whole first Brillouin zone. In Fig.2, we show the spectra of the bound states for the interaction V = 2t, 3t, 5t and 8t, respectively. The energies of the bound states decrease with increasing V . At V = 2t, the calculated existence intervals are 0.58π < |K x | ≤ π.
From the residues of the Green function, we can also extract the bound state wave functions. For example, the bound state with a particle on the right of a hole is written as where C = 1 − t 2 x V 2 is the normalized constant. In Fig.3, we show the probability distribution of the wave function for K = 0 and V = 8t.
The mean size of the bound state is calculated as x . In Fig.4, we show the mean size of the bound state as a function of the nearestneighbor interaction V at K x = 0. The stronger is the interaction, the smaller is the size and the closer do the particle and hole bind.

B. Bound states in two dimension
Before presenting the numerical results, we discuss the symmetries of the effective Hamiltonian H ef f in Eq. (7) in details. After a gauge transformation ) and θŷ =arctan( 1 3 tan( Ky 2 )), we could remove the phase factors in the Hamiltonian and get where σ x and σ y denote the nearest neighbors along the x and y directions, respectively, and tx = t 9cos 2 ( Kx 2 ) + sin 2 ( Kx 2 ) and tŷ = t 9cos 2 ( Ky 2 ) + sin 2 ( Ky 2 ) are effective hoppings along the x and y directions.
When K x = ±K y , tx = tŷ, H ′ ef f has symmetries of D 4 group. According to the irreducible representations of this group, the bound states can be classified and labeled with s, p x(y) and d x 2 −y 2 wave, respectively. Among them, the s wave belongs to an identical representation A 1 of D 4 group, the degenerate p x and p y waves form a two dimensional irreducible representation E of D 4 group, and the d x 2 −y 2 wave belongs to an irreducible representation B 1 of D 4 group.
When K x = ±K y , tx = tŷ, the symmetry reduces to D 2 , a subgroup of D 4 . For simplicity, we still label the four bound states with s, p x(y) and d x 2 −y 2 . Differently, here all the s and d wave belong to identical representations A 1 of D 2 group, while p x (p y ) wave belongs to an irreducible representation B 3 (B 2 ) of D 2 group. The degeneracy of p x and p y waves is lifted.
With ρ ′ = (0, 0), (±1, ±1), Eq. (9) reduces to five linear equations in two dimensional case. The free Green function G 0 ( ρ, ρ ′ ) can be expressed with elliptic integrals [32,35]. Although no brief solutions of < ρ|G|ρ ′ > could be found, we may factorize the particle-hole bound state equations and analyze the existence conditions for the bound states along the symmetric lines of K x = ±K y = K [35]. Considering the asymptotic behaviours of elliptic integrals, we get the thresholds of the interaction V as follows: V cr,s (K) = tx for s wave, V cr,p (K) = 2π 2π−4 tx for p x(y) wave and V cr,d (K) = 2π 8−2π tx for d x 2 −y 2 wave. From these thresholds, the existence region of every bound state is obtained respectively. For s wave, we have When V > 3t, the s wave bound state may be found in the whole Brillouin zone. The existence region of p wave bound states is When V > 3π π−2 t, the existence region extends to the whole Brillouin zone. For d wave, existence region is expressed as When V > 3π 4−π t, the d wave existence interval is the whole Brillouin zone. Now we present the numerical results in two dimensions. At K = (0, 0), four bound states could be found when V ≥ V cr,d (0) = 10.97t. The wave functions with s−, p x − and d x 2 −y 2 −wave symmetry for V = 12t are shown in Fig. 5. With the decrease of V , the highest d x 2 −y2 -wave, the degenerate p x(y) -waves and the lowest s-wave merge into the continuous spectrum one by one. Finally, all the bound states disappear when V ≤ V cr,s (0) = 3t. In Fig. 6, we illustrate the variations of the bound state energies with the changes of V .
We then search for bound state solutions along the line of K x = K y = K. As shown in Fig. 7, there are four bound states for all K ∈ [−π, π] at V = 12t. At V = 6t, s-wave exists in the whole region, while p x(y)and d x 2 −y 2 -wave states appear at 0.52π < |K| ≤ π and 0.69π < |K| ≤ π respectively. Decreasing V further, we find that the regions of the bound states shrink, and all the bound states disappear when V < V cr,s (±π) = 1t. Similar results are obtained along the line of K x = −K y = K.
Away from the lines of K x = ±K y , the degeneracy of p x and p y bound states is lifted, as mentioned before. At a given V , different bound states have different existence regions in K space. Consequently, the number of bound states vary in the (K x , K y ) space. As an example, we show the number distribution of the bound states for V = 6t in Fig. 8.

C. Discussion of observations of the bound states
Inelastic light scattering directly measures the dynamical structure factor S( q, ω), the Fourier transformation of density correlations. Bragg spectroscopy has been proposed to detect quantum phases in optical lattice and successfully applied to measure the excitation spectra (phonons of BEC), the composition of the excitations and the Higgs-type amplitude mode in the superfluid condensate, as well as the particle-hole excitation energies in the Mott-insulator state [36][37][38][39][40][41][42]. When this technique is utilized in the study of ultra-cold polar molecules in the optical lattice, we may expect extra resonance peaks corresponding to the particle-hole bound states lying outside the particle-hole continuum.
Compared with the traditional solid state counterparts, the deep Mott insulating state with t/U ≪ 1 and zV < U can be realized by tuning the interaction parameters t, U and V in ultra-cold dipolar molecules in the optical lattices. In deep optical lattice the hopping is approximately evaluated as t = (4/ √ π)E r (V opt /E r ) 3/4 exp[−2(V opt /E r ) 1/2 ], with V opt = 18 ∼ 20E r the optical lattice depth and E r the lattice recoil energy [43]. Under the condition of as a ⊥ ≪ 1, where a s is the s-wave scattering length and a ⊥ = mω ⊥ with ω ⊥ transverse trapping frequency in one dimension or trapping frequency of the z direction in two dimensions, the on-site interaction is estimated as U 1D = 2 π ω ⊥ a s /l 0 and U 2D = √ 8π 2 as maz(2π)l 2 0 with l 0 = ( Er Vopt ) 1 4 a π and a the lattice constant [43,44].
Taking Bose molecule 23 Na 7 Li for example, we may estimate the typical values of t, U and V for observing the particle-hole bound state. The molecule prepared in the ground state has a permanent electric moment d = 0.58D [45]. With a lattice constant a ∼ 0.5µm, a s /a ∼ 0.01 and the transversal tapping frequency ω ⊥ ∼ 2π × 10 4 Hz, we have t ∼ 1nK, U 1D ∼ 25nK and U 2D ∼ 50nK. The nearest-neighbor interaction V can be tuned as V ∼ 10nK with an applied electric field.

IV. SUMMARY
In summary, we have investigated bound states of particle-hole pair resulting from the dipole-dipole interaction between polar molecules in the optical lattice. For a large enough dipole-dipole interaction, two degenerate bound states, which correspond to a particle on the left and the right of a hole, are shown to exist in one dimension. While in two-dimensional case, four bound states, with s-, p x(y) -and d x 2 −y 2 -symmetry respectively, are found along the lines of K x = ±K y . Away from the lines of K x = ±K y , the degeneracy between p x and p y waves is lifted. With decreasing the nearest-neighbor interaction V , the energies of the bound states increase and merge into the particle-hole continuum gradually. The wave functions, the dispersion relations and the existence regions of the bound states are studied in details. For a given nearest-neighbor V , the number of bound states is different in different regions of K space and a number distribution of bound states is given for the nearestneighbor V = 6t in two dimensions. The possible experimental observation of the particle-hole bound states is also discussed.
Electron-hole bound state excitations (excitons) have been extensively studied for many years and very recently, the excitonic Bose-Einstein condensates have been realized experimentally in semiconductors [46,47]. We hope our study on the particle-hole bound states in bosonic systems would enrich our understanding of elementary excitations in quantum many-body systems and stimulate more efforts on the bound state phenomena in various systems such as magnets, superconductors and atomic systems.