Multi-mode Bose-Hubbard model for quantum dipolar gases in confined geometries

TL;DR

This paper develops a multi-mode Bose-Hubbard model to study quantum phases of ultracold polar molecules in confined geometries, capturing quantum fluctuations and phase transitions at the linear-zigzag instability with limited computational resources.

Contribution

The authors derive a generalized multi-mode Bose-Hubbard model for dipolar gases in confined geometries, incorporating quantum fluctuations and analyzing phase diagrams at the zigzag instability.

Findings

Quantum fluctuations significantly influence ground state properties.

Aspect ratio of transverse trap frequencies controls quantum phases.

Model accurately captures key features with limited computational effort.

Abstract

We theoretically consider ultracold polar molecules in a wave guide. The particles are bosons, they experience a periodic potential due to an optical lattice oriented along the wave guide and are polarised by an electric field orthogonal to the guide axis. The array is mechanically unstable by opening the transverse confinement in the direction orthogonal to the polarizing electric field and can undergo a transition to a double-chain (zigzag) structure. For this geometry we derive a multi-mode generalized Bose-Hubbard model for determining the quantum phases of the gas at the mechanical instability taking into account the quantum fluctuations in all directions of space. Our model limits the dimension of the numerically relevant Hilbert subspace by means of an appropriate decomposition of the field operator, which is obtained from a field theoretical model of the linear-zigzag instability. We determine the phase diagrams of small systems using exact diagonalization and find that, even for tight transverse confinement, the aspect ratio between the two transverse trap frequencies controls not only the classical but also the quantum properties of the ground state in a non-trivial way. Convergence tests at the linear-zigzag instability demonstrate that our multi-mode generalized Bose-Hubbard model can catch the essential features of the quantum phases of dipolar gases in confined geometries with a limited computational effort.

Figures

Equations

$$U_d(\mathbf{r})=\frac{p^2}{r^3}-\frac{3(\mathbf{p}\cdot \mathbf{r}{)}^2}{r^5}.$$

$$V_{\mathrm{trap}}(y,z)=\frac{1}{2}M({\omega }_y^2{y}^2+{\omega }_z^2{z}^2),$$

$$V_{\mathrm{opt}}(x)=V_L{\cos }^2(\pi x/a),$$

$$\mathcal{H}=$$

$$+\frac{1}{2}\int \mathrm{d}{}^{3}r\int \mathrm{d}{}^3{r}^{\prime }{\Psi }^{†}(\mathbf{r}){\Psi }^{†}({\mathbf{r}}^{\prime })U(\mathbf{r}-{\mathbf{r}}^{\prime })\Psi ({\mathbf{r}}^{\prime })\Psi (\mathbf{r}),$$

$$U(\mathbf{r})=U_d(\mathbf{r})+U_g(\mathbf{r}),$$

$${\theta }_0(z)=\frac{1}{\sqrt{\sqrt{\pi }{\sigma }_z}}\exp \left(-\frac{z^2}{2{\sigma }_z^2}\right),$$

$$H_{\text{loc}}^{(j)}=-\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial y_j^2}+\frac{1}{2}M{\omega }_y^2{y}_j^2+U_{\mathrm{pin}}^{(j)},$$

$$U_{\mathrm{pin}}=$$

$$\frac{1}{2}M{\omega }_y^2{y}_j^2+U_{\mathrm{pin}}^{(j)}\simeq \frac{1}{2}\left(M{\omega }_y^2-\frac{p^2}{a^5}{\mathcal{M}}_1\right)y_j^2+\frac{p^2}{2a^7}{\mathcal{M}}_2{y}_j^4,$$

$${\mathcal{M}}_{q=1,2}=\frac{\left(2^{3+2q}-1\right)\Gamma (q+\frac{3}{2})}{q!4\Gamma (\frac{3}{2})}\zeta (3+2q),$$

$${\omega }_y^{(c)}=\sqrt{{\mathcal{M}}_1{p}^2/(M{a}^5)}.$$

$$\Psi (\mathbf{r})=\sum _{j,m,n}{w}_j(x){\varphi }_m(y){\theta }_0(z)a_{j,m},$$

$$H_{BH}=H^x+H^y+H^{xy},$$

$$H_m^x$$

$${ϵ}^x=\int dx{w}_j(x)\left[-\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial x^2}+V_L{\cos }^2\left(\frac{\pi x}{a}\right)\right]w_j(x),$$

$$J^x=\int dx{w}_j(x)\left[\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial x^2}-V_L{\cos }^2\left(\frac{\pi x}{a}\right)\right]w_{j+1}(x),$$

$$U_m^x=\int d^3{r}_1\int d^3{r}_{2}U({\mathbf{r}}_1-{\mathbf{r}}_2)w_j^2(x_1)w_j^2(x_2){\varphi }_m^{\ast 2}(y_1){\varphi }_m^2(y_2){\theta }_0^2(z_1){\theta }_0^2(z_2).$$

