Ising Antiferromagnet in the 2D Hubbard Model with Mismatched Fermi Surfaces

TL;DR

This paper investigates the phase diagram of a 2D Hubbard model with spin-dependent anisotropic hopping, revealing the emergence of Ising antiferromagnetic order at infinitesimal interactions and providing critical temperature predictions using advanced quantum Monte Carlo methods.

Contribution

It provides the first unbiased predictions of critical temperatures for Ising antiferromagnetism in this model using multiple state-of-the-art quantum Monte Carlo techniques.

Findings

Ising antiferromagnetic order appears at infinitesimal repulsive interactions.

Critical temperatures are predicted for various interaction strengths.

Results are relevant for ultracold atom experiments with spin-dependent lattices.

Abstract

We study the phase diagram of the two-dimensional repulsive Hubbard model with spin-dependent anisotropic hopping at half-filling. The system develops Ising antiferromagnetic long-range order already at infinitesimal repulsive interaction strength in the ground state. Outside the perturbative regime, unbiased predictions for the critical temperatures of the Ising antiferromagnet are made for representative interaction values by a variety of state-of-the-art quantum Monte Carlo methods, including the diagrammatic Monte Carlo, continuous-time determinantal Monte Carlo and path-integral Monte Carlo methods. Our findings are relevant to ultracold atom experiments in the p-orbital or with spin-dependent optical lattices.

Tables

Table 1. Table 1: A comparison of the QMC methods used. For the Worm and LCT-QMC methods the sign-positive regimes are mentioned in the table. They scale linearly and cubically in the system volume, respectively, and both linearly with the inverse temperature. DiagMC simulations work directly in the thermodynamic limit. In practice, open boundary conditions are used in the Worm simulations.

Method	Anisotropy	Filling	Interaction
Worm Prokof’ev et al. (1998); Pollet et al. (2007)	$\alpha =0$	arbitrary	arbitrary
DiagMC Van Houcke et al. (2010); Gukelberger (2015)	arbitrary	arbitrary	$U\lesssim 4t$
LCT-QMC Iazzi and Troyer (2015); Wang et al. (2015a)	arbitrary	half filling	arbitrary

Equations

$$\hat{H}=$$

$$H_{\alpha =0}$$

$$\chi (Q)$$

$$H_{\uparrow }^j=\frac{1}{2\pi }\int dx{u}_{\uparrow }^j{K}_{\uparrow }^j(\nabla {\theta }_{\uparrow }^j(x){)}^2+\frac{u_{\uparrow }^j}{K_{\uparrow }^j}(\nabla {\varphi }_{\uparrow }^j(x){)}^2,$$

$${\rho }_{\uparrow }^j(x)={\rho }_0-\frac{1}{\pi }\nabla {\varphi }^j(x)+{\rho }_0\sum _{p\ne 0}{e}^{i2p(\pi {\rho }_{0}x-{\varphi }_{\uparrow }^j(x))},$$

$$M_{\mathrm{st}}=\sum _{\mathbf{r}=(x,y)}(-1{)}^{x+y}({\hat{c}}_{\mathbf{r},\uparrow }^{†}{\hat{c}}_{\mathbf{r},\uparrow }-{\hat{c}}_{\mathbf{r},\downarrow }^{†}{\hat{c}}_{\mathbf{r},\downarrow }).$$

$${\chi }_{\uparrow \uparrow }^0(\mathbf{Q})$$

$$\Gamma (Q{)}_{\sigma ,k;{\sigma }^{\prime },k^{\prime }}=U{\delta }_{\sigma ,-{\sigma }^{\prime }}+\mathcal{O}(U^2).$$

$$\lambda$$

$$T_c$$

Full text

I Introduction

