Dimerized Mott insulators in hexagonal optical lattices

We study bosonic atoms in optical honeycomb lattices with anisotropic tunneling and find dimerized Mott insulator phases with fractional filling. These incompressible insulating phases are characterized by an interaction-driven localization of particles in respect to the individual dimers and large local particle-number fluctuations within the dimers. We calculate the ground-state phase diagrams and the excitation spectra using an accurate cluster mean-field method. The cluster treatment enables us to probe the fundamental excitations of the dimerized Mott insulator where the excitation gap is dominated by the intra-dimer tunneling amplitude. This allows the distinction from normal Mott insulating phases gapped by the on-site interaction. In addition, we present analytical results for the phase diagram derived by a higher-order strong-coupling perturbative expansion approach. By computing finite lattices with large diameters the influence of a harmonic confinement is discussed in detail. It is shown that a large fraction of atoms forms the dimerized Mott insulator under experimental conditions. The necessary anisotropic tunneling can be realized either by periodic driving of the optical lattice or by engineering directly a dimerized lattice potential. The dimers can be mapped to to their antisymmetric states creating a lattice with coupled p-orbitals.

Superlattices with different site offsets and barrier heights can be created by superposing lattices with different wavelengths [4-8, 14-16, 28-31] or by employing polarizationdependent light potentials [11][12][13]17]. The unit cells of these lattices consist of multiple sites, adding a new degree of complexity to the system. If the sites within the unit cell have different energy offsets, normal Mott insulators with a population imbalance emerge [17,[28][29][30]32], where the atom are localized on individual sites. For the case without site-offsets, anisotropic tunneling couplings lead to non-trivial insulator phases with fractional filling in between the conventional Mott insulator phases [33][34][35][36][37]. Here, the particles localize on individual unit cells with a vanishing superfluid order parameter but are still delocalized within each unit cell. Originally, respective experiments have been proposed for one-dimensional superlattices, but so far, these phases with fractional filling have not been observed. For this the reason is two-fold. First, the phase diagram suggests a very small fraction of atoms to occupy the dimerized phase in a confined system. Second, a clear signature for discriminating the fractionally filled phase from the conventional Mott insulator was missing.
In this work, we theoretically study bosonic atoms in honeycomb optical lattices with adjustable tunneling matrix elements. The latter is achieved either by a periodic driving of a honeycomb lattice or by engineering directly a dimerized potential [14][15][16]. In addition to normal Mott insulating phases with integer filling, we find insulating phases with half-integer filling where the particles are delocalized on dimers. The dimers are naturally defined by the biatomic unit cell of the honeycomb lattice, allowing for non-trivial fractional-filling phases [33][34][35][36][37]. The phase diagram is studied by means of the strong-coupling expansion approach similar to Ref. [34] as well as by the cluster Gutzwiller mean-field method giving accurate results for honeycomb lattices [38]. This method grants the great advantage of the access to the excitation spectrum. The sites within a dimer are coupled by the tunneling matrix element J1, while the coupling between neighboring dimers is associated with the smaller tunneling J2. For J2 J1 a quasi-quadratic lattice structure of dimers is formed. (b) In the Mott insulator phase with integer filling ρ = 1 each atom is localized at a lattice site. (c) In the dimerized Mott insulator phase with ρ = 1/2, 3/2, ... the atoms are still delocalized within individual dimers, while the superfluid order parameter of the lattice vanishes. (d) By lattice shaking the tunneling matrix elements along the dimers can be tuned to negative values, where each dimer wave function resembles a porbital state. In the superfluid state a stripe order is formed.
We show that the characteristic local excitations allow distinguishing experimentally between the conventional and this dimerized Mott insulator state. The excitation gap gives an estimate of the required temperatures for observing this quantum phase. Furthermore, we simulate two-dimensional lattice planes with harmonic confinement using realistic experimental parameters. We find that the dimerized state is formed by a large fraction of atoms (> 70%) and thus is well observable in the proposed experiment.