$$V_m^x=\int d^3{r}_1\int d^3{r}_2{U}_d({\mathbf{r}}_1-{\mathbf{r}}_2)w_j^2(x_1)w_{j+1}^2(x_2)|{\varphi }_m(y_1){\varphi }_m(y_2){|}^2{\theta }_0^2(z_1){\theta }_0^2(z_2),$$

$$P_m^x=\int d^3{r}_1\int d^3{r}_2{U}_d({\mathbf{r}}_1-{\mathbf{r}}_2)w_j(x_1)w_{j+1}(x_1)w_j(x_2)w_{j+1}(x_2)|{\varphi }_m(y_1){\varphi }_m(y_2){|}^2{\theta }_0^2(z_1){\theta }_0^2(z_2),$$

$$H_j^y$$

$${ϵ}_m^y=\int dy{\varphi }_m(y{)}^{\ast }\left[-\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial y^2}+V_{\mathrm{trap}}(y)\right]{\varphi }_m(y),$$

$$J_{mn}^y=\int dy{\varphi }_m(y{)}^{\ast }\left[\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial y^2}-V_{\mathrm{trap}}(y)\right]{\varphi }_n(y),$$

$$U_{l,m,p,q}^y=\int d^3{r}_1\int d^3{r}_{2}U({\mathbf{r}}_1-{\mathbf{r}}_2)w_j^2(x_1){\theta }_0^2(z_1)w_j^2(x_2){\theta }_0^2(z_2){\varphi }_l^{\ast }(y_1){\varphi }_m(y_2){\varphi }_n^{\ast }(y_1){\varphi }_q(y_2).$$

$$H^{xy}=\sum _{j,m}{H}_{j,m}^{xy}$$

$$H_{j,m}^{xy}=\frac{1}{2}{\sum _{{\ell }_1,{\ell }_2,{\ell }_3}}^{\prime }{\sum _{n_1,n_2,n_3}}^{\prime }{V}_{m,n_1,n_2,n_3}^{j,{\ell }_1,{\ell }_2,{\ell }_3}{a}_{j,m}^{†}{a}_{{\ell }_1,n_1}^{†}{a}_{{\ell }_2,n_2}{a}_{{\ell }_3,n_3},$$

$$V_{m,n,q,r}^{j,{\ell }_1,{\ell }_2,{\ell }_3}=\int d^3{r}_1\int d^3{r}_2{U}_d({\mathbf{r}}_1-{\mathbf{r}}_2)w_j(x_1)w_{{\ell }_1}(x_1)w_{{\ell }_2}(x_2)w_{{\ell }_3}(x_2){\varphi }_m^{\ast }(y_1){\varphi }_{n_1}^{\ast }(y_2){\varphi }_{n_2}(y_2){\varphi }_{n_3}(y_1){\theta }_0^2(z_1){\theta }_0^2(z_2).$$

$$\Delta n_j=⟨\delta n_j⟩,$$

$$\delta n_j^M=\sum _m(n_{j,m}-⟨n_{j,m}⟩{)}^2.$$

$$\varphi =\sum _{j,m}⟨a_{j,m}^{†}(a_{j+1,m}+a_{j-1,m})⟩.$$

$$S_x(q_x)=\frac{1}{N^2}\sum _{j,l=1}^N{\mathrm{e}}^{ı(j-l)q_{x}a}\sum _m\left(⟨n_{j,m}{n}_{l,m}⟩-⟨n_{j,m}⟩⟨n_{l,m}⟩\right).$$

$$\Phi =\sum _{j,m}⟨a_{j,m}^{†}{a}_{j,m}^{†}{a}_{j+1,m}{a}_{j+1,m}+H.c.⟩.$$

$$S_y\left(\frac{\pi }{a}\right)=\frac{1}{N^2}\sum _{j\ne l}(-1{)}^{j-l}⟨y_j{y}_l⟩$$

Full text

I Introduction

Dipolar bosonic gases offer a laboratory for studying the interplay of finite-range interactions and quantum fluctuations Lahaye ; Jin . The study of their dynamics in optical lattices, moreover, allows one to realize and characterize strongly-correlated states of ultracold matter Lahaye ; Jin ; BlochRMP ; Ferlaino2016 . The essential features of the quantum phases of ultracold dipoles in optical lattices are believed to be captured by the so-called extended Bose Hubbard Model EBH . This model reduces to the single-band Bose-Hubbard model for vanishing dipolar coupling, which at commensurate densities exhibits the Mott-Insulator to Superfluid quantum phase transition Bose-Hubbard ; BlochRMP ; FisherAndFisher . For finite strengths of the dipolar interactions, in addition, it includes a finite-range interaction term that favours the appearance of diagonal long-range order Goral ; Menotti ; Batrouni2013 ; Deng2013a ; Deng2013b ; Batrouni2014 .

