Bond-located spin density wave phase in the two-dimensional (2D) ionic Hubbard model

We investigate the quantum phase transitions in the half-filled ionic Hubbard model on a two-dimensional (2D) square lattice using the variational cluster approach (VCA). We present explicit evidence for the tendency toward a novel intermediate phase in this model. This phase is characterized by bond-located magnetization. For weak Coulomb repulsion U, the system is a band insulator, and then undergoes a transition to the intermediate phase at the first-phase boundary U=Uc1. As U is increased beyond the second transition point Uc2, there occurs a Mott insulator accompanied by a long-range antiferromagnetic (AF) order. The bond-located spin density wave competes with the antiferromagnetism while the charge-density modulation exists all the way due to the staggered potential Δ.


Introduction
Studies of the electron correlation-induced transition from a band insulator (BI) to Mott insulator (MI) have been of considerable interest in recent years. Such a transition has been realized in a variety of extended versions of the Hubbard model [1,2], which can contribute to understanding the fundamental difference between these two insulators. The Hubbard model plays an important role in understanding the effects of correlations for its ability to drive a transition to a charge-gapped MI, where the valence electrons are localized by the interaction U . MI is characterized by a long-range antiferromagnetic (AF) order with all sites singly occupied. In contrast, there is a modulation of charge density in a BI with one sub-band fully filled and the other one empty. Besides the random case [3], the periodic external potential can also lead to a BI. Consequently, abundant attention has been paid to the ionic Hubbard model (IHM) [4]- [6], an extended version of the Hubbard model with a staggered on-site potential . Several types of phases such as the BI, MI, bond-ordered insulating phase and metallic phase can be described within this simple model via tuning the strength of parameters [7]- [9].
The IHM was originally proposed for the description of the neutral-ionic transition in charge-transfer organic salts [5] and then potentially used to interpret the enhanced response of ferroelectric perovskites [6]. The Hamiltonian of IHM is given as where t is the hopping amplitude between the nearest-neighbor sites belonging to the sublattice A and B, respectively, and c † i,σ (c i,σ ) creates (annihilates) an electron with spin σ at lattice site i. Electrons experience an on-site Coulomb repulsion U and different on-site energies at different sublattices, that is ε A = + and ε B = − . The model shows two limiting cases. (i) In the noninteracting limit U = 0, the Hamiltonian can be easily diagonalized in momentum space and the system is a BI with two bands separated by a gap 2 . (ii) In the strong-coupling limit U ( , t), the Hamiltonian can be mapped onto an effective Heisenberg model where the ground state is an MI, with a superexchange coupling J eff = 4t 2 U/(U 2 − 2 ). Therefore, a transition from the BI to MI must occur at some intermediate U .
In the extensively studied one-dimensional (1D) IHM [7,8], [10]- [12], whether two transition points or one separate these two phases has given rise to much controversy. Most studies agree on a two-point scenario in which an insulating bond-ordered phase intervenes between the BI and MI. The bosonization studies [7], effective in the weak coupling, demonstrated that the system goes from the BI to a bond-order wave (BOW) phase at the 3 first critical U = U c1 associated with the collapse of the charge gap. As U is further increased, a Kosterlitz-Thouless-type transition between the BOW and the MI at U = U c2 takes place, involving vanishing of the spin gap. However, some density-matrix renormalization-group (DMRG) calculations [10] have suggested a nonvanishing charge gap at U = U c1 , unable to precisely identify the second transition. Other DMRG [11] results, a careful finite-size scaling engaged, have come to a conclusion coincident with that by the bosonization method.
Recently, the IHM in higher dimensions has attracted much interest. The nature of the intermediate phase separating the BI and MI has given rise to plenty of controversy. Phase diagrams with a finite-width metallic stripe between the MI and BI have been presented using the single-site dynamical mean-field theory (DMFT) [9,13]. It has been pointed out that these results are obtained restricted to paramagnetic solutions and the metallic phase can only exist if antiferromagnetism is destroyed by frustration. Other DMFT solutions [14], including longrange AF order that leads to a gap in the spectrum, have implied the absence of a metallic phase. Although the finite temperature determinant quantum Monte Carlo (DQMC) works [15,16] were in favor of an intermediate metallic phase, the influence of the AF correlations on the zero temperature phase diagram remains not quite clear. An intermediate insulating phase with a possible bond order has been suggested by Kancharla and Dagotto [17]. They used the cellular DMFT (CDMFT), a cluster extension of the DMFT, to determine a phase diagram having similar features to the ID phase diagram.
In this work, we use the variational cluster approach (VCA) [18] to clarify the open questions about the ground state properties of the 2D IHM. An unusual intermediate phase featured with a bond-located spin density wave (BSDW) order is explicitly shown to exist in this model for the first time. This phase has been studied in the 1D extended Hubbard model in the presence of site-off-diagonal interactions [19,20]. The transition from the BI to MI is a twostep transition with the BSDW phase intervening between them. We elucidate the properties of the BSDW phase such as the bond pattern and amplitude of the order parameter which were not clarified in the cellular DMFT study. The charge-density wave (CDW) order coexists with the AF order and BSDW order because of the externally staggered potential , whereas the latter two compete with each other. We study the symmetry-breaking phases with the help of VCA, which has been applied successfully to many systems. The magnetic and superconductive properties of the ground state in the 2D Hubbard model [21]- [24], and the spiral phase of the half-filled Hubbard model on an isotropic lattice [25], have been studied within this approach. This paper is organized into the following four sections. Section 2 gives a brief description of the VCA method. The AF order (staggered magnetization) and CDW order (staggered charge density) are demonstrated in section 3. In section 4, we present the detailed information of the intermediate bond-ordered phase. Finally, section 5 draws conclusions based on our findings.