The Hubbard model plays an important role in condensed matter research combining a plethora of physical phenomena such as the Mott insulator transition, magnetism and (un)conventional superconductivity. Despite this richness, exact analytical solutions have only been found in a few special cases, including the one-dimensional model Essler et al. (2005), the atomic, and the non-interacting limit. Furthermore, the Hubbard model is numerically tractable in infinite dimensions by dynamical mean-field theory Georges et al. (1996), on bipartite lattices at half-filling by determinantal quantum Monte Carlo methods, and in the infinite-$U$ limit on ladder geometries with the density matrix renormalization group Liu et al. (2012). Going away from either of these special limits poses tremendous challenges to our theoretical understanding. Reference LeBlanc et al. (2015) summarizes the presently known results of the Hubbard model from a wide range of numerical algorithms.

In this paper, we add another parameter regime which can be exactly solved numerically: we study the phase diagram of the repulsive Hubbard model with spin-dependent anisotropic hopping by three different kinds of unbiased quantum Monte Carlo (QMC) algorithms. As we will see, some of these methods work only in certain parameter regimes but can be more efficient when they are applicable. The Hamiltonian reads

[TABLE]

Specifically, we consider the spin-dependent anisotropic hopping amplitude $t_{x\uparrow}=t_{y\downarrow}=t$ and $t_{x\downarrow}=t_{y\uparrow}=\alpha t$, where $\alpha\in[0,1]$ is a tuning parameter. The hopping is stronger along x(y) direction for spin up(down) fermions, shown in Fig. 1(a). It leads to a spin-dependent nematic distortion of the Fermi surface in the reciprocal space shown in Fig. 1(b). Physically, cold atomic systems Bloch et al. (2008); Esslinger (2010) may be well suited to study this system. The Hubbard model has been realized years ago Jordens et al. (2008); Schneider et al. (2008) and with fermionic microscopes antiferromagnetic correlations have been measured Boll et al. (2016); Parsons et al. (2016); Cheuk et al. (2016), which can now extend over the entire system size and realize a Heisenberg antiferromagnet Mazurenko et al. (2016). The hopping anisotropy can either be realized with spin-dependent optical lattices, or due to the anisotropic shape of the Wannier function on the $p$-band of an optical lattice. The last term of (1) denotes an onsite repulsive interaction with $U>0$. We focus on magnetic order of the model (1) in the half-filled case on a square lattice.

Diagonalizing the single-particle part of the Hamiltonian, one has dispersions $\epsilon_{k}^{\uparrow}=-2t\cos(k_{x})-2\alpha t\cos(k_{y})$ and $\epsilon_{k}^{\downarrow}=-2\alpha t\cos(k_{x})-2t\cos(k_{y})$. The noninteracting bandwidth is thus $W=4(1+\alpha)t$. Spin nematicity explicitly breaks the spin SU(2) symmetry and removes the divergence in the density of states at the Fermi energy. However the perfect Fermi surface nesting with wave vector $\boldsymbol{Q}=(\pi,\pi)$ is still preserved. As a consequence, the longitudinal spin susceptibility $\chi^{\sigma\sigma}(\boldsymbol{Q})=\frac{1}{N}\sum_{k}\frac{f(\epsilon_{k}^{\sigma})-f(\epsilon^{\sigma}_{k+Q})}{\epsilon_{k}^{\sigma}-\epsilon_{k+Q}^{\sigma}}$ still diverges at zero temperature while the transverse spin susceptibility $\chi^{+-}(\boldsymbol{Q})=\frac{1}{N}\sum_{k}\frac{f(\epsilon_{k}^{\uparrow})-f(\epsilon^{\downarrow}_{k+Q})}{\epsilon_{k}^{\uparrow}-\epsilon_{k+Q}^{\downarrow}}$ saturates to a constant value (Here, $f(\cdot)$ is the Fermi-Dirac distribution). Therefore, a weak-coupling analysis predicts Ising antiferromagnetic (AF) order already at infinitesimally strong repulsive interaction.

The strong coupling limit of the model (1) was studied in the context of $p$-orbital Mott insulators Wu (2008); Zhao and Liu (2008). It reduces to a spin-$1/2$ XXZ model with Ising anisotropy, which favors an antiferromagnetic Ising ground state. At intermediate interaction strength, the system exhibits a crossover from a weak-coupling spin-density-wave state to the strong-coupling AF Ising state, similar to the case of the 3D half-filled Hubbard model. However, since the Ising state breaks only a discrete $\mathbb{Z}_{2}$ symmetry, it has a finite critical temperature, even in two dimensions.

