Phase diagram of hard core bosons with anisotropic interactions

Lattice hard core bosons have been studied theoretically for decades, as a minimal model of a strongly interacting system displaying a superfluid phase at low temperature, as well as a quantum phase transition between superfluid and crystalline ground states. Moreover, significant insight has been gained from the investigation of this model into the possible occurrence of a supersolid phase, and on the role played by point defects such as vacancies and interstitials in its stabilization [1, 2, 3, 4, 5, 6, 7, 8, 9].
For a long time, lattice models were regarded as only of academic interest, unable to provide a realistic microscopic description of any actual physical system. However, this state of affairs has drastically changed over the past two decades, as optical lattice technology makes it now feasible, using ultracold atoms and molecules, to create artificial many-body systems accurately realizing the physics contained in these models [10, 11, 12, 13, 14]. Therefore, one may on the one hand, compare theoretical predictions and experimental detection with a degree of accuracy not attainable in experiments with naturally occurring physical system (e.g., solid helium), thereby probing even subtler aspects of the theory; in addition, the possibility of engineering interactions among cold atoms[15, 16, 17, 18, 19] not occurring in ordinary condensed matter (at least not as the dominant interactions) paves the way to the possible observation of novel, exotic phases of matter [20, 21, 22].
For example, experimental advances in the spatial confinement and cooling of large assemblies of atoms with finite electric and magnetic dipoles, make it now possible to study many-body systems whose interaction is primarily dipolar, i.e., anisotropic (attractive or repulsive depending on the relative direction of approach of two particles) and, while not strictly long-ranged, nonetheless extending to distances significantly beyond nearest neighbors, at experimentally attainable densities. These two features of the pairwise interaction make it in principle possible to stabilize different crystalline and/or superfluid phases, breaking rotational symmetry.

The physics of a system of spin zero particles possessing a dipole moment, with all dipole moments aligned along the same direction, has been extensively explored, both theoretically [23, 24] and experimentally (see, for instance, Ref. [25]). In three dimensions, a supersolid phase has been observed in computer simulations of continuous systems [26, 27]. In two dimensions, the situation is more complex; possible supersolid phases arising from speculated “microemulsion” scenarios, in systems with dipole moments aligned perpendicularly to the plane of particle motion [28], have been shown not to occur [29]. On the other hand, evidence of supersolid phases has been reported in computer simulations of dipolar bosons on the triangular lattice [30], whereas on the square lattice no supersolid phases can be stabilized [31], even if dipole moments are tilted with respect to the direction perpendicular to the plane of particle motion, a prediction confirmed for continuous systems [32]¹¹1In three dimensions, as well as lower dimensions if the tilting angle exceeds a critical value, the dipolar interaction becomes purely attractive along specific directions, making the system unstable against collapse. In these cases, the dipolar interaction must be supplemented by a short-range hard core repulsion, whose presence is customarily assumed in standard theoretical studies, to ensure thermodynamic stability..
An interesting theoretical issue that can be addressed, is the relative importance of the long range of the interaction in stabilizing specific, exotic thermodynamic phases, e.g., the supersolid. For example, the phase diagram of lattice dipolar bosons of spin zero, with dipole moments aligned perpendicularly to the plane of motion, is nearly the same as that with only nearest-neighbor interactions, both on the square [4, 31] and on the triangular lattice [9, 30]. In particular, interstitial supersolid phases exist on the triangular lattice, at particle density 1/3 (2/3), while they are absent at half filling on the square one. Remarkably, vacancy and interstitial supersolid phases can be stabilized on the square lattice at density 1/4 (3/4) if only nearest and next-nearest neighbor interactions are included [33, 34].
In this paper, we report results of a theoretical investigation of the phase diagram of a system of hard core bosons on the square lattice, interacting via an anisotropic nearest-neighbor potential. We allow the interaction to be repulsive in one direction and attractive in the other. We perform computer simulation of the system using state of the art computational technology, and map out the complete finite temperature phase diagram, which we find to be qualitatively very similar to that of a system of hard core dipolar bosons with aligned dipole moments, tilted with respect to the perpendicular to the plane. The interplay of attraction and repulsion along different directions leads to the occurrence of three phases, a superfluid (SF) one, as well as two crystalline phases at half filling, specifically a checkerboard (CB) and a striped (ST) phase. No supersolid phase is found, either at exactly half filling, consistently with fundamental theoretical arguments [35], or by doping either the CB or ST phase with vacancies or interstitials. A conventional first order quantum phase transition separates both crystalline phases from the superfluid at zero temperature. Basic features of the finite temperature phase diagram are also consistent with what observed in other lattice hard core boson models.
This paper is organized as follows: in section 2 we describe the model of interest and briefly review the computational technique utilized; in section 3 we present our results, outlining our conclusions in section 4.

