Pair hopping in systems of strongly interacting hard-core bosons

TL;DR

This study uses quantum Monte Carlo simulations to explore how pair-hopping processes influence the phase diagram of strongly interacting hard-core bosons on a square lattice, revealing their role in destabilizing Mott phases.

Contribution

It introduces the effect of pair hopping in the hard-core boson model and demonstrates its impact on the stability of Mott phases, a novel aspect not previously analyzed.

Findings

Pair hopping destabilizes Mott phases at certain fillings.

Different Mott phases melt at different pair-hopping amplitudes.

Pair hopping influences the competition between superfluid, insulating, and supersolid phases.

Abstract

We have used the Stochastic Series Expansion quantum Monte Carlo method to study interacting hard-core bosons on the square lattice, with pair-hopping processes supplementing the standard single-particle hopping. Such pair hopping arises in effective models for frustrated quantum magnets. Our goal is to investigate the effects of the pair hopping process on the commonly observed superfluid, insulating (Mott), and super-solid ground-state phases in the standard hard-core boson model with various interaction terms. The model is specifically motivated by the observation of finite dispersion of 2-magnon bound states in neutron diffraction experiments SrCu$_2$(BO$_3$)$_2$. Our results show that the pair hopping has different effects on Mott phases at different filling fractions, "melting" them at different critical pair-hopping amplitudes. Thus, it appears that pair hopping may have an important role in determining which out of a potentially large number of Mott phases (stabilized by details of the charge-diagonal interaction terms) actually survive the totality of quantum fluctuations present.

Equations

$$\begin{array}{rl}\mathcal{H}=& -t\sum _{⟨i,j⟩}(a_i^{†}{a}_j+h.c.)-t_p\sum _{\square }(a_i^{†}{a}_j^{†}{a}_k{a}_l+h.c.)\\ & +V\sum _{⟨i,j⟩}{n}_i{n}_j+V_d\sum _{⟨i,k⟩}{n}_i{n}_k-\mu \sum _i{n}_i\end{array}$$

$$⟨n⟩=\frac{1}{N}\sum _i{n}_i.$$

$${\rho }_s=\frac{{\partial }^{2}f(\varphi )}{\partial {\varphi }^2},$$

$${\rho }_s=({\omega }_x^2+{\omega }_y^2)/\beta ,$$

$${\omega }_{\alpha }=(N_{\alpha }^{+}-N_{\alpha }^{-})/L,(\alpha =x,y).$$

$${\rho }_t=({\omega }_{t,x}^2+{\omega }_{t,y}^2)/\beta ,$$

$${\rho }_{tp}=({\omega }_{tp,x}^2+{\omega }_{tp,y}^2)/\beta .$$

$$\begin{array}{rl}{\omega }_{t,\alpha }& =(N_{t,\alpha }^{+}-N_{t,\alpha }^{-})/L,\\ {\omega }_{tp,\alpha }& =(N_{tp,\alpha }^{+}-N_{tp,\alpha }^{-})/L,\end{array}$$

$${\omega }_{\alpha }={\omega }_{t,\alpha }+{\omega }_{tp,\alpha },(\alpha =x,y)$$

$${\rho }_s\ne {\rho }_t+{\rho }_{tp}.$$

$$S(\vec{k})=\frac{1}{N^2}\sum _{\vec{r}}{e}^{i\vec{k}\cdot \vec{r}}C(i,j),$$

$$C(i,j)=⟨n_i{n}_j⟩.$$

Full text

I Introduction