The above considerations continue to hold in the fully anisotropic case $\alpha=0$. In this limit, the kinetic part is purely one-dimensional – i.e. the motion of a spin is limited to a row or a column of the 2D lattice – whereas the density-density interactions on each site connect the two spin species and make the system effectively two dimensional. Therefore, as we will show with unbiased worldline QMC simulations, the system still possesses a finite critical temperature.

References Feiguin and Fisher (2009, 2011); Gukelberger et al. (2014) studied the model (1) in the doped attractive case in search of an elusive Bose liquid and exotic superfluid states. It turns out that close to half filling the most stable phase is an incommensurate density wave state, which is related to the AF Ising state of the repulsive model discussed above upon a particle-hole transformation. It was also remarked in Ref. Gukelberger et al. (2014) that in the fully anisotropic limit the particle number of each column and row is separately conserved. This allows one to perform unbiased worldline QMC calculations by mapping the fermions to quantum spins, similar to what was done in Xu et al. (2015), and which served as a benchmarking tool for the diagrammatic Monte Carlo calculations.

The model (1) is different from the one studied in Ref. Berg et al. (2012), which studied onset of antiferromagnetism in a similar Fermi surface geometry. However, Ref. Berg et al. (2012) doubles the fermion species to avoid the fermion sign problem. In addition, the current study focuses on the half-filled case where the AF Ising state is strongly enhanced due to the commensurate filling.

The organization of the paper is as follows. In Sec. II we summarize the Monte Carlo methods used in this paper and comment on their advantages and disadvantages. In Sec. III we report on results obtained by various QMC calculations, where Sec. III.1 contains results on the fully anisotropic case of the model (1) and Sec. III.2 results for general anisotropic cases. Section IV summarizes our main findings and discusses their implications for future experimental and theoretical studies.

II Methods

In this section we summarize the three different quantum Monte Carlo techniques used to study the model Eq. (1): Path Integral Monte Carlo simulations with worm-type updates (Worm), diagrammatic Monte Carlo simulations (DiagMC), and continuous-time determinantal Monte Carlo simulations (LCT-QMC). Table 1 summarizes their main features and allows one to quickly read off the method of choice. In their domain of applicability all three methods yield unbiased results on the physical observables. Whenever there is an overlap in their application range we have checked that they give consistent results. In the subsections below we explain in more detail the specifics of all three methods for the anisotropic Hubbard model.

II.1 Path-integral Monte Carlo (Worm)

In the fully anisotropic limit, which is where the Worm algorithm can be applied, the model Eq. (1) reduces to

[TABLE]

The hopping is one-dimensional, implying that for each row (column) the number of up (down) particles is conserved. By translational invariance we expect that all or none of these symmetries are simultaneously broken. As a consequence of the 1D character, individual rows and columns can be mapped onto hard-core bosons at any density through the celebrated Jordan-Wigner transformation Cazalilla et al. (2011); Giamarchi (2003), which in turn allows us to use path-integral Monte Carlo simulations with worm-type updates Prokof’ev et al. (1998), here in the implementation of Ref. Pollet et al. (2007). Spin densities and density-density correlations functions, which we measure in order to identify the phase transition, are not affected by the Jordan-Wigner transformation and identical for the original fermions and the simulated hard-core bosons. For ease of the Jordan-Wigner transformation, we use open boundary conditions. This comes at the price of greater finite size effects through the influence of the boundary terms, which is however minor in light of the mapping to a positive expansion for all filling factors and the linear scaling of the Worm algorithm with system size and inverse temperature. For $\alpha\neq 0$ the Worm algorithm has a sign problem leading to an exponential scaling in the system volume and inverse temperature.

II.2 Diagrammatic Monte Carlo (DiagMC)