In three dimensions the anisotropic nature of the dipolar interaction is reflected in the properties of the Bose-Hubbard coefficients and can be analysed by orienting the dipolar structure by means of an external field Ferlaino2016 . When the motion is confined on a plane, instead, the mutual dipolar interaction can be made effectively isotropic and repulsive by orienting the dipoles perpendicularly to the plane itself. In this regime crystalline structures can emerge from the competition between the external confinement and the particles repulsion Astrakharchik:2007 ; Buechler:2007 ; Jin .

One exemplary situation is the linear-zigzag instability. This instability can be observed by tuning the frequency of the transverse trap, confining the dipoles along the array, and is illustrated in Fig. 1 for a chain of dipoles in an optical lattice. For an incompressible chain the transition is continuous and the classical order parameter is the transverse displacement Fishman ; Astrakharchik ; Altman , while the quantum linear-zigzag transition is of the same universality class as the Ising model in transverse field Shimshoni . When the chain is compressible, instead, the classical transition becomes of weak first order Cartarius , while the corresponding quantum behaviour is yet unexplored. In these respects the model we consider is peculiar, since the compressibility results from the interplay between interactions and quantum fluctuations and can be thus tuned by changing the lattice depth of the transverse confinement. Furthermore, previous literature pointed out that quantum fluctuations in the transverse directions can substantially modify the effective interaction the dipoles experience along the axis Sinha ; Deuretzbacher ; Recati ; Sowinski . The description of the structural instability, therefore, requires the development of a suitable model which describes spatial selforganization in the transverse direction, while the dipoles density is periodically modulated along $x$ and quantum fluctuations in all directions of space are appropriately taken into account.

In this work we systematically derive a multi-mode extended Bose-Hubbard (EBH) model which is particularly apt to describe the phase diagram deep in the linear chain as well as close to the linear-zigzag instability. Our model is derived by identifying a suitable basis for the transverse excitations, which is obtained using the field theoretical description of the linear-zigzag instability Shimshoni ; Silvi ; Podolsky . The dynamics takes into account the anisotropic nature of the dipolar interaction by calculating the integrals defining the EBH coefficients in three dimensions, thus including the fluctuations in the three-dimensional space.

This manuscript is organized as follows. In Section II we report the detailed derivation of the multi-mode extended Bose-Hubbard model. In Section III we determine the phase diagrams of small systems for different aspect ratios of the transverse confinement using exact diagonalization. Moreover, we test the convergence of our basis choice at the linear-zigzag instability. The conclusions are drawn in Sec. IV while the Appendix reports calculations complementing the material of Sec. II.

II Derivation of the multi-mode Bose-Hubbard model

We consider a gas of identical dipolar molecules with mass $M$ and dipolar moment ${\bm{p}}$, interacting via the dipolar potential $U_{d}(\bm{r})$ (with $\bm{r}=(x,y,z)$ the distance between the centers of mass of two molecules):

[TABLE]

An external electric field, moreover, aligns the dipoles along the $z$ direction. The molecules form an array along the $x$-axis due to the tight confinement of an external harmonic trap,

[TABLE]

where the trap frequencies $\omega_{y}$ and $\omega_{z}$ are chosen such that $\omega_{y}\ll\omega_{z}$. The frequency $\omega_{y}$ is assumed to take value close to the critical value $\omega_{y}^{(c)}$, at which the linear-zigzag instability occurs in the mean-field model Astrakharchik ; Shimshoni .

The dipoles are ultracold and obey the Bose-Einstein statistics. They also interact via $s$-wave van-der-Waals collisions and occupy the lowest bands of an optical lattice along the $x$ direction,

[TABLE]

where $V_{L}$ is the lattice depth and $a$ the lattice constant. Their state is described in second quantization by means of the bosonic field operators $\Psi(\bm{r})$, $\Psi(\bm{r}^{\prime})^{\dagger}$, with $[\Psi(\bm{r}),\Psi(\bm{r}^{\prime})^{\dagger}]=\delta^{(3)}(\bm{r}-\bm{r}^{\prime})$ and is governed by Hamiltonian $\mathcal{H}$, which reads

[TABLE]

where $V_{\rm tot}(\bm{r})=V_{\rm trap}(y,z)+V_{\rm opt}(x)$. The interaction potential is the sum of the dipolar and of the contact interaction:

[TABLE]

where $U_{g}(\bm{r})=g\delta^{(3)}(\bm{r})$ describes the $s$-wave scattering contribution, with $g=4\pi\hbar^{2}a_{S}/M$ and $a_{S}$ the $s$-wave scattering length.

II.1 Mode expansion of the bosonic field operator