Variational cluster approximation (VCA)
The VCA is a self-consistent cluster approximation for correlated lattice models, which can be constructed within the self-energy-functional theory (SFT) [26]. VCA can be seen as one special case of the SFT, which includes the short-range correlations within a cluster and treats the long-range correlations exceeding the cluster scale on the mean-field level. Taking advantage of an embedding approximation, one can work directly in the thermodynamic limit and describe phase transitions even though the actual numerical treatment can only be done for a system of finite size. The SFT allows us to approximate the physics of an infinite system via an exact solution of a small isolated cluster H , so prime called 'reference system'. The essence of the SFT is to express the grand potential as a functional of the self-energy , where F[ ] is the Legendre transform of the Luttinger-Ward functional [27], and G 0 = (ω − µ − t) −1 the bare Green's function. The subscript t indicates the explicit dependence of the grand potential on the hopping terms t, or generally, all one-body terms. The notation Tr A = T ω,α A αα (iω) is used, where the summations run over the diagonal variables α, Matsubara frequencies iω = i(2n + 1)π T with integer n and temperature T . For a Hamiltonian consisting of the one-body terms t and interaction part U , the self-energy functional F[ ] is independent of t [28] and satisfies −( The grand potential satisfies Consider a so-called 'reference system' H that shares the same local interaction part U as H and a modified one-body part t . The reference Hamiltonian H is obtained by decoupling H on an infinite lattice into identical clusters and switching off the inter-cluster hopping, with the same local interaction part as H and a modified intra-cluster one-particle part t . Theoretically, t can include all the single-particle terms to optimize the results, but in practice only those related to the problem under consideration should be chosen. In this paper, the fields for AF and BSDW are included. Because the main effects are sufficiently taken into consideration by a local [29] or shortranged [30] self-energy at least for the high lattice dimension, we can take the self-energy of H , = (t ), as an ansatz for H . Taking advantage of the universality of F[ ], we have where, and G denote the grand potential and Green's function of H . The reference system must be finite and can be exactly solved for any t we choose. The functional F[ ] can be eliminated through combining equations (2) The thermodynamic grand potential can be obtained by performing a search for the stationary point on a restricted space S of trial self-energies (t ) ∈ S. This space is spanned by varying the properly chosen one-particle parameters t . In the actual calculation, at zero temperature T = 0, the numerical integration over frequencies can be replaced by the sum over the single-particle excitation energies to ensure high accuracy [31], and where β = 1/T , R is the contribution from the poles of the self-energy, which can be eliminated through equation (4).
where V is the hopping between adjacent clusters, i, j are site indices within a cluster, and K is a superlattice wave vector. For a wave vector k within the original reciprocal lattice space, we have G(k, ω) = 1 L i j e −ik(i− j) G i j (k, ω), where L is the number of sites within a cluster.