The model of our interest is the well-known lattice hard core Bose Hamiltonian, expressed as follows:

	$\displaystyle\hat{H}=$	$\displaystyle-$	$\displaystyle t\sum_{{\bf r}}\biggl{(}\hat{a}^{\dagger}_{{\bf r}+\hat{\bf x}}\ \hat{a}_{{\bf r}}+\hat{a}^{\dagger}_{{\bf r}+\hat{\bf y}}\ \hat{a}_{{\bf r}}+h.c.\biggr{)}+$		(1)
		$\displaystyle+$	$\displaystyle\sum_{\bf r}\biggl{(}V_{{\bf r},\hat{\bf x}}\ \hat{n}_{\bf r}\ \hat{n}_{{\bf r}+\hat{\bf x}}+V_{{\bf r},\hat{\bf y}}\ \hat{n}_{\bf r}\ \hat{n}_{{\bf r}+\hat{\bf y}}\biggr{)}+$
		$\displaystyle-$	$\displaystyle\mu\sum_{\bf r}\hat{n}_{\bf r},$

where the sum runs over all the sites of a square lattice of $N=L\times L$ sites, with periodic boundary conditions. Here, $\hat{\bf x}$ and $\hat{\bf y}$ are the unit vectors in the two crystallographic directions, $\hat{a}^{\dagger}_{\bf r}$, $\hat{a}_{\bf r}$ are the standard Bose creation and annihilation operators, $n_{\bf r}\equiv\hat{a}^{\dagger}_{\bf r}\ \hat{a}_{\bf r}$ is the occupation number for site ${\bf r}$, $t$ is the particle-hopping matrix element and $\mu$ is the chemical potential. The (nearest-neighbor) interaction potential is defined as follows:

$$V_{{\bf r},\hat{\bf x}}=V,\ \ \ V_{{\bf r},\hat{\bf y}}=\lambda\ V,$$

(2)

$\lambda$ being a real number. The Hamiltonian (1) is defined in the subspace of many-particle configurations in which no lattice site can be occupied by more than one particle, which expresses the hard core condition. Thus, the total number of particles $N_{P}\equiv\sum_{\bf r}\hat{a}^{\dagger}_{\bf r}\hat{a}_{\bf r}$ can take on any integer value from 0 to $N$. The particle density (or, filling) is defined as $\rho\equiv(N_{P}/N)$, i.e., $0\leq\rho\leq 1$.

The phase diagram of this model has been extensively investigated in the isotropic $\lambda=1$ case; it is known that, for $-2t<V<2t$, the system displays a single homogeneous fluid phase, undergoing a SF transition at low temperature. On the other hand, if $V>2t$ a first order quantum phase transition occurs at temperature $T=0$ between a SF and a CB crystal at half filling (i.e., $\rho=1/2$) [4]. For $V<-2t$, on the other hand, only two phases occur with $\rho=0$ or 1, phases of intermediate filling (some of them SF) only stabilized by external agents, such as disorder [36].
The $\lambda\neq 1$ case, on the other hand, is relatively unexplored, mainly because for a long time a possible experimental realization did not seem possible or straightforward. However, an interaction with that kind of anisotropy can be obtained, in dipolar systems, simply by tilting the direction of alignment of dipole moments with respect to the direction perpendicular to the plane of the lattice (see, for instance, Ref. [31] for details). While the dipolar interaction extends significantly beyond nearest neighboring lattice sites, in this work we are interested in assessing which features of the phase diagram arise simply as a result of anisotropy in the short range part of the interaction. Of particular interest is the case $\lambda<0$, i.e., the interaction is repulsive (attractive) in the $x$ ($y$) direction.

In this work, we systematically investigated the finite temperature phase diagram of Eq. 1, as a function of the three parameters $T$, $\mu$ and $\lambda$, by means of first principle numerical (Monte Carlo) simulations. We used the worm algorithm in the lattice path-integral representation [37], specifically its implementation described in Ref. [38]. In order to characterize the various phases of the system we computed the superfluid density $\rho_{S}(T)$ using the well-known winding number estimator [39], as well as the static structure factor

$$S({\bf Q})=\frac{1}{N^{2}}\ \biggl{\langle}\biggl{|}\sum_{\bf r}\hat{n}_{\bf r}\ e^{i{\bf Q}\cdot{\bf r}}\biggr{|}^{2}\biggr{\rangle}$$

(3)

where $\langle...\rangle$ stands for thermal average. Presence of crystalline long-range order is signaled by a finite value of $S({\bf Q})$ for some specific wave vector(s). In particular, ${\bf Q}=(\pi,\pi)$ is the wave vector associated to checkerboard order, while ${\bf Q}=(\pi,0)$ (and not $(0,\pi)$ because of the way we have defined the interaction, Eq. 2) signals stripe order, in both cases at half filling. Numerical simulations have been carried out on square lattices of different sizes ($L=32$ being the largest) in order to gauge the importance of finite size effects.