In order to derive a convenient multi-mode EBH model we use a suitably-chosen mode expansion. We first assume that the molecules are tightly bound at the minima of the optical lattice and we perform the single-band approximation. We thus denote by $w_{j}(x)$ the real-valued Wannier function at site $j$ for the motion of a particle of mass $M$ moving along $x$ and experiencing the potential $V_{\rm opt}(x)$. The motion along the $z$ axis is assumed to be in the ground state of the harmonic oscillator at frequency $\omega_{z}$ with wave function $\theta_{0}(z)$:

[TABLE]

and $\sigma_{z}=\sqrt{\hbar/(M\omega_{z})}$. The motion along $y$ is instead decomposed into the basis $\{\phi_{m}(y)\}$ which diagonalizes an effective local Hamiltonian along $y$ according to a procedure first developed in Ref. Silvi . This effective Hamiltonian includes the harmonic oscillator in the $y$ direction as well as the effective potential along $y$ due to the dipolar interactions. At the linear-zigzag structural transition the Hamiltonian describes an effective $\varphi^{4}$ model on a lattice, where the transition point at fixed linear density is given by the transverse trap frequency $\omega_{y}^{(c)}$. In detail, at site $j$ the Hamiltonian reads

[TABLE]

where $U_{\rm pin}^{(j)}$ is the local component of the dipole-dipole interaction $U_{\rm pin}$ describing the motion along $y$ assuming the dipoles are pinned at the classical equilibrium positions in the $x-z$ plane. Specifically,

[TABLE]

and $U_{\rm pin}=\sum_{j}U_{\text{pin}}^{(j)}+\sum_{[j\neq l]}H_{\text{int}}^{(j,l)}$. The effective $\varphi^{4}$ model is found close to the structural instability, where $|y_{j}|\ll 1$, when discarding the coupling with the axial modes due to the term $U_{\rm comp}=U_{d}-U_{\rm pin}$. In this limit the potential of Eq. (6) can be cast in the form Shimshoni ; Silvi :

[TABLE]

while the relevant terms of the sum $\sum_{[j\neq l]}H_{\text{int}}^{(j,l)}$ are $H_{\text{int}}^{(j,j+1)}=p^{2}/(2a^{5})\mathcal{N}_{1}\left(y_{j}+y_{j+1}\right)^{2}$ and describe an effective nearest-neighbour interaction. For completeness, we report the explicit form of the dimensionless coefficients: $\mathcal{N}_{1}=(9/4)\zeta(3)$, and

[TABLE]

with $\zeta(\ell)$ Riemann’s zeta function and $\Gamma(z)$ the Gamma’s function Abramowitz . From potential (7) one directly determines the mean-field critical frequency $\omega_{y}^{(c)}$, at which the chain becomes mechanically unstable Astrakharchik ; Shimshoni ; Silvi :

[TABLE]

We note that the basis $\{\phi_{m}(y)\}$ is found by numerically diagonalizing Hamiltonian $H_{\text{loc}}^{(j)}$ at site $j$, without performing any Taylor truncation of the local potential $U_{\rm pin}^{(j)}$. For $\omega_{y}\gg\omega_{y}^{(c)}$ we checked that it is well approximated by the eigenbasis of the harmonic oscillator at frequency $\omega_{y}$. For $\omega_{y}\simeq\omega_{y}^{(c)}$, instead the eigenbasis significantly differs from the oscillator eigenstates Silvi .

Using these prescriptions we decompose the field operator as

[TABLE]

where $a_{j,m}$ is the bosonic operator which annihilates a particle at site $j$ and in the local quantum state $|m\rangle$, and $[a_{j,m},a^{\dagger}_{\ell,n}]=\delta_{j,\ell}\delta_{n,m}$.

II.2 Multi-mode Bose-Hubbard model

The multi-mode EBH model $H_{BH}$ for our study is obtained by substituting Eq. (9) in the field operators of Hamiltonian (4), by integrating out the position variables and by keeping only nearest-neighbor interactions. The resulting EBH model exhibits a number of terms of different origin, which are conveniently identified by writing $H_{BH}$ as

[TABLE]

where the three terms describe the axial and transverse motion, as well as their mutual interaction, respectively.

II.2.1 The axial motion

The axial motion can be cast into the sum over the transverse bands labeled by the quantum number $m$, $H^{x}=\sum_{m}H^{x}_{m}$, with

[TABLE]

where we used a notation which can be put in direct connection with the EBH model of Ref. Sowinski for the single band case ($m=0$). Here, $n_{j,m}=a_{j,m}^{\dagger}a_{j,m}$ denotes the particle number at site $j$ and with quantum number $m$. The first three terms on the right-hand side (RHS) are the onsite energy, with $\epsilon^{x}$ the single-particle energy in the lattice, the hopping along the axis scaled by the hopping coefficient $J_{x}$, and the onsite interaction including the contribution of the dipolar term. Their explicit form is

[TABLE]

All other terms are solely due to dipole-dipole interaction and are the dipole blockade, whose strength is scaled by the coefficient $V^{x}_{m}$, the pair-hopping term, scaling with $P^{x}_{m}$, and the density-dependent tunnelling, proportional to $T^{x}_{m}$. These latter coefficients depend on the transverse quantum state $m$ and read