The diagrammatic Monte Carlo (DiagMC) method evaluates Feynman diagrammatic expansions by means of a stochastic process that samples sums over diagram topologies and internal variables on equal grounds Prokof’ev and Svistunov (2007, 2008). Our implementation for the Hubbard model Gukelberger et al. (2014); Gukelberger (2015), which is based on diagrams with bare propagators $G_{0}$ and interactions $U$, is not directly applicable within a magnetically ordered phase. Therefore, we detect a continuous phase transition to AF order by monitoring the divergence of the magnetic susceptibility on approaching the critical temperature. This aspect is different from the other two Monte Carlo methods, which are not formulated in the thermodynamic limit. To this end we sample the self-energy $\Sigma_{\sigma}(k)$ and the irreducible scattering vertex in the particle-hole channel $\Gamma^{ph}_{\sigma\sigma^{\prime}}(Q,k,k^{\prime})$ for fixed total four-momentum $Q=(\boldsymbol{Q},i\Omega_{m}=0)$ with $\boldsymbol{Q}=(\pi,\pi)$ the AF ordering vector. According to the Bethe-Salpeter equation

[TABLE]

the susceptibility $\chi(Q)$ diverges when the largest eigenvalue of the kernel $-\chi^{0}(Q)\Gamma(Q)$ reaches unity. The above should be read as a matrix equation for the generalized susceptibility $\chi^{\sigma\sigma^{\prime}}(Q;k,k^{\prime})$ in spin and four-momentum space. Furthermore, the particle-hole bubble $\chi^{0}_{\sigma\sigma^{\prime}}(Q;k,k^{\prime})=G_{\sigma}(k+Q/2)G_{\sigma}(k-Q/2)\delta_{\sigma,\sigma^{\prime}}\delta(k-k^{\prime})$ is the diagonal product of two one-particle propagators and the one-particle propagators in turn are calculated from the self-energy via Dyson’s equation.

With DiagMC the system is directly simulated in the thermodynamic limit, but the diagrammatic series for the irreducible quantities $\Sigma$ and $\Gamma$ must be restricted to orders $n\leq N_{*}$ because the sign of a fermionic series vanishes factorially with diagram order $n$. All DiagMC results must therefore be extrapolated in the cutoff parameter $N_{*}\to\infty$. The uncertainty in this extrapolation is typically the dominant contribution to the error bars and the extrapolation may be impossible when the series does not converge quickly enough. This happens frequently if the interaction is too strong, e.g. $U\gtrsim W$. For models like the half-filled Hubbard model where determinantal QMC methods do not suffer from the sign problem, the sign-problem-free method will generally yield smaller error bars than DiagMC under comparable computational efforts. The main advantage of DiagMC is that it can equally well be applied away from half filling, where simulations with other QMC methods are often unfeasible due to a severe sign problem. Additionally, the comparison of finite-size extrapolations (e.g. from path-integral or determinantal QMC) with finite-order extrapolations from DiagMC yields a very nontrivial crosscheck that all systematic errors in the different methods are under control.

II.3 Continuous-time determinantal Monte Carlo (LCT-QMC)

We employ the continuous-time quantum Monte Carlo method scaling linearly in $\beta$ (LCT-QMC) Iazzi and Troyer (2015); Wang et al. (2015a) to study the model (1) at general anisotropies on finite lattices. The LCT-QMC methods perform continuous-time interaction expansion of the partition function and evaluate each expansion as a matrix determinant. Thanks to recent progress on the fermion sign problem Huffman and Chandrasekharan (2014); Li et al. (2015); Wang et al. (2015b); Wei et al. (2016) these matrix determinants can be shown to be nonnegative. There is no sign problem in the simulation despite the mismatched Fermi surfaces: The crucial conditions are half filling and the presence of bipartite lattices. The implementation of the LCT-QMC simulation is similar to the recent study of the mass-imbalanced Hubbard model Liu and Wang (2015). As the signature of the phase transition we measure the staggered magnetization square according to the Wick’s theorem in the LCT-QMC simulations.

Compared to the path-integral Monte Carlo method of Sec. II.1, the drawback of the LCT-QMC algorithm is that it scales cubically with the system size. We are therefore limited to system sizes $L\leq 24$ for the LCT-QMC results. The advantage, however, is that one is able to study also finite anisotropy ratios and systems with periodic boundary conditions can be simulated without further constraints because the method does not rely on the Jordan-Wigner mapping.