The interplay between competing interactions, enhanced quantum fluctuations due to reduced dimensionality and external fields in interacting lattice bosons result in a rich array of novel quantum states of matter that have been intensely studied over the past several decades Batrouni et al. (1995); Batrouni and Scalettar (2000); Frey and Balents (1997); Hébert et al. (2001); Sengupta et al. (2005); Jiang et al. (2012); Schmidt et al. (2006); Freericks and Monien (1994); Kühner and Monien (1998); Capogrosso-Sansone et al. (2008); Chen et al. (2008). In recent years, experimental advances have allowed the realization of these bosonic phases, such as the superfluid (SF), Bose-Einstein condensation (BEC), Mott insulator (MI) and density modulated crystalline phases with different ordering wave vectors, in a variety of physical systems such as optical lattices with ultracold atoms Jaksch et al. (1998); Greiner et al. (2002); Morsch and Oberthaler (2006); Islam et al. (2015), quantum magnets and excitons and polaritons Kasprzak et al. (2006); Snoke (2002) in semiconductor quantum wells. These systems are now opening up new frontiers in the study of strongly-correlated quantum many-body systems.

Quantum magnets, in particular, have long served as a versatile testbed for interacting lattice bosons in a controllable manner. The low-lying magnetic excitations, magnons, obey Bose-Einstein statistics and are an almost ideal realization of lattice bosons Giamarchi et al. (2008). The discovery of Bose-Einstein condensation in insulating magnets such as \ceTlCuCl_3 Nikuni et al. (2000); Rüegg et al. (2003), \ceBaCuSi_2O_6 Sasago et al. (1997); Rüegg et al. (2007); Jaime et al. (2004) and \ceNiCl_2-4SC(NH_2)_2 Yin et al. (2008); Paduan-Filho et al. (2009) heralded the search for novel quantum phases of interacting bosons in quantum magnets.Often, in quantum magnets, geometrical frustration induces quantum phases and phenomena that are not observed in their non-frustrated counterparts, e.g., dimensional reduction at a quantum critical point in \ceBaCuSi_2O_6 Jaime et al. (2004) and magnetization plateaus in \ceSrCu_2(BO_3)_2 on the geometrically frustrated Shastry-Sutherland lattice Miyahara and Ueda (1999); Kageyama et al. (1999); Koga and Kawakami (2000). Understanding the nature and mechanism of formation of the plateaus in \ceSrCu_2(BO_3)_2 has been the subject of intense experimental and theoretical studies during the past two decades Dorier et al. (2008); Corboz and Mila (2014); Miyahara and Ueda (2000). The ground state of the compound is comprised of orthogonal dimer singlets within the weakly coupled two-dimensional planes. In an external magnetic field, field-induced triplons constitute the lowest magnetic excitations. Theoretical modeling and neutron scattering experiments show that strong geometric frustration significantly suppresses the delocalization of triplons Kageyama et al. (2000) and prevents the onset of field-induced BEC of triplons that is commonly observed in other dimer quantum magnets such as \ceTlCuCl_3. The magnetization plateaus are understood as periodic arrangements of the triplons in regular patterns at commensurate fillings. However, the mechanism of triplon rearrangement into crystal orderings remain an open question.

Several different models have been proposed in the past to describe the magnetization profile of \ceSrCu_2(BO_3)_2, treating the field-induced triplons as hard-core bosons, with varying degrees of success. These include long range interactions and correlated nearest neighbor hopping of triplons, among others Miyahara and Ueda (2003); Momoi and Totsuka (2000); Dorier et al. (2008); Mila et al. (2008). In Ref. Kageyama et al. (2000), neutron scattering experiments performed by Kageyama et al. on \ceSrCu_2(BO3)_2 showed that while isolated triplons are localized, bound pairs of triplons exhibit pronounced dispersion, although the cost of pair formation is high. This may provide a potential mechanism for the rearrangement of the triplons into periodic patterns observed at the magnetization plateaus and has motivated us to explore the role of dynamically generated triplon pairs in modifying the field driven properties of interacting triplons. Our goal is not to derive an exact microscopic model of \ceSrCu_2(BO3)_2 and provide a quantitative explanation for the magnetization plateau formation therein. Instead, we want to isolate the effects of dispersive bound pairs of triplons in a generic quantum magnet with multiple competing interactions through the introduction of a new effective Hamiltonian, and investigate the dynamics that arise from such a lattice model.