Phase diagram without consideration of bond-ordered phase
We first discuss the BI and MI only, without regard to the bond-ordered phase. The phase diagram can be compared with the results by the single-site DMFT, which is not good at dealing with the bond-ordered phase. To study the MI associated with the AF order, one can use a symmetry-breaking Weiss field that permits the AF order in H as a variational parameter. The Weiss field for AF is expressed as where M is the strength of the staggered magnetic field and Q = (π, π) is the AF wave vector. The numerical calculations are performed on L = 2 × 3, 2 × 4 and 2 × 6 clusters. In figure 1, we display the staggered magnetization m AF = n ↑ − n ↓ and staggered charge density m CDW = n B − n A as a function of U for = 1.0 (the unit of energies is taken as the hopping constant t). We can see that the CDW order m CDW decreases with increasing U and jumps at the phase boundary, illustrating that the transition is of first order. The CDW order is non-vanishing everywhere, because the system is invariant with respect to translation by two lattice sites. The sudden increase in staggered magnetic magnetization at U = U c implies that the system experiences a transition from a paramagnetic BI to an AF MI. That the critical point U c slightly increases as the cluster is enlarged reveals the convergence behavior of the order parameters.  These results are in quantitative agreement with the DMFT result [14], which includes a longrange AF order.

Bond-located spin-density wave (BSDW) phase
The main finding of this paper is the existence and property of the intermediate phase. The phase we focus on is the bond-ordered phase. We have tried both the bond-located charge-density wave (BCDW), i.e. normal BOW, and spin density wave and a multitude of possible lattice distortion patterns. The BSDW symmetry-breaking pattern displayed in figure 2 is favored from an energy point of view. Even though the BCDW, rather than the BSDW, proves to separate the BI and MI in the 1D IHM, the phase diagram can be quite different in the 2D IHM. As shown in figure 3, the BCDW favors formation of a spin-singlet connected by a strong bond. The spin density and charge density are both homogeneous in this phase. On the other hand, the BSDW (here we only discuss the z-component of the spin) describes a staggered magnetization situated on bonds between adjacent sites. Both of the BSDW and MI are distinguished by gapped charge excitations and gapless spin excitations in the 1D IHM. Since the true long-range spin order is prohibited by quantum fluctuations in one dimension [32] and Ising criticality is absent in the spin sector [33], the BSDW is just the MI in the 1D IHM. The remarkable difference between the 1D and 2D systems lies in the absence and existence of long-range spin order. In the 2D IHM, the MI possesses a true long-range AF order, which competes with the BCDW and possibly suppresses it. Meanwhile, the BSDW can exhibit distinctive spin order and it is no longer just the MI. Therefore, we can expect a competition between the BSDW and AF, which is analyzed in the following. Now we take not only the AF variational parameter M but also the BSDW variational parameter ε into consideration. The corresponding Weiss field for BSDW is where δ i, j is +1, −1 for the solid and dashed lines shown in figure 2, which correspond to the strong (weak) and weak (strong) nearest-neighbor bonds for spin up (down), respectively. The symbol ε is the intensity of bond order. We find the MI and BSDW do not coexist. If ε = 0, the ground state energy is not lowered by allowing AF. The AF solution is obtained by setting ε = 0 and vice versa. The grand potential (ε) − (0) as a function of variational parameter ε for 2 × 4 and 2 × 6 clusters is displayed in figure 4. The stationary point of the grand potential is located at ε = 0 in the BI, moving to a finite ε when the BI-BSDW transition takes place. The optimal value of ε slightly decreases with increasing cluster size, depending weakly on the cluster size.
To determine which phase is preferred, we make a comparison of the ground state energy per  site E 0 = + µ N e , N e denoted as the electron density. The BSDW is the dominant instability for U c1 < U < U c2 , because its energy gain is more than the MI over this parameter range, and the MI is favorable for U > U c2 (see figure 5).
The BSDW is featured with the order parameter (along the x-direction), shown in figures 6 and 7, m BSDW = 2 k,σ σ sin(k x )c † k+π,σ c k,σ . The phase diagrams obtained on the 2 × 4 cluster are in overall agreement with the 2 × 6 cluster. Clearly, the width of BSDW shows a trend of shrinking with increasing . Whether there is a critical value c beyond which the BSDW vanishes is hard to determine within the VCA. For comparison, we also calculate m BSDW for = 0, namely the pure Hubbard model (not shown here). Most studies are in favor of an antiferromagnetically ordered ground state in the 2D Hubbard model. In our study, for = 0, the AF order parameter m =0 AF is nonzero for U > 0. However, the BSDW solution has a lower energy than the MI solution for very weak strength of U . It is hard to investigate AF and the competition between AF and other instability for very weak U in any real-space cluster approach, since there is no well-defined reciprocal space. As a result, the nesting, which leads Metallic phase is not involved in this paper. VCA is not as effective as the DMFT or CDMFT to deal with a metal due to the lack of bath which is used as parameter order. Although the transition from an insulator to a metal was determined by observing the density of states at the Fermi level ρ 0 [35], finite ρ 0 corresponding to a metallic phase. Such a definition of metal is much broader. The Landau picture requires that the spectral function at the Fermi level should be of delta-function shape, and therefore only nonzero ρ 0 is not sufficient to define a perfect metal. We have calculated the charge gap and Drude weight using the exact diagonalization with twisted boundary condition, which will be discussed elsewhere. The charge gap only closes at a single point, which also indicates that the intermediate phase is insulating.
Different phases are characterized not only by order parameters but also by distinct behaviors of the spectral function. Taking advantage of A(k, ω) = −2 lim η→0 † Im G(k, ω + iη), we evaluate the spectral function from (0, π) to (π, 0), as shown in figure 8. At weak U , the uncorrelated BI is characterized by nondispersive peaks, having a momentum-independent gap. In the BSDW phase, the minimum of the gap is located at k = (0, π) or (π, 0) instead of k = (π/2, π/2), as a signature of the spontaneously dimerized phase. For strong interaction, the ground state is a correlated MI and the spectral function does not behave the way it behaves in the BSDW. The low-energy peaks show inward dispersion and the gap reaches its minimum value at k = (π/2, π/2).