III Results

In this section we first present our results for the fully anisotropic case, followed by the results for the more general case. The unit of energy is set by the hopping $t=1$ unless explicitly noted otherwise.

III.1 The fully anisotropic model

Below we use bosonization arguments to get an intuitive and analytical understanding of the phase diagram at zero temperature, followed by quantum Monte Carlo simulations addressing the phase transition at finite temperature. We will see that the ground state is always gapped and ordered in spin space, whereas at finite temperature a $\mathbb{Z}_{2}$ transition between a normal liquid and an antiferromagnet is found. Unless otherwise specified, we limit ourselves to the half-filled case.

III.1.1 Bosonization considerations of the ground state

Thanks to the one-dimensional nature of the hopping, each row and column can be bosonized separately. Following the notation and the formulas of App. D in the standard book (Ref. Giamarchi (2003)) we write the harmonic action for row $j$ as

[TABLE]

where $u_{\uparrow}$ is a velocity and $K_{\uparrow}$ the dimensionless Luttinger parameter. The fields $\nabla\phi$ and $\nabla\theta$ are proportional to the sum and the difference of right and left movers, respectively. For a column $\bar{j}$ a similar expression can be written down with the replacements ${\uparrow}\leftrightarrow{\downarrow}$, $x\leftrightarrow y$ and $j\leftrightarrow\bar{j}$. We still need to investigate the Hubbard term, which couples the spin densities on intersecting rows and columns, and take care of the filling factor. The density in bosonized form is

[TABLE]

with $\rho_{0}=1/2$ at half filling.

Introducing the charge $\phi^{j\bar{k}}_{\rho}=(\phi^{j}_{\uparrow}+\phi^{\bar{k}}_{\downarrow})/\sqrt{2}$ and the spin $\phi^{j\bar{k}}_{\sigma}=(\phi^{j}_{\uparrow}-\phi^{\bar{k}}_{\downarrow})/\sqrt{2}$ fields we get a non-oscillating term $\cos(\sqrt{8}\phi^{j\bar{k}}_{\rho})$ resulting from the Hubbard interaction, as well as a term $\cos(\sqrt{8}\phi^{j\bar{k}}_{\sigma})$. If we assume that translational invariance is not broken, then the fields for all $j$ and $\bar{j}$ are the same, and the cosines become relevant in both sectors; i.e., similar to the 1D Hubbard model with spin and repulsive interactions at half filling the charge sector is always massive at zero temperature. Its gap can be exponentially small $\ln\Delta\sim-1/\sqrt{U}$ in the weak-coupling regime (cf. Eq. (10) below). However, in contrast to the 1D Hubbard model, the spin sector cannot remain a spin liquid because of the 2D nature of the lattice (which we see in the bosonization via the presence of the second cosine term). The system therefore orders into an Ising antiferromagnet in order to lower its energy. Away from half filling, similar arguments can be applied leading to incommensurate spin density waves, in line with the weak-coupling and DiagMC results of Ref. Gukelberger et al. (2014) for the attractive case.

III.1.2 Monte Carlo results for the Ising transition at finite temperature

At finite temperature one expects a phase transition between a normal liquid and an Ising antiferromagnet with critical exponents belonging to the 2D classical Ising universality class. In order to test this, we performed large scale Monte Carlo simulations using the Worm algorithm and computed the expectation value square $\langle M_{\rm st}^{2}\rangle$ and fourth power $\langle M_{\rm st}^{4}\rangle$ of the staggered magnetization $M_{\rm st}$,

[TABLE]

Finite size scaling theory predicts, in leading order, that the curves $\langle M_{\rm st}^{2}\rangle L^{2\beta/\nu}$ intersect in a single point. Here, $L$ is the linear system size of the system, $\beta$ is the critical exponent for the order parameter which is $\beta=1/8$ for the 2D classical Ising model, and $\nu=1$ is the critical exponent for the correlation length. This is shown in Fig. 2, where we see that the system sizes $L=8$ and $L=16$ are too small to be taken into account in the finite size analysis. For system sizes $L=24$ and larger we get curves that intersect, within error bars, in almost a single point when the staggered magnetization squared is multiplied with the correct power of the system size, $L^{1/4}$, in agreement with the critical exponents of the classical 2D Ising class.