I. PHASE DIAGRAM
There are two perspectives to realize the dimerized Mott insulator phase in experiments with hexagonal lattices. First, we propose to use an optical honeycomb lattice generated by three running laser beams as in Refs. [11][12][13]. A honeycomb lattice features a two-atomic unit cell, where the intra-cell bond and inter-cell bonds have different orientations, which allows addressing them independently by lattice shaking techniques. By modulating the relative phases of the beams, the lattice is periodically accelerated on an elliptical orbit [23,25]. Employing Floquet theory, one can obtain an effective timeaveraged Hamiltonian, where the tunneling matrix elements are modified by a Bessel function depending on the driving parameters [23][24][25][26][27]. This allows engineering two different and tunable tunneling matrix elements J 1 and J 2 in the vertical and horizontal direction (see Fig. 1a). When J 1 is larger than J 2 , the honeycomb lattice separates into a dimerized square lattice of double-wells coupled by the reduced matrix element J 2 . Second, a dimerized lattice can also be obtained in the setup [14][15][16], where the dimerized honeycomb lattice sketched in Fig. 1 is created due to the interference of two collinear phase-shifted laser beams. Here, an additional shaking is not required but the control of the tunneling parameters is not independent. In the first case it is, e.g., possible to achieve negative tunneling matrix elements.
The phase diagram for this setup is shown in Fig. 2 for J 1 = 10J 2 as a function of the chemical potential µ and the tunneling energy J 1 in units of the repulsive on-site interaction U . Both the Mott insulator phases (MI) and the dimerized Mott insulator phases (DMI) are characterized by a vanishing superfluid order parameter ψ = 0. In the remaining regions of the phase diagram with ψ = 0 the particles are in a superfluid state. In the Mott insulator phase, the particles are localized at individual lattice sites as depicted in Fig. 1b. In contrast, the dimerized Mott insulator phase with fractional filling ρ = 1/2, 3/2, ... is characterized by a delocalization of particles within the dimers (Fig. 1c). In the limit of fully separated dimers (J 2 → 0), the local ground state on a dimer for filling ρ = 1/2 reads |s = 1 √ 2 (|1, 0 +|0, 1 ), where |n L , n R denotes the occupation of left and right dimer sites. In the dimerized insulator phase the local particle-number fluctuations (∆n) 2 are large, whereas they are strongly suppressed in a conventional Mott insulator, which can be used to distinguish between the phases. Thus, the single-site resolution available in a quantum-gas microscope experiment would reveal a random distribution on individual lattice sites but a fixed integer occupation when summing up both dimer sites. When the tunneling matrix element J 1 is negative, which can be achieved with the lattice shaking technique, the antisymmetric state |a = 1 √ 2 (|1, 0 −|0, 1 ) becomes the dimer ground state resembling a p-orbital. Hence, we can understand the setup as a square lattice of p-orbitals. In the superfluid phase, a negative tunneling matrix element leads to an alternating sign of the superfluid order parameter, which can be mapped onto a classical spin model and is therefore referred to as antiferromagnetic coupling. The phases align due to the positive dimer coupling J 2 > 0 such that the sign of the dimer wave function is the same along the respective bonds ( Fig. 1d), which minimizes the tunneling energy. This alignment leads to a stripe order of the superfluid order parameter. In the Mott state where the superfluid order parameter vanishes, this alignment persists in the nearest-neighbor correlations. A suitable gauge transformation maps the symmetric (ferromagnetic) onto the antisymmetric ground state which is not frustrated. Therefore, we restrict ourselves to the symmetric case for J 1 > 0 in the following. We should stress however that the aforementioned equivalence means that the dimerized lattice resembles p-orbital-like physics even without antiferromagnetic driving of the lattice.
The solid lines in the phase diagram in Fig. 2 are computed using a cluster mean-field approach with a cluster size of 12 sites, as described in Sec. V. Due to the structure of the dimerized state, conventional mean-field theory is not able to capture the dimerized Mott insulator phases. As a second approach, we apply strong-coupling perturbative expansion to third order (dashed lines) which allows analytical results for the phase diagram. This approach is detailed below in Sec. IV, where we also give the explicit expression for the second-order perturbation (dotted lines) as a reasonable approximation. Both methods use the tight-binding Bose-Hubbard Hamiltonian Here, the brackets i, j denote sites i and j on the same dimer connected via J 1 , while sites i, j are nearest-neighbors sites connected via J 2 on different dimers. The repulsive onsite interaction is denoted by U with the particle number operatorn i =b † ib i .

