Anisotropy crossover in the frustrated Hubbard model on four-chain cylinders

TL;DR

This study maps the phase diagram of a frustrated Hubbard model on four coupled chains, revealing various magnetic and metallic phases as interchain coupling varies, using advanced numerical methods.

Contribution

It introduces a hybrid real-momentum-space DMRG approach and compares it with variational cluster approximation for analyzing the phase diagram.

Findings

Identification of antiferromagnetic insulator at strong interactions

Discovery of incommensurate spin-density wave and metallic phases at intermediate couplings

Characterization of phases via magnetic, dielectric, and correlation properties

Abstract

Motivated by dimensional crossover in layered organic ${\kappa}$ salts, we determine the phase diagram of a system of four periodically coupled Hubbard chains with frustration at half filling as a function of the interchain hopping ${t_{\perp}/t}$ and interaction strength ${U/t}$ at a fixed ratio of frustration and interchain hopping ${t'/t_{\perp}=-0.5}$. We cover the range from the one-dimensional limit of uncoupled chains (${t_{\perp}/t=0.0}$) to the isotropic model (${t_{\perp}/t=1.0}$). For strong ${U/t}$, we find an antiferromagnetic insulator; in the weak-to-moderate-interaction regime, the phase diagram features quasi-one-dimensional antiferromagnetic behavior, an incommensurate spin-density wave, and a metallic phase as ${t_{\perp}/t}$ is increased. We characterize the phases through their magnetic ordering, dielectric response, and dominant static correlations. Our analysis is based primarily on a variant of the density-matrix renormalization-group algorithm based on an efficient hybrid-real-momentum-space formulation, in which we can treat relatively large lattices albeit of a limited width. This is complemented by a variational cluster approximation study with a cluster geometry corresponding to the cylindrical lattice allowing us to directly compare the two methods for this geometry. As an outlook, we make contact with work studying dimensional crossover in the full two-dimensional system.

Equations

$$\mathrm{H}=$$

$$-t^{\prime }\sum _{\text{llangle}\mathbf{r},{\mathbf{r}}^{\prime }\text{rrangle}\sigma }{c}_{\mathbf{r}\sigma }^{†}{c}_{{\mathbf{r}}^{\prime }\sigma }^{}+U\sum _{\mathbf{r}}{n}_{\mathbf{r}\uparrow }^{}{n}_{\mathbf{r}\downarrow }^{}$$

$$\epsilon (\mathbf{k})$$

$${\epsilon }_0(k_x)$$

$${\epsilon }_1(k_x)={\epsilon }_3(k_x)$$

$${\epsilon }_2(k_x)$$

$$\Omega [\Sigma ]=F[\Sigma ]-\mathrm{tr}\log (-G_0^{-1}+\Sigma ),$$

$$\Omega [\Sigma ]={\Omega }^{\prime }[\Sigma ]+\mathrm{tr}\log (-{G_0^{\prime }}^{-1}+\Sigma )-\mathrm{tr}\log (-G_0^{-1}+\Sigma ),$$

$$\frac{\delta \Omega [\Sigma ]}{\delta {\mathbf{t}}^{\prime }}=\frac{\delta \Omega [\Sigma ]}{\delta \Sigma }\frac{\delta \Sigma }{\delta {\mathbf{t}}^{\prime }}=0.$$

$$\overline{S}=\frac{1}{N}\sum _{i}S(i),$$

$$S_S(\mathbf{k})=\frac{1}{N}\sum _{\mathbf{r}{\mathbf{r}}^{\prime }}{e}^{\mathrm{i}\mathbf{k}(\mathbf{r}-{\mathbf{r}}^{\prime })}S(\mathbf{r},{\mathbf{r}}^{\prime })$$

$$S(\mathbf{r},{\mathbf{r}}^{\prime })=⟨S_z(\mathbf{r})S_z({\mathbf{r}}^{\prime })⟩-⟨S_z(\mathbf{r})⟩⟨S_z({\mathbf{r}}^{\prime })⟩$$

$$P=\sum _{x\sigma }x⟨\Psi |n_{x\sigma }|\Psi ⟩,$$

$${\mathcal{X}}_e=\frac{P}{E{L}_x},$$

$$S(\mathbf{r},{\mathbf{r}}^{\prime })$$

$$C(\mathbf{r},{\mathbf{r}}^{\prime })$$

$$D_{\mathrm{y}\mathrm{y}}(\mathbf{r},{\mathbf{r}}^{\prime })$$

$$S_z(x,y)$$

$$n(x,y)$$

$${\Delta }_{\mathrm{y}}^{†}(x,y)$$

Full text

I Introduction

How dimensional crossover between one-dimensional and two-dimensional behavior takes place in antiferromagnetic spin systems and in fermionic systems with Hubbard interactions is a question that has long been of interest Bourbonnais and Caron (1991); Boies et al. (1995); Affleck and Halperin (1996); Clarke and Strong (1997), especially in the context of ladder and anisotropic cuprate systems Dagotto and Rice (1996); Rice (1996); Nagata et al. (1998); Kojima et al. (2001) and in Bechgaard salts Giamarchi et al. (2004); Vescoli et al. (1998). Recent experimental investigations of the universality class of the Mott transition in layered organic charge-transfer salts Kagawa et al. (2005, 2009); Furukawa et al. (2015); Abdel-Jawad et al. (2015) have reopened the question of dimensional crossover, but with several new ingredients. The half filled systems are quasi-two-dimensional, anisotropic, and frustrated. The relative strength of these ingredients as well as of the relative interaction strength can be tuned experimentally by changing the molecular composition of the salts, by applying physical pressure, and by carrying out chemical substitution. In particular, a Mott metal-insulator transition (MIT) takes place which can be analyzed very precisely by measuring the conductivity as a function of temperature and ambient pressure Kanoda and Kato (2011).

On the theoretical side, a relatively simple model that is thought to capture the essential features of dimensional crossover is the two-dimensional, anisotropic, frustrated Hubbard model at half filling. Note that for the layered organic $\kappa$ salts, effective models have been postulated that are based on the anisotropic triangular lattice Kandpal et al. (2009); Nakamura et al. (2012), whereas, for the Bechgaard and Fabre salts, models based on the anisotropic square lattice have been postulated Jacko et al. (2013); Brown (2015). Let us nevertheless first consider what is known about the isotropic square-lattice case with isotropic frustration. Geometrical frustration can most easily be introduced by taking the single-band Hubbard model on a square lattice geometry and adding next-nearest neighbor hopping terms $t^{\prime}$ in either one (equivalent to the anisotropic triangular lattice) or in both diagonal directions. Phase diagrams for both cases have been obtained using the cluster dynamical mean-field theory (CDMFT) Kotliar et al. (2001) in Ref. Kyung and Tremblay (2006), the variational cluster approximation (VCA) Potthoff et al. (2003) in Refs. Nevidomskyy et al. (2008); Yamada et al. (2013), the path-integral renormalization group (PIRG) Kashima and Imada (2001) in Refs. Morita et al. (2002); Mizusaki and Imada (2006), the density-matrix embedding theory (DMET) Knizia and Chan (2012) in Ref. Zheng and Chan (2016), and the determinant quantum Monte Carlo method (DQMC) Blankenbecler et al. (1981) in Ref. Wang et al. (2016). Mott-insulating, metallic, antiferromagnetic, and superconducting phases appear depending on the level of frustration and the interaction strength. The antiferromagnetic phase further subdivides into regions of different magnetic ordering characterized by distinct ordering wave vectors Mizusaki and Imada (2006); Nevidomskyy et al. (2008); Yamada et al. (2013). Despite the enormous efforts made, completely conclusive phase diagrams have yet to be found, and the results obtained by the different methods often differ in the details. Open questions include the exact determination of the various magnetic orderings and the existence and nature of an intermediate phase between the insulating strong-coupling and the conducting weak-coupling regimes of the phase diagrams.

We now return to dimensional crossover, i.e., consider the effect of anisotropy in the couplings in the spatial directions. Here we want to make a clear distinction between dimensionality and frustration, i.e., be able to treat the two effects separately. Thus, in the following, we treat the square lattice geometry with next-nearest-neighbor hopping terms of equal strength in both diagonal directions. By tuning the interchain hopping $t_{\perp}$ the crossover between uncoupled chains over the “quasi-one-dimensional” case of weakly coupled chains to the isotropic two-dimensional model can be investigated Raczkowski and Assaad (2012, 2013); Raczkowski et al. (2015); Lenz et al. (2016). In the presence of frustration, the most straightforward parametrization is to keep the ratio between the value of the frustration, $t^{\prime}$, and $t_{\perp}$ fixed.

A direct motivation for the present work is Ref. Lenz et al. (2016), in which quantum critical behavior was identified using quantum cluster techniques at zero and finite temperature. In particular, the paramagnetic Mott transition was found to change from a discontinuous transition at $T_{c}>0$ at large interchain coupling to a continuous transition at zero temperature for small coupling strengths. In other words, the interchain hopping strength serves as a tuning parameter to realize a quantum phase transition when going below a critical value $t_{\perp\,c}$. This is an interesting starting point for further investigations of quantum critical behavior, which is seen in the layered organic charge-transfer salts, in particular, in the temperature-dependent behavior.

The study of Ref. Lenz et al. (2016), however, has a significant limitation in that only the paramagnetic case was treated, i.e., the possibility of magnetic ordering was explicitly excluded.

The goal of the present study is to go beyond this restriction by utilizing an unbiased numerical method, the density-matrix renormalization-group (DMRG), and by applying a variation of the VCA that allows for antiferromagnetic ordering. Since the variation of the DMRG that we use, which works in a hybrid of real and momentum space, most efficiently treats lattices with cylindrical geometry, we focus on them in both methods. Our study thus addresses the effect of antiferromagnetic fluctuations and benchmarks the VCA against the numerically exact hybrid-space DMRG results. Note that our aim here is not to specifically model the $\kappa$ salts; to do this one would have to treat an anisotropic triangular model. Here we investigate the effect of dimensionality systematically, which is believed to be—together with frustration—one of the driving mechanisms of the crossover in these systems. Due to the limitations of the entropy area law Eisert et al. (2010), it is exponentially difficult to make the cylinders wider. Thus, we generically consider cylinders of width 4 in the following, but also present results for cylinders of width 5 for selected values of the parameters. In particular, we study the ground-state properties as a function of the interaction strength and the anisotropy, fixing the frustrating hopping element ${t^{\prime}=-0.5\,t_{\perp}}$, where $t_{\perp}$ is the interchain hopping.

The low-energy properties of systems of coupled Hubbard chains can be treated in weak coupling using a field-theoretical treatment based on the renormalization group and bosonization Balents and Fisher (1996); Lin et al. (1997); Szirmai and Sólyom (2006); Arrigoni (1996). While Ref. Lin et al. (1997) does work out a detailed weak-coupling ground-state phase diagram for the four-chain case, it does not explicitly treat the half filled case, in which umklapp processes are relevant, and it does not include next-nearest-neighbor hopping, which is relevant for us here, in the band structure; therefore the results are of limited usefulness here. However, we will discuss relevant features of the four-chain band structure with next-nearest-neighbor hopping in the following, with a view to making contact with weak-coupling treatments.

The paper is organized as follows: the model is introduced in Sec. II, and details concerning the DMRG and VCA algorithms used are described in Secs. III and IV, respectively. The results are presented in Sec. V, which is subdivided into Sec. V.1, DMRG results, and Sec. V.2, VCA results. Section VI contains a detailed discussion of the results and their implications and Sec. VII concludes.

II Anisotropic Frustrated Hubbard model

We treat the anisotropic Hubbard model with frustration,

[TABLE]

where ${c^{\dagger}_{{\bf r}\,\sigma}}$, ${c^{\hphantom{\dagger}}_{{\bf r}\,\sigma}}$, and ${n^{\hphantom{\dagger}}_{{\bf r}\,\sigma}}={c^{\dagger}_{{\bf r}\,\sigma}}\,{c^{\hphantom{\dagger}}_{{\bf r}\,\sigma}}$ are the creation, annihilation, and density operators for lattice site ${{\bf r}=(x,y)}$ with spin ${\sigma\in\left\{\uparrow\,,\,\downarrow\right\}}$. As depicted in Fig. 1(a), the amplitudes $t$, ${t_{\perp}}$, and ${t^{\prime}}$ are for hopping in the longitudinal (i.e., intrachain), transverse (i.e., interchain), and diagonal directions; ${\langle{\bf r},{\bf r}^{\prime}\rangle_{\parallel}}$ and ${\langle{\bf r},{\bf r}^{\prime}\rangle_{\perp}}$ denote nearest-neighbor pairs of sites in the longitudinal and transverse directions, respectively, whereas ${\llangle{\bf r},{\bf r}^{\prime}\rrangle}$ denotes next-nearest-neighbor pairs, here along the diagonals. The dispersion relation for the infinite lattice then has the form

[TABLE]

leading to the noninteracting Fermi surface depicted for varying $t_{\perp}$ in Fig. 1(b). We study the model on lattices with periodic boundary conditions in the transverse direction and open boundary conditions in the longitudinal directions, i.e., on a cylinder geometry. Taking the lattice spacing to be unity, we denote the length in the longitudinal direction ${L_{x}}$, the width (in the transverse direction) ${L_{y}}$, and the number of sites $N=L_{x}\,L_{y}$.

Since our purpose is to study the influence of the interchain coupling for a sufficient degree of frustration, we take the intrachain hopping as the unit of energy, i.e., set ${t=1}$, and simultaneously vary the anisotropy, ${t_{\perp}}$, and the frustration, ${t^{\prime}}$, keeping the ratio of the two fixed to a fairly strong value: ${t^{\prime}/t_{\perp}=-0.5}$. We then tune ${t_{\perp}}$ between the limiting cases of uncoupled Hubbard chains at ${t_{\perp}=0}$ and the isotropic case at ${t_{\perp}=1}$. We work at half filling, where we note that the Hamiltonian (II) is symmetric with respect to the sign of ${t^{\prime}/t_{\perp}}$ if a particle-hole transformation is carried out, so that only the absolute value of ${t^{\prime}/t_{\perp}}$ is relevant. We take ${t^{\prime}/t_{\perp}}$ to be negative here in order to retain the form of the band structure in Fig. 2 and to maintain notational consistency with other work.

In the present study, we concentrate on cylinders with a width of four lattice sites, ${L_{y}=4}$. In the infinite-cylinder-length limit, the four bands ${\varepsilon_{j}(k_{x})}$ for transverse momentum ${k_{y}=j\,\pi/2}$ read

[TABLE]

In Fig. 2, the band structure is depicted for ${t_{\perp}=0.1}$, ${t_{\perp}=0.5}$, and ${t_{\perp}=0.9}$ with the Fermi level for the noninteracting case at half filling indicated. We denote the Fermi points within each band by ${\bf k}^{\rm F}_{j}$. For the case of weak interchain coupling, Fig. 2(a), the four bands approach each other, and the longitudinal Fermi momenta, $k^{\rm F}_{j\,x}$, converge towards ${\pm\pi/2}$ for ${t_{\perp}\rightarrow 0}$. In the ${t_{\perp}=0}$ limit, ${k^{\rm F}_{j\,x}=\pm\pi/2}$ for ${j=0,\ldots,3}$, the Fermi surface becomes nested, and intraband umklapp processes with momentum transfer of ${(2\,\pi,0)}$ are allowed within all four bands. For intermediate interchain coupling, Fig. 2(b), the Fermi points deviate from ${\pi/2}$, and intraband umklapp processes are not permitted. Still, interband umklapp processes utilizing all four bands with a total momentum transfer of ${\bf k}^{\rm F}_{0}+{\bf k}^{\rm F}_{1}+{\bf k}^{\rm F}_{2}+{\bf k}^{\rm F}_{3}={(2\,\pi,2\,\pi)}$ remain possible. Note that this is possible because at half filling ${\sum_{j=0}^{3}k^{\rm F}_{j\,x}=2\,\pi}$, assuming that all bands are partially filled. The scattering vectors involved, ${q_{a}:\,k^{\rm F}_{1}\rightarrow-k^{\rm F}_{3}}$ and ${q_{b}:\,k^{\rm F}_{0}\rightarrow-k^{\rm F}_{2}}$, are depicted as gray dashed and solid arrows in Fig. 2(b). Note that other interband umklapp processes involving all four Fermi bands are also allowed; however, we only show the scattering vectors $q_{a}$ and $q_{b}$ since our DMRG results (see Sec. V.1) indicate that the relevant processes have transverse momentum $\pi$. For strong interchain coupling, Fig. 2(c), the ${k_{y}=0}$ band becomes flat and completely filled, only three bands remain active, and no umklapp processes are allowed.

The transition between the latter scenarios shown in Figs. 2(b) and 2(c), i.e., the exact point where ${\varepsilon_{0}(k_{x})}$ has a maximum at ${k^{\rm F}_{0\,x}=\pi}$, marks a Van Hove singularity. In the two-dimensional case, this transition corresponds to the point where the noninteracting Fermi surface undergoes a topological change, changing from an open to a closed surface, as depicted in Fig. 1(b). Note that the critical interchain hopping strength for the width-$4$ cylinder is ${t_{\perp}\approx 0.707}$, somewhat displaced from that of the two-dimensional case, where the singularity is at ${t_{\perp}\approx 0.622}$.

III DMRG methods

We apply the DMRG within a hybrid-space representation that is composed of a real-space representation in the longitudinal and a momentum-space representation in transverse direction, as was recently introduced for the fermionic Hofstadter model Motruk et al. (2016) and the two-dimensional Hubbard model Ehlers et al. (2017). In this representation, the transverse momentum ${k_{y}}$ is conserved; this can be utilized to speed the DMRG algorithm up significantly. The retention of a real-space representation in the longer spatial direction makes it possible to avoid the undesirable volume-law scaling of the entropy Ehlers et al. (2015) that severely restricts the applicability of the pure momentum-space DMRG Xiang (1996). Our implementation Ehlers et al. (2017) is based on the matrix product state (MPS) formulation of the DMRG Schollwöck (2011). The utilization of conserved transverse momentum, total charge, and total spin quantum numbers enables us to push the dimension of the virtual bonds of the MPS, $m$, to up to ${35\,000}$ states; a detailed analysis of the performance of the hybrid-space DMRG for the two-dimensional Hubbard model is contained in Ref. Ehlers et al. (2017).

Despite the enhanced performance of the hybrid-space formulation, the DMRG can still suffer from convergence difficulties. First, the convergence of the DMRG, both in the real-space and in the hybrid-space formulations, is exponentially costly in system width due to the entropy area law Eisert et al. (2010). Therefore, we study the system at a relative small width of ${L_{y}=4}$, i.e., four coupled chains. For isolated parameter points, we will show results for ${L_{y}=5}$ in order to roughly determine how the behavior changes as the number of chains is increased. Calculations for ${L_{y}=5}$, for which convergence to the qualitatively consistent ground state only occurs for ${m\gtrsim 20\,000}$, are prohibitively expensive to carry out for the full parameter space treated here.

Second, the algorithm sometimes gets trapped in a metastable state, which then, for calculations of finite accuracy, converges as the apparent ground state of the system. Such behavior is often found for parameter values close to a discontinuous transition. In order to minimize such problems and define the phase boundaries as accurately as possible, we use a variation of the standard DMRG ground-state search: we start DMRG calculations for two sets of model parameters located on either side of, but sufficiently far from, a discontinuous phase transition. We initialize both calculations by following the standard ground-state search procedure, in which the model parameters are kept fixed and $m$ is increased during the initial sweeps. Once sufficient convergence has been attained for a particular starting point, we fix $m$, continue sweeping, and alter the model parameters in small steps towards the other phase. As the system goes through a discontinuous phase transition, we typically obtain a hysteresis effect and can then track the ground states of both phases into the parameter regime of the other phase. By determining the exact point at which the energies of the two states cross, we obtain the phase boundary with maximal accuracy.

As an example, Fig. 3 shows calculations at fixed ${t_{\perp}=0.7}$ as a function of ${U}$. In this case, we find two discontinuous phase transitions at the critical interaction strengths ${U_{c\,1}}$ and ${U_{c\,2}}$. The first two panels, Figs. 3(a) and 3(b), show the energy-level crossings for the first and second transitions, respectively. These crossings correspond to finite jumps in the double occupancy ${D=\sum_{{\bf r}}\langle n_{{\bf r}\,\uparrow}n_{{\bf r}\,\downarrow}\rangle}$, as can be seen in Fig. 3(c). This example shows that the hysteresis effect is present in the double occupancy as well as in the energy. Thus, we can determine the point of the transition and the height of the jump in the double occupancy accurately. Results from standard DMRG calculations, carried out separately for every distinct value of $U$, are depicted as blue circles for comparison.

The procedure described above has both an advantage and disadvantages: the advantage is that convergence problems close to the phase boundary are diminished because the critical initial phase of the DMRG algorithm takes place in a more stable region of the parameter space. One disadvantage is that, since $m$ is fixed during the second phase of the calculation, the results cannot be extrapolated in the truncation error without further effort. Another is that the CPU time needed for a single calculation is comparably high because a large number of sweeps in which the maximum $m$ is kept must be carried out. In our calculations, we have used a fixed MPS bond dimension of ${m=10\,000}$ after the initial sweeps, which was large enough to converge the calculations and keep runtimes at an acceptable level.

In order to validate our results qualitatively, we have carried out additional calculations in which we used the standard DMRG ground-state search at all points of parameter space. For these calculations, we have used a maximum bond dimension between ${m=20\,000}$ and ${m=35\,000}$, have varied the chain length, and have extrapolated the results in the truncation error.

IV variational cluster approximation

The variational cluster approximation (VCA) Potthoff et al. (2003) is a quantum cluster technique Maier et al. (2005), which is used to study strongly correlated electron systems with local interactions and to investigate phases with broken symmetries Dahnken et al. (2004). Its validity has been successfully demonstrated for the one-dimensional Hubbard model Potthoff et al. (2003); Balzer et al. (2008), and it has been applied to two-dimensional Hubbard systems Dahnken et al. (2004); Aichhorn et al. (2006), allowing for both, the investigation of magnetically ordered phases Dahnken et al. (2004); Nevidomskyy et al. (2008); Yamada et al. (2013) and of the Mott transition Balzer et al. (2009); Lenz et al. (2016).

The VCA is based on Potthoff’s self-energy functional theory Potthoff (2003a, b, 2012), in which the grand-canonical potential of the infinite system is written as a functional of the self-energy $\Sigma$,

[TABLE]

where $F[\Sigma]$ denotes the Legendre-transformed Luttinger-Ward functional Luttinger and Ward (1960); Potthoff (2006), and ${G_{0}=(\omega+\mu-\mathbf{t})^{-1}}$ is the noninteracting single-particle Green’s function. At the physical self-energy $\Sigma_{\text{phys}}$, this functional is stationary as a function of the self-energy, i.e., ${\delta\Omega[\Sigma]/\delta\Sigma|_{\Sigma_{\text{phys}}}=0}$. Since ${F[\Sigma]}$ is usually not known, one has to resort to approximations to determine the self-energy so that this stationary condition is fulfilled.

In the VCA, a reference system consisting of decoupled clusters is used to calculate the self-energy functional. This reference system, described by a Hamiltonian ${H^{\prime}}$, has the same interaction terms as the original system and the cluster self-energy can be computed exactly, for example, using exact diagonalization. Since ${F[\Sigma]}$ is universal in the sense that it only depends on the interaction part of the Hamiltonian, both the original and the reference systems have the same Luttinger-Ward functional. Therefore, the self-energy functional can be expressed in terms of the reference system as

[TABLE]

where all cluster quantities are denoted by a prime. Through variation of the one-body parameters ${\mathbf{t}^{\prime}}$ of the reference system, the stationary point of the functional can be determined because

[TABLE]

The approximation of the VCA thus lies in restricting the variational space of all possible self-energies in Eq. (5) to a limited set of cluster-representable self-energies $\Sigma^{\prime}$.

In order to account for phases with broken symmetry such as an antiferromagnetic phase, additional Weiss fields are added to the cluster Hamiltonian, and their field strength is used as one of the variational parameters to find the stationary point of ${\Omega[\Sigma]}$. In this way, short-range correlations within the cluster are treated exactly, but longer-ranged correlations are treated on a mean-field level Aichhorn et al. (2006). In order to properly describe the discontinuous (paramagnetic) Mott transition, it was shown that it is necessary to include noninteracting bath sites in the reference cluster Balzer et al. (2009). This is permitted within self-energy functional theory because additional couplings to noninteracting sites do not change the Luttinger-Ward functional, and the clusters thus are a valid reference system.

V Results

In the following, we present our results for the DMRG, Sec. V.1, and the VCA, Sec. V.2. The DMRG results cover a detailed phase diagram in the $t_{\perp}$ and $U$ plane and determine the static magnetic and electric properties of the phases. For the VCA, we focus on four different values of $t_{\perp}$, for which we determine the position and nature of the metal-insulator transitions and compare them to the DMRG results.

V.1 DMRG results

We now use the tracking method described in Sec. III to locate phase-transition lines in the $t_{\perp}$–$U$ plane for cylinders with ${16{\times}4}$ sites. We find two phase-transition lines in the regime between ${t_{\perp}=0}$ and ${t_{\perp}=1}$, both of which are due to level crossing scenarios in the ground state for the entire range of $t_{\perp}$. The resulting phase diagram is depicted in Fig. 4. Note that the transition lines are determined only for the ${16{\times}4}$ lattice. In order to unambiguously determine the order of the phase transitions, a finite-size scaling analysis would be needed. Due to the high computational costs, we refrain from doing this here. However, as discussed further below in this section, we perform a finite-size scaling analysis for the electrical susceptibility, and investigate in detail the behavior of the spin structure factor. The behaviors of both quantities support that the lines of level crossings depicted in Fig. 4 will indeed remain discontinuous phase transitions in the infinite-length limit.

The phase diagram in Fig. 4(a) features a metallic phase, an antiferromagnetic phase (AF), a quasi-one-dimensional antiferromagnetic phase (Q1D-AF), and an incommensurate spin-density-wave phase (ISDW) that separates the metallic and the antiferromagnetic phases. Both transition lines are characterized by finite jumps in the double occupancy, ${\Delta D}$, Fig. 4(b), along the lines indicating the presence of discontinuous phase transitions. The metallic phase extends from ${U=0.0}$ to ${U\approx 6.0}$ in the isotropic case ${({t_{\perp}=1})}$, becomes narrower with decreasing ${t_{\perp}}$, and vanishes at ${t_{\perp}\approx 0.6}$. Note that close to the latter point, at ${t_{\perp}\approx 0.707}$, a Van Hove singularity marks the transition between three and four active bands in the band structure for width-$4$ cylinders; see Sec. II. The ISDW phase bounds the metallic phase and is present for transverse hopping strength ${0.4\lesssim t_{\perp}\leq 1.0}$. These DMRG results for the phase diagram are compared to transitions found by the VCA at ${t_{\perp}=0.5}$, $0.8$, and $1.0$. As can be seen in Fig. 4, discontinuous transitions from a metallic to an AF phase are found for large values of $t_{\perp}$, which agree very well with the locations of the transitions identified by the DMRG. At smaller values of $t_{\perp}$, a continuous transition is seen, which is difficult to capture using the DMRG. These results and the comparison to the DMRG phase diagram will be discussed in detail in the following and in Sec. V.2.

Figure. 5 shows the behavior of the double occupancy and the block entropy at the two transitions between the metallic and the ISDW phases and between the ISDW and AF phases at fixed ${t_{\perp}=0.7}$ for width-4 cylinders with different lengths and for the width-5 cylinder of length 16. In Fig. 5(a), the double occupancy is shifted by ${0.03U}$ to compensate for the linear growth in $U$ and to emphasize the jumps between the phases. For width-4 cylinders, the jump between the ISDW and AF phases seems to be stable as a function of the system length. At the transition between the ISDW and the metallic phase, the jump is present for all system lengths, but shrinks somewhat with increasing $L_{x}$. Thus, from the behavior of the double occupancy alone, we cannot rule out that this transition becomes continuous in the infinite-length limit. In order to confirm this, one would need to perform a scaling analysis, which cannot be done reliably with the data set available.

The discontinuous behavior is also marked by a discontinuity in the entropy. Figure 5(b) shows the averaged von Neumann entropy,

[TABLE]

where ${S(i)=-\mathrm{tr}(\rho_{i}\ln\rho_{i})}$ is the von Neumann, entanglement entropy Amico et al. (2008) calculated for a bipartite partitioning of the system at the $i$th MPS bond and $\rho_{i}$ the corresponding reduced density matrix of one of the subsystems. (Note that the state obtained within the DMRG is a pure state, usually an approximation to the ground state of the system.) Again, for width-4 cylinders, the jump is stable for the transition between the ISDW and AF phases, but the data does not allow for a definite statement on the nature of the transition between the metallic and ISDW phases.

For the ${16{\times}5}$ lattice, the signatures of both transitions are evident in the double occupancy and in the entropy. For width 5, a scaling in system length is not possible due to the high computational costs.

We classify the magnetic ordering of the AF, Q1D-AF, and ISDW phases by examining the static structure factor

[TABLE]

of the static spin correlations

[TABLE]

with ${S_{z}(x,y)=}\,{c^{\dagger}_{x\,y\,\uparrow}\,c^{\vphantom{\dagger}}_{x\,y\,\uparrow}}-{c^{\dagger}_{x\,y\,\downarrow}\,c^{\vphantom{\dagger}}_{x\,y\,\downarrow}}$. Since we have open boundary conditions in the $x$ direction, we use particle-in-a-box eigenmodes for the transformation in that direction to approximate momentum modes, as in Ref. Hager et al. (2005). Figure 6 depicts the structure factor with panel (a) showing the result for a single chain and panels (b)–(e) following a cut through the phase diagram at constant ${U=4.0}$ with increasing interchain coupling. For weak ${t_{\perp}=0.075}$, Fig. 6(b), ${S_{S}({\bf k})}$ has one peak at ${k_{x}=\pi}$ in all four transverse momentum branches, but the form of ${S_{S}({\bf k})}$ only varies weakly with transverse momentum, $k_{y}$. The shape of ${S_{S}({\bf k})}$ in the longitudinal direction has almost exactly the same form as that of the one-dimensional case, Fig. 6(a), especially for the ${k_{y}=\pi}$ branch. For larger ${t_{\perp}}$, the two-dimensional antiferromagnetic phase is characterized by a very strong peak in ${S_{S}({\bf k})}$ at ${{\bf k}=(\pi,\pi)}$, which contrasts with much weaker peaks at ${k_{x}=\pi}$ in the other $k_{y}$ channels, as can be seen in Fig. 6(c) for a representative point at ${t_{\perp}=0.5}$.

The $U$-dependence of the height of the peak at ${{\bf k}=(\pi,\pi)}$ is displayed in Fig. 7. For small to moderate ${t_{\perp}\lesssim 0.3}$, the height of the peak at ${{\bf k}=(\pi,\pi)}$ grows as $U$ is increased. For all ${t_{\perp}\lesssim 0.3}$, the growth follows that of the one-dimensional case, ${t_{\perp}=0.0}$, up to ${U\approx 1.5}$. Above ${U\approx 1.5}$, the peak heights continue to follow the one-dimensional case for smaller $t_{\perp}$ (${t_{\perp}\lesssim 0.1}$), but for ${0.1\lesssim t_{\perp}\lesssim 0.3}$, the peak height grows significantly more strongly with $U$ than the one-dimensional case. This is an indication of the crossover to a two-dimensional antiferromagnetic phase. The crossover line determined in this way is marked by the blue dashed line in Fig. 4(a).

The ISDW phase, Fig. 6(d), is characterized by an incommensurate ordering wave vector ${{\bf Q}^{*}=(\pm q^{*},\pi)}$, where the longitudinal wave vector of the peak, $q^{*}$, depends on ${U}$ and ${t_{\perp}}$. This incommensurate structure corresponds to antiferromagnetic spin correlations in real space with a sine-shaped modulation in the longitudinal direction, where the envelope has a wavelength of ${2\,\pi/(\pi-q^{*})}$. Figure 8 shows plots of ${S_{S}(k_{x},\pi)}$ within the ISDW phase for varying $L_{x}$, ${t_{\perp}}$, and ${U}$. In the metallic phase, depicted for ${t_{\perp}=0.9}$ and ${U=4.0}$ in Fig. 6(e), the structure factor ${S_{S}({\bf k})}$ lacks a distinct peak in all transverse momentum sectors.

Figure 8(a) shows that the wave vector of the incommensurate peak, $q^{*}$, is stable as a function of system length $L_{x}$ for $L_{x}$ larger than the corresponding wavelength. All but one of the curves shown are calculated with open boundary conditions in the $x$ direction. The additional curve is calculated on a $16\times 4$ lattice with periodic boundary conditions (the largest periodic lattice size for which we obtained good convergence). While the finite-size effects and the placement of momentum points are different for periodic than for open boundary conditions (in particular, the momentum points are more widely spaced), the values of ${S_{S}(k_{x},\pi)}$ and the approximate position of the peak are consistent with the open-boundary-condition results, showing that the incommensurate peak structure is not an artifact of the boundary conditions.

As ${t_{\perp}}$ and ${U}$ are varied, Figs. 8(b) and 8(c), $q^{*}$ changes, ranging between ${3\,\pi/4}$ and $\pi$, with $q^{*}$ moving towards $\pi$, i.e., with the wavelength of the modulation becoming longer, as the AF phase is approached, and $q^{*}$ moving towards ${3\,\pi/4}$ (shorter modulation wavelength) as the metallic phase is approached. Note that this movement of the peak position, $q^{*}$, can only be observed for systems of sufficient length; if $L_{x}$ is too small, the wavelength can get locked in at a particular discrete value due to the boundary conditions. This effect causes the large discrepancy in the double occupancy in Fig. 5(a) between system lengths $16$, $24$, and $32$.

In order to roughly gauge the effect of cylinder width, we display the spin structure factor for the ${16{\times}5}$ cylinder in Fig. 8(d). Note that the antiferromagnetic ordering is frustrated by the odd system width, and, for the periodic transverse boundary conditions used here, transverse momentum ${k_{y}=\pi}$ is not available, so that we take the closest point, ${k_{y}=4\,\pi/5}$. As can be seen, the incommensurate peak characterizing the ISDW phase is still present, but less distinct for ${U=3.4\text{--}4.2}$. At ${U=4.4}$, the peak is at ${k_{x}=\pi}$, as expected for the antiferromagnetic phase, and, at ${U=3.2}$, the form is as expected for the metallic phase. Note that the peaks are also weakened due to the small system length; compare to Fig. 8(a).

The static spin structure factor can also be used to characterize the transitions between the ISDW and the metallic and between the metallic and the AF phases. The peak height and position changes smoothly within the ISDW phase, but changes abruptly at the transition points toward the metallic and antiferromagnetic phases; compare Fig. 6(d) with Figs. 6(e) and 6(c), respectively. These discontinuous transitions are most evident in the behavior of the peak height of the structure factor, as plotted for width-4 cylinders of different lengths in Fig. 9. It can be seen that, although the position of the transition between the metallic and the ISDW phases changes, the size of the jump in the peak height at the transition remains stable as a function of $L_{x}$. For the transition between the AF and the ISDW phase, the position of the peak is stable and the height of the jump increases with system length. In the AF phase, the peak height increases with cylinder length. However, we note that the scaling is quasi-one-dimensional because we scale with system length only. In one dimension, long-range antiferromagnetic order cannot occur; the correlations can decay at most with a power law Giamarchi (2003). Increase of the peak height with system size will only occur if the power-law correlations fall off more slowly than $1/x$ (with logarithmic scaling occurring at $1/x$). This picture is consistent with the scaling of the AF peak in Fig. 9, which seems to grow sublinearly with cylinder length. In the ISDW phase, there is no significant scaling of the peak height with cylinder length, especially for the two larger $L_{x}$ values. Thus, if power-law decay is present, it must be more rapid than $1/x$, and we cannot differentiate between such a decay and exponential decay with a relatively long decay length, i.e., between a weakly critical and a weakly gapped phase.

For sufficiently weak interaction strength, one might expect that longitudinal wave vector $q^{*}$ of the peak of the spin structure factor would be determined by wave vectors characteristic of relevant scattering processes in a weak-coupling picture. Such scattering wave vectors are determined by the ${U=0}$ band structure, Eq. (2), as described in Sec. II, and are depicted in Fig. 2. Figure 10 shows the position of the maximum of ${S_{S}(k_{x},\pi)}$ as a function of $t_{\perp}$ for ${U=1.0}$ and ${U=1.6}$. For small $t_{\perp}$, the maximum is at transverse momentum $\pi$. Within the ISDW phase, the single maximum splits into two maxima symmetric to $\pi$, which move towards ${(1\pm 1/4)\pi}$ as $t_{\perp}$ is increased, and finally, in the metallic phase, the positions of the maxima remain stable at around ${(1\pm 1/4)\pi}$. For comparison, scattering wave vectors ${q_{a}=k^{\rm F}_{0\,x}+k^{\rm F}_{2\,x}}$ and ${q_{b}=k^{\rm F}_{1\,x}+k^{\rm F}_{3\,x}}$, defined as in Fig. 2(b), are depicted. As can be seen, the DMRG results for both ${U=1}$ and ${U=1.6}$ are in good agreement with the scattering wave vectors given by the ${U=0}$ band structure, including the signature of the Van Hove singularity at ${t_{\perp}\approx 0.7}$.

Since translational invariance is broken in the longitudinal direction by the boundary conditions, the structure of the phases can be seen in the behavior of local quantities, in particular, the deviation of the local charge density, ${\Delta\,n(x)}={L_{y}^{-1}\sum_{y}\langle n(x,y)\rangle-1}$, from its average value ${n=1.0}$ and the deviation of the local $z$-spin moment squared, ${\Delta\,S_{z}^{2}(x)}={L_{y}^{-1}\sum_{y}\langle S_{z}^{2}(x,y)\rangle}-{S_{z}^{2}(x,y)}$, from its average value ${S_{z}^{2}=N^{-1}\sum_{\bf r}\langle S_{z}^{2}(\bf r)\rangle}$ as a function of $x$, the position in the longitudinal direction. In the AF phase, as depicted in Fig. 11(a), spatial fluctuations in both the charge and spin densities are strongly suppressed, with the end effects from the open boundaries damping out very quickly. Thus, the dominant antiferromagnetic correlations strongly suppress all spatial fluctuations in the charge and spin density. Note that the average value ${S_{z}^{2}=0.825}$ of the $z$-spin moment squared is relatively large, i.e., the polarization of the local spin moment suppresses fluctuations. In the ISDW phase, Figs. 11(b) and 11(c), ${\Delta\,{S_{z}^{2}}}$ exhibits a wave pattern consistent with the incommensurate wave vector $q^{*}$ found in the spin structure factor. The rapidity of the decay of the fluctuations away from the ends depends on the location within the ISDW phase; for larger ${U}$ the SDW is more pronounced, and the modulations reach more deeply into the bulk of the system. Note that the charge-density fluctuations are comparably small in the ISDW phase, albeit not as small as in the AF phase. In the metallic phase, Fig. 11(d), fluctuations in the spin density are relatively small, but still present, and fluctuations in the charge density are significantly larger than in the other phases, and have a high-frequency component that extends to the center of the lattice.

In order to investigate the metallic nature of all three phases, we have calculated the electric susceptibility $\mathcal{X}_{e}$ for different system lengths $L_{x}$. By applying an electric field $E$ in the longitudinal direction, we can measure the polarization

[TABLE]

which, in turn, allows us to obtain the electrical susceptibility

[TABLE]

assuming that $E$ is chosen small enough so that the response is in the linear regime. (See Ref. Aebischer et al. (2001) for a detailed description of the method.) Figure 12 depicts the scaling of $\mathcal{X}_{e}$ with the inverse system length ${1/L_{x}}$ for different ${U}$ and ${t_{\perp}=0.7}$. In the metallic phase, $\mathcal{X}_{e}$ grows proportionally to ${L_{x}^{2}}$, as can be seen from its approximately linear behavior with a slope of $-2$ on the log-log scale. Such a scaling is characteristic of quasi-one-dimensional metallic behavior Aebischer et al. (2001). In contrast, in the AF phase $\mathcal{X}_{e}$ clearly saturates as $L_{x}$ becomes larger (i.e., ${1/L_{x}\to 0}$). In the ISDW phase, $\mathcal{X}_{e}$ is significantly larger than in the AF phase, but still saturates as $L_{x}$ becomes larger, thus indicating that the ISDW phase is, indeed, insulating. Note that, since the determination of the electrical susceptibility becomes progressively more difficult with decreasing interaction strength, this analysis cannot be continued along the transition lines towards smaller $t_{\perp}$.

Finally, we address the question of whether $d$-wave pairing correlations are enhanced in any part of the phase diagram. Since the system is quasi-one-dimensional, we, in general, expect at most power-law decay of the correlation functions with distance and must compare the strength of the decay in the different channels. We thus calculate the equal-time spin, charge, and pair correlation functions,

[TABLE]

with the local spin density, local charge density, and pair creation operators

[TABLE]

Figure 13 shows the decay of the correlation functions as a function of distance in the longitudinal direction, ${l_{x}=|x-x^{\prime}|}$. In the AF and ISDW phases, Figs. 13(a) and 13(b), respectively, the spin correlations are markedly dominant, as expected. In the ISDW phase, the modulations of the spin correlations are compatible with the wave pattern seen in the local spin fluctuations in Fig. 11(b). In the metallic phase, Fig. 13(c), all correlations decay at approximately the same rate. Thus, there is no clearly dominant correlation. Note, however, that we have not systematically investigated the dependence of the relative strength of the correlations on the model parameters within the metallic phase.

V.2 VCA results

In this section we apply the VCA to the Hamiltonian of Eq. (II) on a four-leg tube geometry, Fig. 14(a), to obtain results that can be compared with the hybrid-space DMRG results on width-4 cylinders. Our reference system consists of a ${2{\times}4}$ cluster with periodic boundary conditions in the $y$ direction, which is coupled to four bath sites, as depicted in Fig. 14(b). In order to find the stationary point of the self-energy functional, we use the hybridization $V$ between interacting cluster sites and bath sites 111 As shown in Ref. Balzer et al. (2009), using the hybridization $V$ as a variational parameter is essential to allow for the coexistence of insulator and metal in a region around the (paramagnetic) MIT. The insulating and metallic solutions then differ in $V$., the chemical potential $\mu^{\prime}$ of the cluster 222 Only when using $\mu^{\prime}$ as a variational parameter, the determination of the electron density $n$ is thermodynamically stable: The calculation of $n$ via the Green’s function and via a functional derivative of the self-energy functional leads to the same result Aichhorn et al. (2006)., the chemical potential $\mu$ of the system 333 Following the procedure introduced in Ref. Balzer and Potthoff (2010), we Legendre transformed the self-energy functional with respect to $\mu$ to be able to set a electron density to that of half filling. The corresponding chemical potential $\mu$ can then be determined via the variational principle., and the strength of an antiferromagnetic Weiss field on the bath sites 444 We have checked that choosing a Weiss field that acts on the cluster sites instead leads to qualitatively same results for the AF phase at ${t_{\perp}=1}$. Choosing the Weiss field to act on the bath sites is in the style of cluster dynamical mean-field theory, where the antiferromagnetism also enters via the bath Kyung and Tremblay (2006). as variational parameters.

Figure 15 depicts the energy and double occupancy for the different solutions found by the VCA for interchain couplings ${t_{\perp}=0.5}$, ${t_{\perp}=0.8}$, and ${t_{\perp}=1.0}$, together with DMRG results for ${16{\times}4}$ cylinders. Since the VCA only captures correlations within the cluster exactly, the ISDW phase found in the DMRG cannot be seen within the VCA. Within the VCA, at interaction strengths below a critical value ${U^{\phantom{*}}_{c}}$, a paramagnetic (PM) metal is realized; at larger interaction strengths, the system is found to be an antiferromagnetic insulator.

A paramagnetic insulating solution is found at large interaction strengths (black line in Fig. 15), which, however, is always higher in energy than the antiferromagnetic insulator and is therefore not realized. Note that we find a crossing of the energies of the two paramagnetic solutions at $U_{c}^{*}$, which is marked by the black dashed lines in Figs. 15(b) and 15(c). The crossing is associated with a jump in the double occupancy at $U_{c}^{*}$. This, together with the substantial coexistence region of the paramagnetic insulating and metallic solutions around $U_{c}^{*}$, shows that the relatively small number of bath sites used in our cluster is indeed sufficient to capture the discontinuous transition in the paramagnetic channel. In contrast, using the VCA without bath sites not only leads to a continuous paramagnetic metal-insulator transition with a much smaller $U_{c}^{*}$, but also yields a ${U^{\phantom{*}}_{c}}$ that is significantly smaller (e.g, ${U_{c}=3.8}$ for ${t_{\perp}=0.8}$).

Comparing the VCA data to the DMRG results, we find that, for ${t_{\perp}=0.8}$ and ${t_{\perp}=1.0}$, the energies of the antiferromagnetic insulator in the VCA and in the DMRG calculation on the ${16{\times}4}$ cylinder nearly coincide, and the agreement between the double occupancies in this phase is also very good. The phase transition between the antiferromagnetic insulator and the paramagnetic metal in the VCA occurs within the ISDW phase seen by the DMRG, i.e., ${U^{\rm DMRG}_{c\,1}<U^{\phantom{*}}_{c}<U^{\rm DMRG}_{c\,2}}$. This is consistent to within the limitations of the methods, since, as discussed before, the intermediate ISDW phase cannot be found by the VCA for the cluster sizes treated. Nevertheless, it is interesting to see that both a transition from an AF insulating phase to a PM metal as well as discontinuous behavior are found by both methods. For smaller interaction strengths, the energies of the metallic phases within the VCA and the DMRG are again closer to each other. Note that the degree of numerical agreement has been examined here for particular (different) fixed sizes of the DMRG lattice and the VCA cluster; a more complete analysis would compare results only after separate finite-size scaling for each method.

The double occupancy of the paramagnetic metal is a little higher than in the DMRG. This overestimation of the strength of the paramagnetic phase might be due to the fact that two cluster sites must share one bath site: For the isotropic unfrustrated two-dimensional Hubbard model, Balzer et al. found that either the metallic or the insulating phase is preferred, depending on whether an even or an odd number of bath sites is coupled to each cluster site Balzer et al. (2009). We expect that the choice of the number of bath sites would have a similar effect in the present case. In order to check this assumption, one would have to treat a cluster that includes at least eight bath sites. The fact that we are using a band-Lanczos solver at zero temperature and a variational space of dimension 4 means that treating the 16-site cluster would, unfortunately, have a much higher computational cost, precluding a systematic investigation of this question.

Whereas the DMRG and the VCA are in good quantitative agreement for the system or cluster sizes used, respectively, for large values of $t_{\perp}$, qualitative as well as quantitative differences are present at intermediate transverse hopping, ${t_{\perp}=0.5}$. At relatively strong coupling, the VCA ground state remains an antiferromagnetic insulator, just as in the DMRG analysis. A comparison of the energies of the AF solution of the VCA and the (AF) DMRG ground-state energy down to the transition to a ISDW again yields good agreement. However, while the DMRG finds the ISDW phase for moderate to low $U$, the VCA finds a paramagnetic metal at small $U$, yielding a metal-insulator transition at ${U_{c}\approx 2.9}$. This transition point is comparable to the value of ${U\approx 2.3}$ for the transition from the AF insulator to the ISDW obtained by the DMRG. Since the VCA on the treated cluster size is fundamentally unable to find an ISDW phase, the AF-to-ISDW transition cannot be seen. Nevertheless, it is interesting to see that the VCA identifies a phase transition in a similar parameter region as the DMRG.

Finally, Fig. 16 shows how the situation changes for an even smaller value of the interchain hopping, ${t_{\perp}=0.2}$. In the VCA, an antiferromagnetic solution is still realized for ${U>1.5}$, Fig. 16(a), but, in contrast to the previous cases at larger ${t_{\perp}}$, the energies of the AF insulator and PM metal approach each other smoothly as $U$ is decreased, Fig. 16(b), and a level crossing cannot be observed. At ${U=1.5}$, both solutions have almost the same energy ${(|E_{\text{PM-M}}-E_{\text{AF}}|<10^{-5})}$, and the double occupancy exhibits no measurable jump at the transition point, as can be seen in Fig. 16(c). For ${U<1.5}$, an AF solution cannot be found within the VCA. For the same $t_{\perp}$, the DMRG results show that the AF ordering in the form of a distinct peak in ${S_{S}({\bf k})}$ at ${{\bf k}=(\pi,\pi)}$ smoothly vanishes around the same interaction strength, ${U\approx 1.5}$; see Figs. 4(a) and 7. Since, within the VCA, it is not possible to introduce a Weiss field that captures the quasi-one-dimensional order without introducing some kind of order in $y$ direction, the Q1D-AF cannot be reproduced. Still, the VCA finds a vanishing AF order and consequently a crossover to a paramagnetic phase.

VI Discussion

We now discuss the relation between our results and those of other studies. For systems of a finite number of coupled Hubbard chains with varying interchain coupling strength, a number of calculations using bosonization and renormalization group techniques in the weak-coupling limit have been carried out Balents and Fisher (1996); Lin et al. (1997); Szirmai and Sólyom (2006); Arrigoni (1996). The case of four periodically coupled chains was treated in detail in Ref. Lin et al. (1997). However, the half filled system, for which additional umklapp processes are relevant, was not handled explicitly. In addition, next-nearest-neighbor hopping, which generates frustration and is essential to our study here, was also not included. Thus, we know of no weak-coupling studies that would be directly applicable to the system studied here. Such calculations would be extremely useful, in particular, in clarifying the role of the umklapp processes within the ISDW phase and in determining the phases present in the ${U\rightarrow 0}$ limit. However, such a study would go beyond the scope of the present work.

For the isotropic case, a number of studies are available that investigate the phase diagram of the two-dimensional Hubbard model as a function of the level of frustration. For the value of the frustration treated here, ${{t^{\prime}}=-0.5}$, VCA calculations of Nevidomskyy et al. Nevidomskyy et al. (2008), in which no bath sites were included, found a transition from a superconducting phase to a ${(\pi,\pi)}$ antiferromagnet, with an intermediate region of phase coexistence for ${3.0\leq U\leq 4.0}$. For the same level of frustration, Mozusaki and Imada reported a transition from a metallic phase to an antiferromagnetic insulating phase separated by an insulating, nonmagnetic, intermediate phase between ${U\approx 4.0}$ and ${U\approx 6.0}$ in their PIRG study Mizusaki and Imada (2006). In another study using the VCA without bath sites, Yamada et al. found either a transition between a metallic phase and an insulating antiferromagnetic or between a metallic phase and an insulating paramagnetic phase Yamada et al. (2013), depending on the size of the cluster considered. The above studies also treated stronger levels of frustration, at which they found different antiferromagnetic phases: for ${t^{\prime}\gtrsim 0.8}$, all three studies find a “striped” phase, ${{\bf q}=(\pi,0)}$, and the PIRG and, depending on the cluster size, the VCA, find antiferromagnetic ordering with periodicity of two lattice sites, ${{\bf q}=(\pi,\pi/2)}$, in the region between the ${(\pi,\pi)}$ and the ${(\pi,0)}$ phases.

In Ref. Lenz et al. (2016), the VCA with bath sites was used at zero temperature and the CDMFT at finite temperatures to study the anisotropic model restricted to the paramagnetic case. These calculations yield a zero-temperature phase diagram containing two phases: a paramagnetic metallic phase at small $U$ and a paramagnetic insulating phase at large $U$ for all ${t_{\perp}>0.0}$. In particular, the transition between the two phases was found to occur at finite values of the Coulomb repulsion for all $t_{\perp}$. Interestingly, the transition was found to be continuous for ${t_{\perp}\lesssim 0.2}$ and discontinuous for ${t_{\perp}\gtrsim 0.2}$.

Dimensional crossover in weakly coupled chains at finite temperature was studied by Raczkowski et al. Raczkowski et al. (2015) using DQMC. For ${U=2.3}$ the authors found a transition between quasi-one-dimensional and two-dimensional antiferromagnetic ordering similar to the one we find.

In the context of the two-dimensional Hubbard model, it is important to address the limitations of our calculations in understanding the discrepancies with other calculations using different methods. The primary limitation of our calculations was that of small system width. Indeed, we focused on systems of width 4, although, as mentioned above, calculations for ${t_{\perp}=0.7}$ were also carried out for width 5. This raises the question of what features of our phase diagram are robust with regard to scaling in the lattice width. In particular, it is possible that some features are remnants of the quasi-one-dimensional nature of the lattice, i.e., of four-chain physics.

The first feature to be considered is the presence of the intermediate ISDW phase. As shown, umklapp processes are available for the width-4 cylinders that involve a wave vector that is compatible with that of the incommensurate structure of the magnetic correlations as found in the DMRG calculations for low $U$. This indicates that the ISDW phase might indeed be a feature of the four-chain system. On the other hand, width-5 cylinders for ${t_{\perp}=0.7}$ also show the signatures of both transitions, metallic phase to ISDW and ISDW to antiferromagnetic phase.

The insulating ISDW phase with incommensurate magnetic ordering found in our DMRG calculations has, at least up to now, not been found in treatments of Hamiltonian (II) in two dimensions. It is possible that this discrepancy is due to the small cluster or lattice sizes used in the other calculations, including the VCA results presented in Sec. V.2 of this work. This argument can be supported by considering other work treating the Hubbard model on another two-dimensional, anisotropic, frustrated lattice, namely the anisotropic triangular lattice. Treating ${12{\times}12}$ anisotropic triangular lattices, Wang et al. found magnetic ordering in the half filled system at different wave vectors away from ${(\pi,\pi)}$ in DQMC calculations Wang et al. (2016). In addition, a recent treatment using dynamical mean-field theory after a local spin-rotating gauge transformation was applied to the system finds a phase with incommensurate magnetic ordering for ${t^{\prime}>0.7}$ for the triangular geometry Goto et al. (2016). Incommensurate spiral magnetic ordering was also found by Tocchio et al. in a variational Monte Carlo calculation Tocchio et al. (2013).

Furthermore, mean-field Hartree-Fock and slave-boson calculations find phase diagrams for the frustrated square lattice that feature incommensurate phases Igoshev et al. (2010, 2013); Timirgazin et al. (2016). In particular, Ref. Timirgazin et al. (2016) finds a MIT with an incommensurate intermediate phase at half filling with a modulation of the ordering wave vector that is consistent with our findings (see Fig. 1 in Ref. Timirgazin et al. (2016)).

The second feature is the discontinuous nature of both transition lines, which for the four-chain case extends over the entire length of the lines. Here we note that the discontinuous nature of the transition between the AF and ISDW phases appears to be stable with regard to scaling in the system length, whereas it is uncertain how the transition between the ISDW and metallic phases scales in the infinite-length limit. However, we cannot perform a rigorous extrapolation of the jump in double occupancy at both transitions with system size and therefore cannot completely rule out either transitions becoming continuous somewhere along the transition lines.

The third feature is that both transition lines curve down, reaching the ${U=0}$ axis at finite $t_{\perp}$; see Fig. 4(a). The lower transition line, i.e., that between the metallic and the ISDW phases, goes to zero $U$ at a value of $t_{\perp}$ close to a Van Hove singularity. Note that for the two-dimensional lattice, the Van Hove singularity associated with a change from open to closed topology of the Fermi surface occurs at ${t_{\perp}\approx 0.622}$, while for the four-chain system, the related Van Hove singularity associated with a four-band to three-band transition occurs at ${t_{\perp}\approx 0.707}$. This three-band–four-band transition point could have consequences in a weak-coupling picture. In particular, it is possible that umklapp processes could lead to a Mott transition at arbitrarily weak $U$ in the four-band regime (${t_{\perp}\lesssim 0.707}$), but not in the three-band regime, leading to a Mott transition as a function of $t_{\perp}$ at weak $U$.

The potential behavior of the upper transition line (between the ISDW and antiferromagnetic phases) upon system-width scaling is unclear. Methods that are not restricted to small finite widths, in particular, the VCA Lenz et al. (2016), obtain lines for the transition from the paramagnetic or antiferromagnetic insulators to the paramagnetic metal that lie above our upper transition line and do not curve down to the ${U=0}$ axis at finite $t_{\perp}$, suggesting that the behavior of our transition line could be a finite-width effect. Since we cannot perform an extrapolation in the lattice width at that point, we avoid further speculation about the possible outcome of such extrapolation. As mentioned above, weak-coupling renormalization-group calculations could help clarify these questions in the weak-$U$ regime. In order to make statements at larger $U$, it would be desirable to be able to accurately treat systems of larger width with the numerical methods; however, significant improvement in the algorithms would be needed to do this.

VII Conclusion

We have investigated the ground-state phase diagram of the frustrated anisotropic Hubbard model for the case of four periodically coupled chains, i.e., width-4 cylinders, at half filling when increasing the interchain coupling with fixed ratio of interchain hopping and frustration, ${t^{\prime}/t_{\perp}=-0.5}$.

Using the hybrid–real-momentum-space DMRG, we have mapped out the phase diagram as a function of the interchain hopping and the interaction strength. In the region of weak on-site interaction and strong interchain coupling, we find a metallic phase, which vanishes at the point where the noninteracting Fermi surface undergoes a topological change at the Van Hove singularity. At higher interaction strengths, we find an antiferromagnetic insulator that smoothly changes from quasi-one-dimensional antiferromagnetic ordering for weak interchain coupling to two-dimensional ordering for stronger interchain coupling. The quasi-one-dimensional antiferromagnetic behavior remains present as the interaction strength becomes weaker and even extends to larger interchain coupling at sufficiently weak interaction strength. The metallic and antiferromagnetic phases are separated by an intermediate incommensurate spin-density-wave state, i.e., a state characterized by dominant incommensurate spin correlations and insulating behavior. For fixed lattice size, we find that the transitions between these phases are discontinuous, as characterized by finite jumps in the double occupancy and the von Neumann entropy and discontinuous changes in the static spin structure factor.

Our numerical calculations were carried out on a lattice with cylindrical topology and cylinders of lengths ranging from 16 to 48 lattice spacing and, primarily, of width 4. Here, the discontinuous nature of the transition between the antiferromagnetic and the incommensurate spin-density wave is stable as a function of the length, while we cannot determine whether the transition between the latter and the metallic phase becomes continuous in the infinite-length limit. We have been able to carry out a minimal check of robustness of the results on lattice width by comparing results for width-4 and width-5 lattices at one value of $t_{\perp}$, ${t_{\perp}=0.7}$ as a function of interaction strength for one length, 16; the signatures of the two transition points found on the width-4 lattices are also present in the width-5 results.

For selected values of the interchain coupling, we apply the VCA within a cluster geometry that corresponds to the cylindrical lattice geometry used in the DMRG calculations so that the VCA results can be compared to the unbiased DMRG results. For intermediate to strong interchain coupling, the VCA obtains a discontinuous phase transition between a metallic conducting and an antiferromagnetic insulating phase with energies and double occupancies that agree well with the DMRG results. However, the restricted cluster sizes in the VCA hinders the treatment of incommensurate spin correlations, and the VCA results hence cannot pick up the intermediate incommensurate spin-density-wave phase. For weaker interchain coupling, the metal-insulator transition becomes continuous in the VCA and occurs around the same value of the interaction strength at which the antiferromagnetic order with distinct ordering wave vector ${(\pi,\pi)}$ vanishes smoothly in the DMRG. Thus, we find good agreement between the VCA and DMRG results in the region where we expect it, at intermediate to strong interchain hopping and moderate to strong interaction strength. The absence of the incommensurate spin-density-wave-phase in the VCA and the discrepancies at weak interchain coupling are understood by the limitations of the VCA clusters used. One central question raised by our study is whether the intermediate incommensurate spin-density-wave phase is pertinent to the two-dimensional system or, indeed, to finite-width cylinders of larger width. On the one hand, we have found that the magnetic structure in this phase is compatible with possibly relevant umklapp scattering processes in width-4 Hubbard cylinders, indicating that the intermediate phase could be specific to four-chain physics. On the other hand, we have found some indications of the presence of an intermediate incommensurate phase in the width-5 cylinder and, as described above, some studies on related two-dimensional systems also find similar intermediate phases. Further investigation of this issue, for example, using the renormalization group and bosonization at weak coupling, or by jointly applying and comparing various numeric methods, as has recently proven to be fruitful for the doped two-dimensional Hubbard model LeBlanc et al. (2015); Zheng et al. (2017), would certainly be very interesting.

Acknowledgements.

We thank M. Raczkowski and F. F. Assaad for helpful discussions, and we acknowledge computer support by the GWDG and the GoeGrid project. This work was supported in part by the Deutsche Forschungsgemeinschaft (DFG) through Research Unit FOR 1807, projects P2 and P7, and the European Research Council (Project No. 617196).