We begin by discussing the physical behavior of the system in its ground state, i.e., at temperature $T=0$. Because we make use of a finite temperature simulation method, this goal is formally accomplished by extrapolating to the $T=0$ limit results obtained at sufficiently low temperatures. In practice, we observe that numerical estimates obtained at temperature $T\sim T_{L}\equiv(t/L)$ are indistinguishable from the extrapolated ones, within the statistical uncertainties of the calculation, i.e., physical estimates obtained at temperature $T_{L}$ can be regarded as essentially ground-state estimates.
As an example of the procedure followed to characterize the phases of the system, Fig. 1 illustrates ground state results for the case $\lambda=-1$ and $V=3t$. In this case, the system displays two distinct phases, a superfluid one (possessing no crystalline order) away from half filling, and a non-superfluid crystalline (striped) one at half filling, with a first order phase transition between the two.
Specifically, the top panel shows the particle density $\rho$ computed as a function of the chemical potential $\mu$ (in units of $t$). The density jump as $\rho$ approaches half filling signals a first-order phase transition, with coexistence of a superfluid phase (middle panel shows the superfluid density $\rho_{S}$) and one that features striped crystalline order (bottom panel shows the static structure factor $S(\pi,0))$. As one can see, no supersolid phase exists, i.e., in which both $\rho_{S}$ and $S(\pi,0)$ are finite. The superfluid density is everywhere finite, except at half filling, whereas the static structure factor is finite only at half filling. This can be contrasted with the physics described in Fig. 2 for $V/t=3$ and $\lambda=0$. In this case, all curves are everywhere continuous, i.e., no phase transition occurs.

Fig. 4 shows two representative cuts of the finite-temperature phase diagram of the system in the ($T,\mu,V$) space. We consider for definiteness the case $\lambda=-1$, for which the ground state at half filling is a stripe solid for $V\gtrsim 2.5\ t$. The left panel has fixed chemical potential ($\mu=t$), while the right panel has fixed interaction strength ($V=3\ t$).

At constant temperature, the striped crystal is separated from the (super)fluid phase by a first-order (quantum) phase transition, shown by solid lines in Fig. 4. On the other hand, at constant $V$ or $\mu$ superfluid order is destroyed by thermal fluctuations in favor of a normal fluid via a Berezinskii-Kosterlitz-Thouless transition, shown by the dotted lines, while the striped crystal melts into a normal fluid via a second-order transition in the Ising universality class (dashed curve in Fig. 4). Altogether, this is a fairly conventional phase diagram, qualitatively identical to that in the $\lambda>0$ sector, with no supersolid phase.

In this work we investigate the phase diagram of a system of hard core bosons on the square lattice, with a short-ranged (nearest-neighbor) interaction that can vary independently along the two crystallographic directions, to include the case of attraction along one direction and repulsion along the other. This kind of anisotropic interaction can be realized experimentally with a system of dipolar atoms or molecules, confined to a 2D optical lattice, if dipole moments are aligned at a varying (tilting) angle with respect to the direction perpendicular to the plane of particle motion.
The main goal of this study are a) to gain understanding of the phases that can arise as a result of the anisotropy of the interaction, b) to gauge the importance of its extended range, by comparing our results to those obtained in similar studies [31] in which the full dipolar interaction is considered. With respect to the case in which the interaction is isotropic, the only difference is the appearance of a different crystalline phase, displaying stripe order. Otherwise, the phase diagram is qualitatively identical to that for isotropic interaction, with a single crystalline phase at half filling, and no supersolid phase. The phase diagram is also qualitatively identical to that observed by allowing the interaction to take on the full dipolar form, i.e., extending beyond nearest neighbors, with varying degree of anisotropy corresponding to different tilting angles [31]. This is consistent with what observed for purely repulsive, isotropic interactions, where the presence of a $1/r^{3}$ dipolar tail does not fundamentally alter the phase diagram, e.g., on the triangular lattice [30]. Altogether, this study reinforces the notion that lattice geometry (i.e., triangular versus square), more than the detailed form of the interaction, is at the root of the occurrence of supersolid phases.
It has been suggested in a recent study [40], that in the presence of anisotropic interactions of dipolar, form supersolid phases may be stabilized, away from half filling, by allowing the interactions to extend to next-nearest neighboring distances. We have not explored this aspect here, but the contention seems plausible and not particularly surprising, given that this was shown to be the case for isotropic repulsive interactions [34] (i.e., the $\lambda>0$ sector of Eq. 2), and the fact that the physics of the system in the positive and negative $\lambda$ sectors is essentially the same, as shown in this study.

This work was supported by the Natural Sciences and Engineering Research Council of Canada. Computing support of ComputeCanada is acknowledged.