Although the linear scaling in the system volume of the Worm algorithm suggests it should be the method of choice in the absence of a sign problem, the Worm algorithm is nevertheless not well equipped to study this Ising transition because the worms are confined to single rows and single columns. The spin-resolved single-particle density matrix is hence one-dimensional and decays exponentially in the gapped phase: The algorithm is in the spin sector not better than a single spin-flip algorithm for a classical 2D Ising model. We have checked for $\beta=8$ that the integrated autocorrelation time increases linearly with $L$ with a very large prefactor. Very close to the transition point, additional critical slowing down takes place with a dynamical exponent $z\approx 2$, just as in the single spin-flip algorithm for a classical 2D Ising model. To give an idea, for $L=32,\beta=9.6$, we find a value around 100 with a binning analysis, where each measured value taken into account in the binning analysis is already an average of 1000 Monte Carlo measurements. Measurements were taken after 1000 Monte Carlo worm updates to compensate for the size of the system. The total calculation lasted several CPU-months per data point and resulted in more than half a million measurements, but that is barely enough. An immediate consequence is that the fluctuations on the Binder cumulant are an order of magnitude worse than the ones in Fig. 2, and are therefore less precise to locate the phase transition. We have also successfully repeated this analysis for $U=4$ (not shown) with $\beta_{c}=6.2(5)$. We leave for future work whether a new algorithm can be devised which combines a spin-cluster algorithm with the Worm algorithm in order to overcome this critical slowing down. Despite the present algorithm’s inefficiencies we have obtained results for larger systems than accessible with any other method. The finite size scaling is further validated by DiagMC results for the thermodynamic limit in Sec. III.2.2.

We also tried a similar analysis for the ground state assuming the universality class of the 3D classical Ising spin model (not shown). For system sizes up to $L=128$ we failed to find a single crossing point: curves for the staggered magnetization squared multiplied with $L^{2\beta/\nu}$ have all the same shape with the steep part shifting parallel to lower values of $U$ with increasing $L$. This is consistent with the ground state being ordered for any $U$, in line with the bosonization arguments. Since the charge gap opens exponentially slowly for low values of $U$ there is of course no chance of observing the ground state in a brute-force numerical approach in the small $U$ limit.

III.2 General anisotropic case

At general anisotropy $\alpha>0$ both fermion species can hop in the 2D plane. Therefore an effective bosonic description can no longer hold. In the following we show that the system has a weak-coupling instability to antiferromagnetic order for all values of $\alpha$. Then we obtain unbiased results for the transition temperature at intermediate interaction using two fermionic QMC methods (LCT-QMC and DiagMC) and cross-check the results. In general one anticipates that in the highly anisotropic case ($\alpha\ll 1$) the critical temperature approaches the one determined by the bosonic Worm calculation in Sec III.1, whereas the critical temperature drops to zero when $\alpha$ approaches unity, restoring the full $SU(2)$ rotational symmetry.

III.2.1 Weak coupling analysis

Particle-hole symmetry of the half-filled model (1) ensures that the Fermi surfaces for the two spin species individually are nested with respect to the AF wave vector $\boldsymbol{Q}=(\pi,\pi)$ independent of the anisotropy $\alpha$. Therefore the longitudinal spin susceptibilities always have a logarithmic divergence

[TABLE]

for $T\to 0$. Here $T_{F}$ denotes the Fermi temperature. In contrast, nesting between $\uparrow$- and $\downarrow$-Fermi surfaces is destroyed by an anisotropy $\alpha\neq 1$, so that the transverse spin susceptibilities $\chi^{0}_{+-}$, $\chi^{0}_{-+}$ saturate to finite values at low temperature. Since there are no first-order pairing instabilities in the particle-particle channel for repulsive interactions, longitudinal (Ising) antiferromagnetism is the only instability to leading order in $U$.