[TABLE]

When the transverse motion is in the ground state, namely, for $m=0$, Hamiltonian $H^{x}_{m}$ reduces to the model studied in Ref. Sowinski . If in addition one discards the pair-hopping and the density-dependent tunneling terms, then $H^{x}_{m}$ corresponds to the so-called extended Bose-Hubbard model, whose phase diagram has been extensively analysed in Refs. Berg ; Batrouni2013 ; Batrouni2014 .

II.2.2 Transverse motion

The EBH term for the Hamiltonian governing solely the motion along $y$ takes the form $H^{y}=\sum_{j}H^{y}_{j}$ and is local in the site $j$. Each term of the sum reads

[TABLE]

where ${\sum_{l,m,n,q}}^{\prime}$ indicates that at least one of the indices $l,m,p,q$ is different from the others. The coefficients are independent of the lattice site $j$ since the Hamiltonian is invariant per discrete translation (for periodic boundary conditions). Here, the eigenenergy $\epsilon_{m}^{y}$ and the tunneling term $J^{y}_{m,n}$ read

[TABLE]

while the interaction term $U^{y}_{l,m,p,q}$ takes the form

[TABLE]

We remark that also the onsite term $U^{y}_{m,m,m,m}=U^{x}_{m}$ contributes in determining the transverse motion. We arbitrarily assigned this term to the axial EBH Hamiltonian $H^{x}$ and did not include it in Eq. (18) in order to avoid double-counting in the resulting EBH Hamiltonian $H_{BH}$, Eq. (10).

II.2.3 Coupling between axial and transverse degrees of freedom

Finally, $H^{xy}$ describes the interaction between excitations along the $x$ and the $y$ direction, it is solely due to the dipolar interaction and can be written as

[TABLE]

where

[TABLE]

and describes a four-vertex type of interaction. We note that ${\sum_{\ell_{1},\ell_{2},\ell_{3}}}^{\prime}$ (${\sum_{n_{1},n_{2},n_{3}}}^{\prime}$) means that at least one of the indices $\ell_{1},\ell_{2},\ell_{3}$ has to be different from $j$ (respectively, at least one of the indices $n_{1},n_{2},n_{3}$ has to be different from $m$). Due to the tight-binding assumption, $\ell_{1},\ell_{2},\ell_{3}=j,j+1$ or $\ell_{1},\ell_{2},\ell_{3}=j,j-1$. The coefficients are found by performing the integral:

[TABLE]

Term $H^{xy}$ contains two physically relevant contributions. One contribution leads to the interaction term $H_{\text{int}}^{(j,j+1)}$ of the $\varphi^{4}$ model. The other is a coupling between axial and transverse modes, which becomes relevant when the chain is compressible Cartarius . When the chain is incompressible, for hard-core bosons and unit filling the action of Hamiltonian $H^{y}+H^{xy}$ can be reduced to an effective $\varphi^{4}$ model and the linear-zigzag transition is of the same universality class of the Ising model in transverse field Shimshoni ; Silvi .

II.2.4 Determination of the Bose-Hubbard coefficients

The coefficients corresponding to the interaction terms in the EBH model (Eqs. (14)-(17), (21), (24)) explicitly depend on the confinement in the $z$ direction, which enters through the wave function $\theta_{0}(z)$ and specifically through the size of the wave function $\sigma_{z}$. We perform the integrals first analitically, by integrating out the $z$ variable in Fourier space, then numerically. The details are reported in the Appendix. All other coefficients, which involve the integrals over two variables, are evaluated numerically.

The dependence on the size of the trap along $z$, where the dipolar interaction is attractive, turns out to be relevant for certain parameter regimes, even if the motion is confined in the orthogonal plane Sowinski ; Recati . In Fig. 2(a) we can observe that increasing the size of the quantum fluctuations along $z$ can change the on-site interaction from being repulsive (positive coefficient) to become attractive (negative coefficient). Figures 2(b) and (c) show that varying $\sigma_{z}$ can substantially modify the strength of the density-assisted tunneling and of the pair tunneling terms, respectively. These results, moreover, highlight that there is an important interplay between the fluctuations along $y$ and $z$ which significantly affects the behaviour of the coefficients of $H^{x}$, and thus could change the phase of a quasi-one dimensional system of dipolar bosons.

III Quantum phases of small systems

We now test the predictions of the multi-mode EBH model we derived by determining the quantum ground state as a function of the various parameters, as specified below. For this purpose we use exact diagonalization and assume periodic boundary conditions along $x$. This procedure limits us to small system sizes, yet it allows us to gain some insight into the possible phases one can observe. Moreover, it allows us to verify that our model reproduces correctly limiting cases analysed in the literature. This also provides us a point of comparison for future more elaborated numerical analysis based on Density Matrix Renormalization Group Silvi . In this work we are specifically interested in determining the phase diagram as a function of (i) the depth of the optical lattice $V_{L}$, (ii) the $s$-wave scattering length, (iii) the transverse frequency $\omega_{y}$, (iv) the strength of the dipole-dipole interactions, (v) the size of the fluctuations along $z$. In this section we discuss the observables, which permit us to identify the quantum phases, and determine the phase diagrams in several limiting cases.