In this paper, we study a system of interacting hard-core bosons with single and pair hopping with nearest-neighbor (nn) and next-nearest-neighbor (nnn) repulsions on a square lattice. Field-induced triplons on the dimers can be faithfully mapped on to hard-core bosons through the Matsubara-Matsuda transformation Matsubara and Matsuda (1956) and the dispersion of bound pairs of triplons translate to pair-hopping processes in the bosonic model, where a pair of hard-core bosons on nn sites hop together to the neighboring sites. While such processes occur within the standard framework of the canonical hard-core boson model with single-particle hopping, the amplitude of the effective process is small. Motivated by the experimental observation in \ceSrCu_2(BO3)_2 (suppressed single triplon dispersion and pronounced triplon pair dispersion), the relative magnitude of the pair hopping process is chosen to be large and considered as an independent term in the Hamiltonian. Our goal here is to explore the effects of finite pair hopping on the various ground state phases of the hard-core boson model, and to investigate if new many-body phases are engineered by the pair hopping process.

The paper is organized as follows. In Sec. II, we introduce our model and define the relevant parameters in the Hamiltonian. Section III describes how the pair hopping process can be incorporated into the Stochastic Series Expansion Quantum Monte Carlo scheme, which involves a straightforward generalization from 2-site bond operators to 4-site plaquette operators. Section IV defines the observables that are measured from the Monte Carlo simulation. Section V presents the main results of our numerical simulations, where we include illustrative phase diagrams for a wide range of Hamiltonian parameters as well as more detailed observable plots and accompanying analyses. Finally, in Sec. VI we discuss the significance of our results; namely how pair hopping modifies the formation and stability of the various Mott phases and its implications to our understanding of \ceSrCu_2(BO_3)_2.

II Model

The Hamiltonian for the model described above is given by

[TABLE]

where $a_{i}^{\dagger}$ and $a_{j}$ are the creation and annihilation operators respectively on sites $i$ and $j$. The $\square$ denotes a four-site plaquette on which our Hamiltonian parameters are defined, with sites we label $i,j,k,l$, as shown in Fig. 1. $n_{i}=a^{\dagger}_{i}a_{i}$ is the number operator at site $i$. $t$ and $t_{p}$ are the single and pair hopping amplitudes, respectively. $V$ and $V_{d}$ are the nn and nnn repulsion, respectively, and $\mu$ is the chemical potential. We work in the hard-core boson limit, i.e., the possible local occupancies are restricted to $n_{i}\in$ {0,1}. A square lattice with periodic boundary conditions of $N=L\times L$ sites is assumed. We set $t=1$ henceforth.

The nnn repulsion term $V_{d}$ in Eq. (1) is necessary for promoting pair formation. In Fig. 1(b), bosons on diagonal sites $\langle i,k\rangle$ incur an energy cost $+V_{d}$, which increases the likelihood of nearest-neighbour pairs occuring ($\langle i,j\rangle$ in Fig. 1(c)). The nn bosons are subsequently able to hop as pairs in proportion to the magnitude of $t_{p}$.

III Stochastic Series Expansion Method

We have used the Stochastic Series Expansion (SSE) Quantum Monte Carlo (QMC) Sandvik (1999); Sandvik and Kurkijärvi (1991) method to simulate the Hamiltonian (Eq. (1)) on finite-size systems. The SSE is a finite-temperature algorithm based on the stochastic evaluation of the diagonal matrix elements of the density matrix, $\exp(-\beta{\cal H})$, in a Taylor series expansion.

The SSE method employs the operator loop update method in sampling the configuration state space for the ground state configuration. The loop update involves the construction of a linked vertex list, where lattice sites are propagated in imaginary time, with diagonal and off-diagonal operators acting between the propagation levels according to a stored operator string. Sites connected by an operator between propagation levels are known as vertices. Configuration updates are achieved by the introduction of a ‘defect’ - a boson occupation inversion in the hard-core limit - into a random leg of a vertex. The defect is then propagated throughout the linked list, until the defect meets its initial introduction site and the loop is closed. The propagation of the defect is stochastically sampled in a manner proportional to the weights of the resulting vertices. After closing the loop, the lattice configuration and operator string are updated to reflect the changes made.