Conclusion
In conclusion, we have conducted research on the phase transitions in the 2D IHM using the VCA based on the exact diagonalization of a finite cluster. We propose that the BI-MI transition  A(k, ω) along the (0, π) to (π, 0) direction for = 1.0 on a 2 × 2 cluster. (a) BI with nondispersive peaks, (b) BSDW having a minimum gap located at k = (0, π) or (π, 0), whereas (c) MI having a minimum gap at k = (π/2, π/2).
is not direct and there exists a BSDW phase between them. This occurs beyond the critical value U c1 > 2 and persists until the second phase boundary U c2 is reached where the MI has more energy gain and becomes the leading instability. The BSDW order parameter for = 1.0 is remarkable only for a finite parameter range and about four orders of magnitude larger than that for = 0. Our results are in qualitative agreement with those derived from the CDMFT [17], and do not contradict the solutions [14] obtained by the single-site DMFT, which is unable to deal with the bond order requiring multiple sites. The appearance of an intervening metallic phase in some single-site DMFT studies [9,13,15,16] is possibly caused by the absence of spatial correlation or the restriction to the paramagnetic sector even when the actual ground state is AF. The metallic phase would survive if there were some frustration, such as the hopping terms exceeding the nearest-neighbor distance in range [36]. Spatial correlations and long-range order play an important role in the metal-insulator transition. The questions about the precise nature of the intermediate phase and the causes for the significant differences among the results obtained by different methods remain open.