III.1 Observables

For a system of few sites we identify whether a phase is compressible by means of the local compressibility $\Delta n_{j}$, which is the expectation value over the ground state of the observable $\delta n_{j}$ and reads

[TABLE]

For a single mode EBH model, $\delta n_{j}=\delta n_{j}^{0}$ with $\delta n_{j}^{0}=n_{j}-\langle n_{j}\rangle$, and a phase is classified as incompressible when $\Delta n_{j}$ vanishes at all sites $j$. For our multi-mode EBH model we use $\delta n_{j}=\delta n_{j}^{M}$, where

[TABLE]

According to this criterion a phase is incompressible when $\Delta n_{j}=0$ at all sites $j$, like in the single-band case.

Off-diagonal order is revealed by the non-vanishing value of the off-diagonal correlations (one-particle correlation function) $\phi$, which we define for the multi-mode EBH model as:

[TABLE]

The dipole blockade, scaling with coefficient $V_{m}^{x}$, favours the formation of a density modulation along $x$, which is signaled by a non vanishing value of the static structure form factor $S_{x}(q_{x})$ at wave number $q_{x}=\pi/a$. For the multimode EBH model we consider the structure form factor

[TABLE]

The value $S_{x}(\pi/a)\neq 0$ signals the formation of a structural order. We denote the phase by super-solid (SS) when this occurs in a compressible phase with non-vanishing off-diagonal correlations. The phase is instead charge-density wave (CDW) when incompressible Goral ; Sowinski .

Additionally, pair tunnelling terms are expected to favour the onset of what has been denoted by pair superfluidity Sowinski , and which shall be signaled by a non-vanishing expectation value of the pair-correlation function, defined as:

[TABLE]

These quantities have been used in the literature to characterize the phases of the one-dimensional EBH model, their expectation value varies with the strength of the dipolar moment, as summarized in Fig. 3, which reproduces the behaviour reported in Ref. Sowinski . For completeness, we mention that the one-dimensional EBH model can also exhibit a topological phase, denoted by Haldane-insulator phase, which is incompressible and characterised by $S_{x}(\pi/a)=0$ Batrouni2013 . The so-called string-order operator $O_{s}(|j-l|)$ signals its appearance Berg ; Batrouni2013 ; Batrouni2014 . We will omit to analyse its expectation value for the small system sizes we consider, since a non-vanishing expectation value is not meaningful.

In addition to this set of observables, we also consider the structure form factor at wave number $q_{y}=\pi/a$, which signals the onset of zigzag order and is defined as:

[TABLE]

where

[TABLE]

and

[TABLE]

is a real matrix whose elements depend on the physical parameters Silvi . When $S_{y}(\pi/a)\neq 0$, the dipoles form a zigzag transverse structure.

III.2 Phase diagrams

We now report phase diagrams for the salient properties of Hamiltonian $H_{BH}$, Eq. (10), evaluated by means of exact diagonalization on a lattice with periodic boundary conditions along $x$ and composed of 4 to 12 sites. The considered number of sites for a given phase diagram depends on the number of transverse modes we need to take into account in order to warrant the convergence of the calculations. In what follows we use the parameters of 85Rb-133Cs bosonic molecules with electric dipole moment of $p_{0}=1.15$ Debye Naegerl , confined by an optical lattice along $x$ at the interparticle distance $a=395$ nm, corresponding to half wavelength of the standing-wave laser, unless otherwise stated. The parameters $a$ and $p_{0}$, moreover, are the units of length and of the dipole moment we will refer to.

III.2.1 Phase diagram of the quasi one-dimensional array

We first consider the limit in which the trapping frequency $\omega_{y}>\omega_{y}^{(c)}$, so that the transverse motion is in the ground state of the transverse oscillator and the model is reduced to a single-band EBH model, described by Hamiltonian $H^{x}_{0}$, corresponding to Eq. (11) with $m=0$. We are interested first in reproducing the results of Ref. Sowinski with our multi-mode EBH model and therefore need to identify the conditions on the trap frequency $\omega_{y}$ for which we reproduce the single- and two-particle correlations and the component of structure form factor $S_{x}(\pi/a)$ as a function of $p$, when we fix the other corresponding parameters. For each value of $p$ we choose $\omega_{y}$ such that the width of the lowest eigenfunction $\phi_{0}(y)$ is equal to the width of the Wannier functions for the given lattice depth, $\int dx\,x^{2}w^{2}(x)=\int dy\,y^{2}\phi_{0}^{2}(y)$. The inequality $\omega_{y}>\omega_{y}^{(c)}$ is fulfilled for $p<3$ Debye.