On 2D square lattices, the SSE loop update algorithm considers operators acting on 2-site bonds, such that vertices are 4-legged: 2 sites before and 2 sites after the action of an operator. To incorporate the pair hopping procedure, one needs to consider operators beyond 2-site bond operators. In particular, we consider operators that act on the 4-site plaquette mentioned in Fig. 1. This means that vertices in the linked list are now 8-legged: 4 sites before and 4 sites after the action of the operator. Conventional diagonal and single particle hopping operators carry over easily to the plaquette case. Our focus of the discussion will be on the incorporation of pair hopping operators in the loop update procedure.

We note that only slight modifications are required to achieve pair hopping in the SSE framework in the context of plaquette operators. Similar to the case of the hard-core boson model with ‘pair rotation’ of two bosons occupying diagonally opposite corners of a plaquette Sandvik et al. (2002), the pair hop operators are introduced to the operator string with the same linked list loop update procedure that is conventional in a 2-site bond operator SSE scheme. An illustration of the procedure is shown in Fig. 2. The crucial insight comes from the fact that pair hop operators are created only from single hop operators in the loop update. Figure 2 shows one way in which an existing single hopping operator in the operator string can be converted to a pair hopping operator through propagation of a single defect. Consequently, this implies that in our SSE framework, a non-zero $t$ in the Hamiltonian of Eq. (1) is necessary for the simulation to incorporate pair boson propagation.

IV Observables

In this section, we define the observables measured in our simulations that will be the basis of our analysis in Sec. V. The average boson density is defined as

[TABLE]

The total superfluid density (stiffness) is given by

[TABLE]

where $f(\phi)$ is the free-energy (or ground-state energy at $T=0$) density in the presence of a phase twist $\phi$. It is evaluated in SSE simulations as

[TABLE]

where $\omega$ is the winding number in the $x$ or $y$ directions, defined as

[TABLE]

$N_{\alpha}^{+}$ is the total number of particle hops in the arbitrarily chosen positive direction of the lattice. This implies a pair particle hop in the positive direction increments $N_{\alpha}^{+}$ by 2. On a conventional Bose-Hubbard model without pair hopping, it is identical to the total number of operators $a_{i}a^{\dagger}_{j}$ in the QMC operator string, if site $j$ is in the positive direction of site $i$.

To quantify the magnitudes of single and pair particle hopping separately, we define the single and pair particle stiffness, $\rho_{t}$ and $\rho_{tp}$, respectively, as

[TABLE]

and

[TABLE]

$\omega_{t,\alpha}$ and $\omega_{tp,\alpha}$ are the net sum of the single and pair particle hops, respectively, for $\alpha=x,y$. Concretely, we define them as

[TABLE]

where $N^{+}_{t,\alpha}$ is the total number of single particle hops in the positive direction, and $N^{+}_{tp,\alpha}$ is the total number of pair hops in the positive direction. A pair hop in the positive direction increments $N^{+}_{tp,\alpha}$ by 2, and vice-versa. From our definitions in Eq. (8),

[TABLE]

i.e., the total winding number is the sum of the single and pair winding numbers. Note that due to the way the various stiffness are defined,

[TABLE]

It should be noted that $\rho_{t}$ and $\rho_{tp}$ serve as useful quantities in measuring the relative contributions of single and pair currents in the system, but do not constitute experimentally measurable observables such as the total stiffness $\rho_{s}$ defined in Eq. (3).

In order to identify the presence of density modulation, or equivalently crystal ordering, we compute the static structure factor

[TABLE]

where $\vec{r}$ is the vector representing the separation of sites $\langle i,j\rangle$, and $\vec{k}=(k_{1},k_{2})$ is the wave vector, where $k_{1},k_{2}\in[0,2\pi]$. $C(i,j)$ is the density-density correlation function Sandvik (1997), defined as

[TABLE]