In the first-order approximation to the Bethe-Salpeter kernel in the longitudinal particle-hole channel the irreducible vertex is replaced by the bare interaction

[TABLE]

This yields an eigenvalue

[TABLE]

which grows logarithmically according to (7) and will hence reach unity for arbitrarily small $U$ at a critical temperature

[TABLE]

which has the typical form of a BCS-type weak-coupling instability ($c$ is the constant prefactor of the logarithmic divergence in (7)).

In summary, a weak-coupling analysis predicts a general low-temperature instability of the Fermi liquid towards Ising-type antiferromagnetic order (cf. the previous section) — except at the isotropic point $\alpha=1$ where longitudinal and transverse channels become degenerate and magnetic order at finite temperature is ruled out by the continuous spin rotation symmetry (as long as the system remains purely two-dimensional). For weak coupling we expect the $T_{c}$ suppression to be confined to a very small region around the isotropic point because at exponentially low temperatures the physics is extremely sensitive to small Fermi surface mismatches. Away from half filling and at $\alpha>0$ the perfect nesting and hence the weak-coupling instability in the particle-hole channel is lifted. Then only second-order instabilities in the particle-particle pairing channel remain, leading to $p$-wave superfluidity in direct correspondence to the attractive case Gukelberger et al. (2014).

III.2.2 DiagMC results

In order to go beyond the weak-coupling analysis we turn to DiagMC simulations, which can address arbitrary anisotropy. As shown in Figs. 3 and 4 we track the leading Bethe-Salpeter eigenvalue for antiferromagnetic order. As the temperature is lowered, the eigenvalue grows and eventually crosses unity, causing a divergence of the AF susceptibility. While a cutoff order $N_{*}=1$ corresponds to a mean-field treatment and strongly overestimates the transition temperature, the eigenvalues for higher cutoffs converge reasonably quickly with a decaying even-odd oscillation: The eigenvalue for each order lies between the values from the next two smaller orders. We take the average of the three largest orders $N_{*}=4,5,6$ as extrapolation to infinite order and give error bars that cover these three finite-order results. For the fully anisotropic model (Fig. 3) we obtain the transition point $\beta_{c}=8.2\pm 1.7$. This is consistent with the Worm result presented above. As expected, the DiagMC error bar is markedly larger than the one obtained with sign-problem-free bosonic QMC. The results at general anisotropy $\alpha=0.75$ (Fig. 4) are very similar to the fully anisotropic case even though the kinetic terms are changed from one-dimensional to two-dimensional. In this case we obtain a slightly larger transition temperature $\beta_{c}=7.5\pm 1.4$.

III.2.3 LCT-QMC results

Using the LCT-QMC method, we obtained the critical temperature at $U/t=3,\alpha=0.75$ by scaling the staggered magnetization according to the 2D Ising critical exponent as is shown in Fig. 5. The estimate is again in agreement with the critical temperature obtained by the DiagMC calculations. Figure 6 summarizes the critical temperature computed at different anisotropic ratios. The critical temperature measured in the unit of the bandwidth $W$ remains high from extreme ($\alpha\sim 0$) to intermediate ($\alpha\sim 0.75$) anisotropy. When measured in the unit of the hopping amplitude $t$ the transition temperature even rises with decreasing anisotropy. However, since $T_{c}$ should drop to zero in the isotropic case, it suggests a quite abrupt change of the critical temperature in the neighborhood of the isotropic point $\alpha\sim 1$. This behavior is reminiscent of the XXZ model Jr et al. (1985), which applies in the strong coupling limit, and also appears natural in the weak-coupling limit (cf. Sec. III.2.1 above).

IV Summary and experimental realization

Breaking of the discrete $\mathbb{Z}_{2}$ spin inversion symmetry results in an antiferromagnetic Ising state with finite critical temperature on a two-dimensional lattice. We presented quantitative predictions for the onset of antiferromagnetic Ising order in a Hubbard model with mismatched Fermi surfaces. Since the model (1) can be implemented using spin-dependent optical lattices or higher orbitals, testing these predictions would be a step stone to further approach the more exotic quantum phases at different filling and interaction Gukelberger et al. (2014) and on different lattice geometries Wu (2008); Zhao and Liu (2008).