II. EXCITATION SPECTRUM
Due to the internal structure of the dimerized Mott insulator state, its fundamental excitations differ strongly from the normal Mott insulator phases. Therefore the excitation spectrum allows distinguishing between the two insulating phases. In the Mott insulator, the lowest excitation is the creation of a particle-hole pair resulting in an empty and a doubly occupied site. The corresponding excitation energy is the additional onsite interaction U . The excitation spectrum of the fractional insulator with ρ = 1/2 is not U -gaped since each dimer offers an empty site, where the particle from a neighboring dimer can tunnel to. Even for dimerized Mott insulators with higher fillings the U -gap vanishes.
In Fig. 3c the ground-state configuration of the dimerized Mott insulator is depicted for two unit cells, each being populated by the symmetric state |s . There are two different fundamental excitations both on the order of 2J 1 . First, the particle-hole excitation E ph , where one particle is excited by hopping to a neighboring dimer as depicted in Fig. 3b. The excitation energy corresponds to the loss of delocalization energy J 1 within the empty dimer reduced by the interaction energy on the doubly occupied dimer. Second, a particle can be excited within the same dimer from the symmetric ground state |s to the antisymmetric state |a associated with the energy E as = 2J 1 (see Fig. 3a).

III. HARMONIC CONFINEMENT
The excitation spectrum as a function of J 2 /J 1 is depicted in Fig. 3d in units of the tunneling energy J 1 for J 1 = 0.2U . The black markers represent the numerical data calculated with the cluster mean-field approach using 16 sites within the insulating phase with ψ = 0. The shaded areas indicate the results of first-order strong-coupling perturbation theory. The antisymmetric excitation E as is not broadened within first-order perturbation and is only affected by higherorder processes. In the limit J 2 → 0, a particle-hole excitation for an insulating phase with n particles per dimer has the energy E ph = E n+1 + E n−1 − 2E n . Thus, the particle-hole energy for the fractional insulator with n = 1 For finite J 2 , first order perturbation (blue shaded area) leads to a delocalization of the particle and the hole. The minimum and maximum of the emerging band are given by where j 2 n + j 2 n+1 ≈ 1 for U J 1 (see Sec. IV for details). The excitation gap E ph ≈ J 1 gives an estimate of the required temperature for the dimerized phase T ≈ E ph /k B , where k B is the Boltzmann constant. Assuming the experimental parameters given in section III and V = 6.5 E rec we obtain T ≈ 20 nK. In general, this value increases with J 1 and decreases with the ratio J 2 /J 1 . Alternatively, a larger value of the scattering length a s /a, where a is the lattice constant, allows for a larger value of J 1 , increasing the excitation gap.
The higher excitations in the energy spectrum are combinations of the two fundamental excitations described above, i.e. two particle-hole excitations 2E ph , two asymmetric excitations 2E as , and a combination of both E ph + E as . The numerical spectrum is distorted by higher-order tunneling processes and the interaction between individual excitations. Due to the finite size of the cluster the band width is reduced, which is in particular noticeable for two particle-hole excitations due to the limited possibilities to delocalize. The U -excitation at E ≈ 5J 1 corresponding to a doubly occupied site lies in the continuous part of the spectrum for most parameters.
In an experiment, the lattice modulation technique proposed in Ref. [39] can be applied for the direct observation of the dimerized Mott insulator phase by addressing these fundamental excitations. The excitation gap on the order of 2J 1 is therefore the characteristic signature of the fractional insulator and could serve as prove for its experimental realization.
In this section, the question is addressed whether the dimerized Mott insulator in the honeycomb lattice can be observed in an experimental setup, where the optical lattice is superimposed by a harmonic confinement, leading to a spatially varying chemical potential µ(r). The relatively small extend of the dimerized Mott insulator phase in the phase diagram might suggest that a dominating part of the atoms are in the superfluid or in the Mott insulator with ρ = 1 coexisting with the dimerized phase. The cluster mean-field approach (see Sec. V) allows the simulation of a two-dimensional lattice of realistic size by iteratively moving the cluster through the lattice [38,40]. This introduces a site-dependent mean-field, where at every iteration the local order parameter is updated until the results converge. If the ratio J 1 /U is considerably smaller than the tip position of the dimerized insulator phase, the results are influenced only to a minor degree by the cluster size. Using six-site clusters, we can determine the extent of the phases in a lattice with harmonic confinement accurately. In Fig. 4 the results are shown for J 1 = 0.1U , J 2 = 0.1J 1 and a chemical potential µ c in the center of the trap. As an example, for 87 Rb with a scattering length of a s ≈ 100 a 0 , this corresponds to a lattice depth of V = 9.5 E rec [12,17], where a 0 is the Bohr radius, E rec = 2 k 2 2m is the recoil energy, k = 2π λ is the wave vector of the lattice beams with a wavelength of λ = 830 nm and m is the atomic mass of 87 Rb. For each lattice site we apply the local density-approximation µ(r i ) = µ c −V trap (r i ) with a harmonic trapping potential V trap and trap frequencies of 32 Hz (a-c) to 35 Hz (d-f).
For µ c = 0.02U , only the dimerized insulator persists (Fig. 4a-c), whereas for µ c = 0.1U we observe the coexistence of normal and dimerized Mott insulator (Fig. 4d-f). The superfluid order parameter ψ shown in Fig. 4b,e vanishes in the insulating phases and increases to a finite value in between forming superfluid rings. For a threshold of Ψ < 0.1, we find that the dimerized Mott insulator phase is occupied by 490 of a total of 680 atoms (a-c) and 400 of 1160 atoms (d-f), respectively. The density in Fig. 4a,d shows the typical wedding cake structure but with half-integer steps. In the dimerized Mott insulator phase with ρ = 1/2, a periodic density modulation along the dimer axis appears. While the average density on a dimer is fixed to ρ = 1/2, the density on the two dimer sites adjust itself according to the gradient of the chemical potential along the dimers. The persistence of the insulating phase despite this strong impact of the gradient illustrates its robustness.
A cut through the trap along the dimer axes is shown in Fig. 4c, f. The density profile (red lines) and the superfluid order parameter (blue lines) clearly indicate the Mott insulator and dimerized Mott insulator regime. This agrees well with the expectation from the phase diagram in Fig. 2 for infinite size (for µ < µ c ). The dimerized Mott insulator plateaus in the density profile show an oscillating behavior along the dimer axis which is caused by the trap as discussed above. The green lines represent the total particle number and indicate a comparatively large occupation of the dimerized phase in the case of Fig. 4f.
The local particle-number fluctuations (∆n) 2 are shown as black lines and demonstrate the expected large value of (∆n) 2 = 0.25 in the dimerized insulator phase. The large particle-number fluctuations in combination with the vanishing order parameter characterize the dimerized Mott insulator phase, whereas in the conventional Mott insulator all fluctuations are suppressed.