Simulations in this study are done at $\beta=L$ to extract ground-state properties, with simulated annealing Kirkpatrick et al. (1983) carried out at the equilibration step of the operator string to ensure convergence of the QMC simulation.

V Results

A representative ground-state phase diagram of the model, Eq. (1), in the parameter space of the pair-hopping amplitude, $t_{p}$, and the strength of next-nearest neighbor interaction, $V_{d}$, at fixed $t$ (single-particle hopping amplitude) and $\mu$ (the chemical potential) is shown in Fig. 3. The nearest neighbor interaction strength, $V$, is chosen as the unit of energy and the Hamiltonian parameters are expressed in units of $V$. In the limit of $t_{p}=0=V_{d}$, Eq. (1) reduces to the canonical Bose Hubbard model where the ground state for the chosen values of $t$ and $\mu$ is a checkerboard solid, with an ordering wave vector $\bf{k}=(\pi,\pi)$. The density of particles is constant at $\langle n\rangle=\frac{1}{2}$ and there is a gap to adding or removing a boson. As the strength of the next-nearest neighbor interaction is increased (at $t_{p}=0$), there is a transition to the superfluid phase at an intermediate value of $V_{d}/V$, where competing nn and nnn interactions suppress any crystallization of the bosons into a density wave. Eventually, for sufficiently strong nnn neighbor repulsion, the ground state enters a supersolid (SS) phase. The wave vector of the underlying solid order (density modulation of the bosons) changes to $(\pi,0)$, reflecting a striped solid. The density of particles deviates from $\langle n\rangle=\frac{1}{2}$ and the additional particles form a superfluid that co-exists with the solid ordering, resulting in a SS ground state. The pair-hopping process enhances the extent of superfluid phase at the cost of the solid orders, suppressing both the checkerboard solid and SS phases completely for sufficiently strong $t_{p}$. The SF phase has contributions from both single particle and pair currents – this is confirmed by finite values of the stiffness for both currents, viz., $\rho_{t}$ and $\rho_{tp}$. The pair current contribution is finite for any non-zero $t_{p}$, with the relative contribution increasing monotonically with $t_{p}$ (as shown by the color density profile in the phase diagram) 111It is important to note that although the Hamiltonian (Eq. (1)) conserves both the particle number and the number of pairs, in the simulations, only the total current is conserved in lieu of the individual $t$ and $t_{p}$ currents. Consequently, one gets $\rho_{t},\rho_{tp}>0$, even in the insulating phases, where the currents should be zero. This is due to Gaussian fluctuations of the individual currents regardless of winding. Hence we observe small, non-zero $\rho_{tp}$ even in the Mott phases. The insulating character of the phases is confirmed by the vanishing of the total stiffness. It should be pointed out that the individual stiffness values become much larger in the non-Mott phases, along with finite value of the total stiffness..

In our model, the transitions between various density wave phases are modified by the appearance of intervening supersolid phases. This is aptly demonstrated in the ground-state phase diagram in the parameter space of the pair-hopping amplitude, $t_{p}$, and the chemical potential, $\mu$, at fixed $t$ , $V=6t$ (the nn repulsion) and $V_{d}=5/3V$ (the nnn repulsion), shown in Fig. 4(a). The pair hopping and chemical potential are expressed in units of $V$. Three distinct Mott insulating lobes are present, corresponding to different filling factors, as the chemical potential $\mu$ is varied. The solid phases are destabilized with an increasing $t_{p}$, as the large pair-hopping amplitude suppresses any crystallization of the bosons into a density wave. This is manifested by the predominantly SF character of the ground state at large $t_{p}$. At sufficiently low $t_{p}$ and $\mu$, the system is in a gapless SF phase, with zero energy cost to the addition of a boson. With increasing $\mu$, there is a transition into a $\frac{1}{4}$ solid phase, which is characterized by a vanishing stiffness and finite gap. The bosons form a density wave with a pattern schematically shown in Fig. 4(b). Increasing $\mu$, the system undergoes a transition to a SS phase, which is characterized by a finite solid order and co-existing superfluid density. Due to $V_{d}>V$, the wave vector of the underlying solid order is $(\pi,0)$, reflecting a striped solid. With further increase in $\mu$, a discontinuous transition drives the ground state to a $\frac{1}{2}$ solid, where the nnn repulsion crystallizes the bosons into alternating stripes, as shown in Fig. 4(c). Finally, another SS phase with $(\pi,0)$ ordering separates the $\frac{1}{2}$ solid and the $\frac{3}{4}$ solid at large values of $\mu$. The boson ordering of the $\frac{3}{4}$ solid is shown in Fig. 4(d).