On the methodological side, we have shown how specific limits of the Hubbard model with mismatched Fermi surfaces on two-dimensional lattices can be brought under full numerical control by using three different quantum Monte Carlo methods. These limits are (i) fully anisotropic spin-dependent hopping at any density with the Worm algorithm, (ii) anisotropic spin-dependent hopping at half filling with LCT-QMC, and (iii) arbitrary anisotropy and density at sufficiently weak interactions using DiagMC. The results of the methods are consistent with each other within their domain of applicability and furthermore supported by analytical weak-coupling arguments.

The main physical result is that we find a discrete $\mathbb{Z}_{2}$ symmetry breaking at a finite critical temperature towards an Ising antiferromagnet in the half-filled model. In view of our numerical results at intermediate interactions as well as the situation in the weak- and strong-coupling limits, a small anisotropy seems to be generically sufficient to create a large critical temperature. One should notice that breaking the spin $SU(2)$ symmetry is crucial to obtain an Ising antiferromagnet with finite critical temperature in 2D. The strictly symmetrical distortion of the spin up and down Fermi surfaces shown in Fig. 1 is however not crucial. Away from half filling the weak-coupling limit predicts incommensurate spin density waves at extreme anisotropy $\alpha=0$ and $p$-wave superfluidity at general anisotropy $0<\alpha<1$. At finite interactions we expect the $\mathcal{O}(U)$ spin density wave instability to dominate over the $\mathcal{O}(U^{2})$ superfluid instabilities in a finite region of the phase diagram at strong anisotropy and around half filling, similar to the situation in the attractive case 111Cf. Ref. Gukelberger et al. (2014). Note that, while there is a direct mapping between results for the half-filled attractive and repulsive models, the doped repulsive model maps into the attractive model in a magnetic field and vice versa Ho et al. (2009)..

Experimentally, the model (1) can be implemented using spin-dependent optical lattices Mandel et al. (2003) or using higher orbitals Müller et al. (2007); Wirth et al. (2010). In both cases, tuning the anisotropic hopping amplitudes of two species of atoms differently will lead to mismatched Fermi surfaces like those shown in Fig. 1(b). Moreover, by using magnetic gradient modulation, Ref. Jotzu et al. (2015) has created a continuously tunable state-dependent optical lattice. In the extreme anisotropic and large $U$ limit the model (2) approaches the 2D Ising model. The critical entropy per particle is $0.30647k_{B}$ Onsager (1944), which is likely to be within reach in current experiment Mazurenko et al. (2016). The presence of an additional trapping potential is likely to reduce the required entropy Jördens et al. (2010) because the metallic wings with incommensurate filling have higher entropy than the trap center.

Detection of the phase transition would be easiest via a spin-resolved density measurement. Forming of the antiferromagnetic pattern would also leave a signature in spin-resolved time-of-flight images. Testing our predictions would be an important step towards approaching the more exotic quantum phases at different filling and interaction strengths Gukelberger et al. (2014) and on different lattice geometries Wu (2008); Zhao and Liu (2008), which are at the limits of numerical control.

V Acknowledgment

We thank Y.-H. Liu for helpful discussions and Wei Tang for providing the critical entropy value of the 2D Ising model. JG is supported by the Swiss National Science Foundation, LW by the Ministry of Science and Technology of China under the Grant No.2016YFA0302400 and the start-up grant of IOP-CAS, and LP by FP7/ERC starting grant No. 306897. The DiagMC calculations were run on the Mammouth cluster of Université de Sherbrooke, provided by the Canadian Foundation for Innovation, the Ministère de l’Éducation des Loisirs et du Sport (Québec), Calcul Québec, and Compute Canada. The LCT-QMC calculations were run on the Tianhe-2 cluster of the National Supercomputer Center in Guangzhou. Simulations and data evaluation made use of the ALPS libraries B Bauer et al. (2011); A.F. Albuquerque et al. (2007).