IV. PERTURBATION THEORY
The phase diagram Fig. 2 as well as the excitation spectrum Fig. 3d can be approximated with the strong-coupling perturbative expansion technique [33,34,[41][42][43]. For this we recast the full Hamiltonian (1) to an effective model with dimer unit cells. As a first step we find the eigenstates |n of n particles in a dimer with energies E n . In contrast to Ref. [34], we can restrict the calculations to the respective symmetric ground states, due to the large ratio J 1 /J 2 . This simplifies the approach significantly and allows us to include perturbations up to third order. The energies for the lowest values of . This approximation is valid as long as the perturbation zJ 2 with the coordination number z = 4 is much smaller than the energy E (1) n − E (0) n , i.e., zJ 2 U, J 1 . The coupling between neighboring dimers is given by the operator where the brackets i, j label neighboring dimers i, j instead of sites and the operatorsd i (d † i ) annihilate (create) a particle on a dimer. More precisely, for the annihilation of a particle on a dimer the operator is given by the projection on the nparticle ground stateŝ whereb acts on one of the two equivalent sites of the dimer. As a compact notation we define the coupling parameter j n = n − 1|b|n . For a conventional single-site lattice model it is j n = √ n, whereas here j n depends on the exact form of the ground states that are functions of J 1 /U . With the above restriction we obtain an effective lattice model Within this dimerized model, all insulating phases are treated on the same footing as product states of dimer ground states |n . Odd fillings n correspond to dimerized insulators and even n to conventional Mott insulators. We now discuss the first lowest three orders of perturbation by inter-dimer tunneling J 2 . In the unperturbed case of J 2 = 0, the creation of a particle (hole) excitation is associated with the energy with the chemical potential µ. In this case, only insulating phases exist and their phase boundaries µ (0) p|h can be obtained from the condition E p|h = 0, where the chemical potential compensates for the energy of one additional or missing particle per dimer.
In first order perturbation, the delocalization of a particle or a hole over the lattice results in lower and upper energy bounds Thus, the energy band for a particle-hole excitations lies be- h± for one and 2E ph± for two excitations indicated in figure Fig. 3d (colored areas). In addition, the asymmetric excitation that goes beyond the effective model (5) is indicated as unbroadened line.
The coupling between the dimers gives rise to superfluid phases surrounding the insulating phases. The phase boundaries shown in Fig. 2 are obtained from the minimum energy of the excitations. In first order perturbation they read For the determination of second-and third-order energy, one has to include processes with amplitudes on the order of J 2 2 and J 3 2 . The phase boundaries in second-order perturbation read The first additional term in (12) accounts for the second-order energy of the insulator that is obtained by processes via a virtual particle-hole excitations. When a particle excitation is present, these bidirectional processes are inhibited on neighboring bonds leading to a factor of 2z = 8. They are partly substituted by processes via the intermediate state where a particle tunnels onto the excitation, which is captured by the second additional term. Analogously, the energy of a holeexcitation relative to the insulator accounts for all possible second-order processes. We compute the phase diagram up to third order and find a good agreement with the cluster mean-field approach detailed below. The pointy tip of the insulator lobes is due to the vanishing energy of particle-hole excitations at the crossing point of the phase boundaries. Here, the perturbation series cannot be limited to finite-order processes.

V. CLUSTER GUTZWILLER THEORY
Our numerical simulations are performed using a cluster mean-field approach [33,38,40,[44][45][46][47]. The cluster approach allows capturing phases with strong short-range correlations such as fractional insulators formed on unit cells covering more than one site. Conventional single-site mean-field approaches such as the Gutzwiller approach are not capable of finding the dimerized insulator phases. Furthermore the cluster mean-field approach is well suited to obtain results for inhomogeneous systems as discussed in Sec. III and gives detailed information on the local excitation spectrum presented in Sec. II.
The cluster Gutzwiller approach decouples a cluster of sites from the rest of the lattice and couples it to a mean-field at its borders. The latter is determined from the exact diagonalization of the cluster and is updated in an iterative process. This allows taking into account local correlations exactly and thereby gives far more precise results than conventional meanfield methods. The improvement is especially pronounced for lattices with a small number of nearest neighbors, such as hexagonal lattices. In the many-particle cluster basis |N the Hamiltonian matrix decomposes in two parts describing the cluster according to Eq. (1) and its boundary. The HamiltonianĤ boundary describes the coupling of sites at the boundary σ of the cluster to sites outside the cluster and readŝ where ν i is the number of mean-field bonds at site i. The superfluid order parameter ψ = b at the boundary is obtained from the innermost site in the cluster. When at least one of the two tunneling matrix elements is negative, the order parameter ψ shows an alternating sign as depicted in Fig. 1d. Note that we apply periodic boundary conditions along one direction of the clusters, which increases the ratio of inner cluster bonds to mean-field bonds. We use clusters of up to 16 sites and restrict the basis |N further using cut-offs on the fluctuations and the number of particles per site (for further details see Ref. [38]). We carefully checked for convergence. The excitation spectrum obtained from the cluster meanfield approach agrees well with the perturbation theory. The reduced band width of the band of two particle-hole excitations is due to the finite size of the cluster, as the two excitations already spread over four dimers on a eight-dimer cluster, which limits the possibilities for delocalization.
The phase diagram matches the predictions of the perturbation theory well for all phases (0 ≤ n ≤ 4). At the tips of the lobes, the perturbation theory is not expected to give precise results due to the vanishing energy of particle-hole excitations. We expect further deviations due to the restriction to the symmetric ground states on each dimer in the perturbation theory approach. The cluster mean-field approach on the other hand is restricted to a finite cluster size. The good agreement between the two methods indicates that the approximations are justified and the true phase boundary is well approximated.

VI. CONCLUSIONS
In conclusion, we have shown that when anisotropic tunneling is introduced to optical honeycomb lattices, dimerized Mott insulator phases with fractional fillings appear. In these phases, the superfluid order parameter vanishes but large particle number fluctuations persist on the individual lattice sites. We have calculated the phase diagram using two different approaches, namely a cluster Gutzwiller approach and the strong-coupling perturbation expansion technique, and found excellent agreement. The former method allows us to study the excitation spectrum, which allows the distinction of normal and dimerized Mott insulator. In a possible experiment with a harmonic confinement the dimerized insulator is formed by a large fraction of the atoms and should therefore be observable. Therefore, optical honeycomb lattice experiments should be well suited to realize and probe the proposed dimerized Mott insulator phase, which to our best knowledge has not been measured experimentally so far. When driving an optical honeycomb lattice the intra-dimer tunneling coupling can be tuned negative to realize p-orbital like physics.
After submission, we became aware of Ref. [48] discussing a similar setup.