Significantly, no new phases – such as additional density wave phases – are stabilized by the introduction of the pair-hopping process. The boson ordering of the solid phases remain unchanged by pair-hopping as well. This is in contrast to the case of the hard-core boson model with ’pair rotation’, which flips two bosons residing on opposite diagonal corners of a plaquette to the other diagonal on the same plaquette. In the mentioned model, even without any diagonal interaction terms, the two-body kinetic term can induce new solid phases Sandvik et al. (2002); Sandvik and Melko (2006).

A key observation from Fig. 4 is that the different Mott lobes are modified differently by the pair hopping process from their counterparts when only $t$ is present. This is nicely illustrated by the observation that the $\langle n\rangle=\frac{1}{2}$ lobe is significantly larger than the $\frac{1}{4}$ and $\frac{3}{4}$ lobes, and persists in a larger range of $t_{p}$ and $\mu$. This has important implications in realistic models with long range interactions. While the $t$-only model may exhibit several plateaus, their extent will be heavily modified by any pair-hopping process, including the possible suppression of some of them. Another interesting feature is that all the transitions into and out of the Mott phases are discontinuous in nature. This is analogous to meta-magnetism in spin models Iaizzi et al. (2017, 2018) and we plan to investigate this further in future studies.

Magnetization plateaus in spin models manifest as boson density plateaus in the boson model. To demonstrate the dynamics of pair hopping on the density plateaus, we plot the full range of observables in Fig. 5, with parameters equivalent to taking a slice of constant $t_{p}=4$ in the phase diagram of Fig. 4(a). For the parameters chosen, we observe the existence of the aforementioned $\langle n\rangle=\frac{1}{4}$, $\frac{1}{2}$ and $\frac{3}{4}$ density plateaus in Fig. 5(a). We note that the $\langle n\rangle=\frac{1}{4}$, $\frac{1}{2}$ plateaus correspond to the $m/m_{s}=\frac{1}{4}$ and $\frac{1}{2}$ plateaus proposed by other studies Miyahara and Ueda (2000). Discontinuities in the first derivative of the density and total stiffness, $\rho_{s}$, indicate discontinuous phase transitions into and out of the three solid phases.

To study the solid ordering in the various plateaus, we plot the structure factor $S(\pi,\pi)$ and $S(\pi,0)$ as a function of $\langle n\rangle$ in Fig. 5(c). A finite $S(\pi,\pi)$ corresponds to a checkerboard boson ordering, while a finite $S(\pi,0)$ corresponds to striped boson ordering. As we have set $V_{d}>V$ in this simulation, the striped ordering out competes the checkerboard ordering and we observe a striped solid at $\frac{1}{2}$-filling, characterized by a peaked $S(\pi,0)$. The total stiffness vanishes in this phase, demonstrating the gapped nature of the striped solid, where it is energetically prohibitive to add another boson.