Figure 3 displays $\phi$, $\Phi$, and $S_{x}(\pi/a)$ as a function of the strength of the dipole moment $p$, in units of $p_{0}$, for 12 sites and at half-filling. The results reproduce the behaviour reported in Ref. Sowinski . In order to highlight the role of pair and density dependent tunneling, in all figures we also give the value obtained by setting $T^{x}=P^{x}=0$. This comparison shows that these terms are essential for the appearance of two-particle correlations, signaling pair superfluidity. This occurs at sufficiently large value of $p$, which in turn scales the corresponding coefficients $T^{x}$ and $P^{x}$.

We now extend this analysis to unit fillings, $\Braket{n_{j}}=1$. This regime was not considered in Ref. Sowinski since the role of pair superfluidity is expected to be small, nevertheless it is relevant for studying the linear-zigzag transition. In order to benchmark this case, we determine the phase diagram in the limit where the lowest transverse band is occupied. Figures 4a-4c display the contour plots of the single particle correlation $\phi$, the two particle correlation $\Phi$ and the structure factor $S_{x}(\pi/a)$, respectively, as a function of the steepness of the confinement along the $z$-axis, $\sigma_{z}$, and of the strength of the dipole moment for the ground state of a lattice composed by 10 sites and 10 particles. The red-coloured region indicates the unstable regime, where the on-site interaction coefficient $U_{0}^{x}$ becomes negative: The border of this region is the line where $U_{0}^{x}=0$. Close to the unstable area, at small values of $U_{0}^{x}$ there is a striped region where the single and two particle correlations, $\phi$ and $\Phi$, vanish. We verified that the local compressibility also vanishes. In the same parameter area, the structure factor $S_{x}(\pi/a)$ is different from zero, see (c); we hence conjecture that the system is in a CDW phase. This conjecture is further supported by the exponential decay of the long range correlations $\phi_{l}=\langle a_{j}^{\dagger}a_{j+l}\rangle$. Outside of this region, in the stable regime and for $p$ sufficiently large the phase is characterized by an exponential decay of $\phi_{l}$ and by $S_{x}(\pi/a)=0$. We identify it with a MI phase. At small values of $p$ as well as close to the border separating with the CDW phase the system is SF. We do not find signatures indicating a PSF phase in the parameter regime we explored: $\Phi$ is different from zero in the region where $\phi\neq 0$ and scales as in a SF phase.

The peculiarity of these results can be better highlighted by considering the phase of the system as a function of $p$ at fixed values of $\sigma_{z}$. For sufficiently small values of $\sigma_{z}$, by increasing $p$ we observe a transition from SF to MI. At sufficiently large value of $\sigma_{z}$ increasing $p$ leads to a transition from SF to a CDW, before the system becomes unstable. Between these two regimes, there seem to be a small interval of values $\sigma_{z}$ where the system goes from SF to MI to SF by increasing $p$.

Figures 5a-5c show the behaviour of $\phi$, $\Phi$ and $S(\pi)$ when varying the strength of the onsite interaction while keeping $p$ constant. The behaviour reported in these plots can be put in direct connection with the ones of Fig. 4 since increasing $g$ partly corresponds to effectively decreasing $p$. Here, we clearly observe that CDW and MI are separated by a discontinuity in the structure form factor and in the single particle correlations, which occurs at the same value of $\sigma_{z}$. In order to determine the properties at the discontinuity we calculated the susceptibility of the ground-state fidelity $\mathcal{F}(\sigma_{z})=|\Braket{\Psi(\sigma_{z})}{\Psi(\sigma_{z}+\delta)}|$ , defined as Fidelity

[TABLE]

The susceptibility is different from zero at the point where the structure factor exhibits a discontinuity. We verified that its value increases with the particle numbers. On this basis we conjecture that this discontinuity signals a quantum phase transition.

III.2.2 The multi-mode model at the linear-zigzag instability

We now report properties of our model at the linear-zigzag instability and for unit filling. Including extra bands is here necessary but it severely limits the computational capability of exact diagonalization, since the dimensionality of the problem rapidly scales up with the number of orbitals. We first check how many states of the local basis shall be considered in order to warrant the convergence of the calculations for $\omega_{y}\sim\omega_{y}^{(c)}$. Figure 6 displays the occupation of the lowest four orbitals ($m=0,1,2,3$) as a function of $\omega_{y}$ for $N=6$ particles in a relatively shallow optical lattice, $V_{L}=6E_{R}$. For $\sigma_{z}=0.1a$ and for the considered values of the trap frequencies $\omega_{y}>\omega_{y}^{(c)}$ we find that only the lowest orbital is relevant, while for $\omega_{y}<\omega_{y}^{(c)}$, also the second orbital is occupied. Recalling that first orbital is even and the second orbital is odd, this change of occupation corresponds to the onset of the zigzag phase. In both cases, 99 % of the population is in the lowest two bands. This result remarkably shows that the basis decomposition we perform warrants a fast convergence even at transverse trap frequencies well below the mean-field critical value.