Compared to the $\langle n\rangle=\frac{1}{2}$ solid, the situation is markedly different for $\langle n\rangle=\frac{1}{4}$ and $\frac{3}{4}$. The $\frac{1}{4}$ solid is stabilized by bosons avoiding both nn ($V$) and nnn ($V_{d}$) repulsive interactions Dang et al. (2008), which is obvious in Fig. 4(b). The $\frac{1}{4}$ solid is then gapped as the addition of one boson incurs energy costs of either $2V-\mu$ or $4V_{d}-\mu$, depending on the neighborhood configuration of the site chosen. On the other hand, the $\frac{3}{4}$ solid manifests as a sequence of alternating fully-filled and half-filled stripes, as shown in Fig. 4(d). It is clear from the figures that the two phases are related by a particle-hole symmetry. Again, the gapped nature of the $\frac{3}{4}$ solid is obvious, as the addition of a boson incurs an energy cost of $4V+4V_{d}-\mu$. The gapped nature of both phases is also evident by the vanishing stiffness shown in Fig. 5.

To characterize the relative magnitudes of single and pair particle flow, we plot $\rho_{t}$ and $\rho_{tp}$ separately in Fig. 5(b). We note that single and pair particle flow co-exist at all values of $\mu$ for non-zero $t$ and $t_{p}$. The two currents reinforce each other in the SF phase, resulting in a total stiffness which is greater than the individual contributions from the single particle and pair currents. One could also have counter-propagation that would cause partial cancellation of the currents and a smaller $\rho_{s}$ than $\rho_{t}$ and $\rho_{tp}$. Here, we observe this effect in the $\frac{1}{2}$ striped solid, where $\rho_{t}=\rho_{tp}$ are non-zero, yet $\rho_{s}$ vanishes. The origin of this counter-propagation is due to trivial local fluctuations in the single and pair currents, such that in the solid phase, two single boson hop fluctuations that break the staggered density pattern is exactly cancelled by a pair boson hop in the opposite direction, and vice-versa. This effectively conserves a vanishing $\rho_{s}$ in the solid, even while $\rho_{t}$ and $\rho_{tp}$ are non-zero.

The stiffness plots exhibit a reflection symmetry about the $\frac{1}{2}$ solid. At small and large filling factors, $\rho_{t}$ is larger than $\rho_{tp}$, despite the fact that $t_{p}=4t$. This is explained as at low fillings, boson occupancy is sparse, making it unlikely for bosons to meet as nn pairs. At large filling factors, the lattice becomes crowded and the presence of pairs of holes, such that boson pairs can hop to fill the holes, become unlikely. This results in a larger $\rho_{t}$ than $\rho_{tp}$, despite the significantly larger pair hopping amplitude $t_{p}$. It is at intermediate filling factors where $\rho_{tp}>\rho_{t}$, in a phase we call ’pair superfluidity’. Intermediate filling factors satisfy the conditions that the lattice is neither too sparse or too crowded, thus being conducive for pair boson hopping. A finite-size scaling analysis of $\rho_{tp}$ in Fig. 6 shows that pair superfluidity is not merely a finite-size effect, and the phase extends to the thermodynamic limit. Additionally, we find that by varying the Hamiltonian parameters for a large range of values (not shown), pair superfluidity is achieved only when we tune $t_{p}\gtrapprox 4t$.

A parameter set that stabilizes a checkerboard solid at $\frac{1}{2}$ filling is shown in Fig. 7. The absence of the $\frac{1}{4}$ and $\frac{3}{4}$ plateaus in the density plot are a result of the lack of simultaneous nn and nnn repulsion, which as mentioned are necessary in the formation of these solid phases. However, a checkerboard solid at $\frac{1}{2}$ filling is still stabilized, characterized by a strongly peaked $S(\pi,\pi)$.