Figure 7 shows the zigzag order parameter $\xi=S_{y}(\pi/a)$, Eq. (28), as a function of $\omega_{y}$ for $N=4,6$ and a steeper optical lattice: $\xi$ increases by decreasing $\omega_{y}$, allowing to identify the zigzag phase. The intersection between the two curves at $N=4$ and $N=6$ suggests the location of the critical point, which is at a smaller value than the mean-field prediction and consistent with the DMRG result of Ref. Silvi .

IV Conclusions

In this work we have derived a multimode EBH model which can naturally describe the effects of the quantum fluctuations of an array of dipoles at the structural transition to zigzag order. Our model takes into full account the three-dimensional, anisotropic nature of the dipolar interaction. Our results show that the frequency of the transverse confinement controls not only the onset of zigzag order, but also determines the quantum phases of the molecules along the chain. The interplay between classical and quantum effects as a function of the transverse confinement is an open question, which will be addressed in future works performing numerical simulations with large numbers of particles. For this purpose the study here presented provides an important benchmark. Moreover, our model could be extended to describe structural transitions of cold polar molecules in arrays of one-dimensional tubes, in the setup analysed in Refs. Kollath ; Knap .

Acknowledgements.

The authors acknowledge discussions with Efrat Shimshoni, André Winter, and Pietro Silvi. They are especially grateful to Rebecca Kraus for the critical reading of this manuscript. Financial support by the German Research Foundation (DFG, GiRyd Priority Programme 1929 ”Giant Interaction in Rydberg Systems”) is gratefully acknowledged.

Appendix A Determination of the Bose-Hubbard coefficients

In this Appendix we derive the effective dipole-dipole interaction in two dimensions by integrating out the motion along the $z$-axis in the integrals needed to evaluate the coefficients for Eqs. (11), (18), and (22). In order to illustrate the procedure we first write these terms in generic form as

[TABLE]

where ${\bm{\rho}}=(x,y)$ and $V_{2D}(\bm{\rho}-{\bm{\rho}^{\prime}})$ contains the integrals in the $z,z^{\prime}$ variables and specifically takes the form:

[TABLE]

Thus, $V_{2D}({\bm{\rho}}-{\bm{\rho}^{\prime}})$ is an effective dipole-dipole interaction in two dimensions. Its form can be simplified by using center-of-mass $Z=(z_{1}+z_{2})/2$ and relative variable $z=z_{1}-z_{2}$. After integrating out the center-of-mass variable, Eq. (31) reads

[TABLE]

where $x=x_{1}-x_{2}$ and $y=y_{1}-y_{2}$. We then determine the integral in Eq. (30) using a convolution method Wall2013 . This consists first in writing the integral in Eq. (30) as

[TABLE]

where we dropped the indices for convenience and introduced $C({\bm{\rho}})$, which is defined as

[TABLE]

with $\mathcal{F}_{\bm{k}}$ the Fourier transform in two dimensions (${\bm{k}}=(k_{x},k_{y})$) and $\mathcal{F}_{\bm{k}}^{-1}$ its inverse. This procedure allows one to calculate the integral by computing a 2D Fourier transform and a 2D integral, instead of integrating a four dimensional integral in real space, thus saving computing time and allowing one to use a finer grid of discretization.

The Fourier transform $\mathcal{F}_{\bm{k}}\left[V_{2D}({\bm{\rho}})\right]$ can be explicitly calculated. We first use the definition of the inverse Fourier transform:

[TABLE]

We further observe that Eq. (32) can be rewritten as

[TABLE]

where $\tilde{A}(k_{z})={\rm e}^{-k_{z}^{2}\sigma_{z}^{2}/2}$ and

[TABLE]

Comparing Eq. (35) and Eq. (36) leads to identity and

[TABLE]

which can be analytically evaluated. It results that $\tilde{V}_{2D}({\bm{k}})=\tilde{V}_{2D}(|{\bm{k}}|)$, and in detail

[TABLE]

where we used $q=\sqrt{k_{x}^{2}+k_{y}^{2}}$ and $\text{erfcx}(x)$ is the scaled complementary error function: $\text{erfcx}(x)=e^{x^{2}}\text{erfc}(x)$ Abramowitz . The expression in Eq. (40) is identical to the one in Ref. Babadi , except for the constant term, which modifies the on-site interaction Erratum . In real-space it reads

[TABLE]

where $K_{0,1}(x)$ are modified Bessel function of second kind Abramowitz and $\rho=|{\bm{\rho}}|$. We note that the last term in Eq. (41) is an effective attractive contact interaction that can substantially modify the onsite coefficient of the multi-mode EBH model.