In this simulation, we notice that despite having $t_{p}/t=4$, as with the simulation in Fig. 5, $\rho_{tp}$ is significantly smaller than $\rho_{t}$ in the SF phases. Importantly, in the SF phase, $\rho_{tp}<\rho_{t}$ for all $\mu$ points, even at intermediate SF filling factors where it was mentioned to be the most favorable for boson pair hopping. This observation is due to the difference in boson ordering for parameters that stabilize a checkerboard and striped solid at half filling. In the SF phase of the former case, bosons will still satisfy a checkerboard ordering as far as possible to minimize nn repulsions. In a checkerboard-like configuration, bosons are largely not occupying nn sites and pair-hopping of bosons then becomes impossible, despite a large $t_{p}$. This results in $\rho_{tp}$ being significantly suppressed. In the latter case, a striped-like ordering in the SF phase implies that bosons are largely paired up, allowing pair hopping of bosons to occur more frequently. This results in more dispersive pair hopping of bosons, and subsequently a larger $\rho_{tp}$. Therefore, a dispersive pair hopping of bosons is stabilized by large $t_{p}$ and $V_{d}$ in the Hamiltonian of Eq. 1. Incidentally, we observe the same effects of currents counter-propagation in the $\frac{1}{2}$ checkerboard solid as we did in the striped solid in Fig. 5, as evident by $\rho_{t}=\rho_{tp}$. The mechanism by which counter-propagation manifests in the checkerboard solid is identical to that of the striped solid.

VI Discussion and Conclusion

Our results provide useful insight into the role of pair hopping process on the ground state phases of a system of interacting hard-core bosons. As discussed earlier, the microscopic origin of magnetization plateaus in the frustrated quantum magnet \ceSrCu_2(BO3)_2 remains incompletely understood. Neutron scattering experiments show that single magnon excitations are almost completely dispersionless. Naively, one might expect this to result in a glassy dynamics in the presence of a magnetic field. Interestingly, the same neutron scattering experiments reveal that bound states of two magnons have pronounced dispersion, although the cost of formation of such pairs is high. This provides a potential mechanism for the delocalization of field-induced triplons necessary for the long ranged ordering of the triplons at the magnetization plateaus. However, important questions remain: does the dispersion of bound pairs retain the stability of the plateaus? Do they result in new processes that are inconsistent with experimental observations in \ceSrCu_2(BO3)_2? While this is observed experimentally, such a pair hopping term has never been incorporated into any model hamiltonian describing \ceSrCu_2(BO3)_2. As such, a rigorous microscopic simulation such as ours studying the effects of such a process is valuable. Our Hamiltonian, Eq. (1), mimics the dispersion of bound pairs as a pair-hopping process. In keeping with the experimental observations, the amplitude of the pair hopping process is chosen to be much greater than the single particle hopping process. The high energy of formation is reflected in the finite near-neighbor repulsion, $V$. Our results demonstrate conclusively that highly dispersive magnon bound pairs are compatible with the formation of magnetization plateaus. However, the exact sequence of plateaus observed in \ceSrCu_2(BO3)_2 is different from the results obtained here. It is highly plausible that one needs to introduce longer-range interactions, beyond what a 4-site plaquette QMC scheme can accommodate, to fully obtain the plateaus observed in \ceSrCu_2(BO3)_2. Hence, while our results do not provide a comprehensive understanding of all plateaus in the experimental system, it provides a plausible explanation for their formation mechanism in the absence of any significant single magnon dispersion. We have also demonstrated that $t_{p}$ is important in governing which plateaus actually survive. In principle one might have a huge number of plateaus for realistic interactions with only $t$, but $t_{p}$ has different effects on different Mott phases and some of them will be destroyed by $t_{p}$ even though they survive in the presence of $t$ only.

In summary, we have investigated the role of a pair hopping process on the ground state phases of interacting hard core-bosons on a 2D square lattice. Our results may provide useful insight into the mechanism of delocalization of field-induced triplons in the frustrated magnet, \ceSrCu_2(BO3)_2, necessary for the formation of periodic patterns observed in the magnetization plateaus.

Acknowledgements.

A.J.R.H. thanks the CN Yang Scholars Programme for financial support. W.G. was supported by the NSFC under Grant No. 11775021 and No. 11734002. A.W.S. was supported by the NSF under Grant No. DMR-1710170 and by the Simons Foundation. P.S. acknowledges financial support from the Ministry of Education via Grant No. MOE2015-T1-001-056. W.G. and P.S. would also like to thank Boston University’s Condensed Matter Theory Visitors Program for support. The simulations in this work were carried out at computing facilities of the Nanyang Technological University and the National Supercomputing Center, Singapore.