Sine-square deformed mean-field theory

We develop a theory that accurately evaluates quantum phases with any large-scale emergent structures including incommensurate density waves or topological textures without {\it a priori} knowing their periodicity. We spatially deform a real-space mean-field Hamiltonian on a finite-size cluster using a sine-squared envelope function with zero energy at system edges. The wave functions become insensitive to the misfit of the lattice and ordering periods. We successfully extract the ordering wave vectors by our deformed Fourier transformation, updating the previous results for hole-doped and spin-orbit coupled Mott insulators. The method further enables the evaluation of a charge gap beyond the mean-field level.

Introduction. How to detect and describe a variety of ordered quantum phases having a long period or a large-scale structure in crystalline solids is one of the major theoretical challenges. Over the last decades, the condensed matter field has focused on these emergent phases with spatially extended periods that exhibit extraordinary properties due to topology. Examples include magnetic skyrmions with swirling spin structures in two dimensions [1][2][3][4] and a chiral soliton lattice in quasi one-dimension (1D) [5][6][7][8][9]. These incommensurate spin textures originate from spin-orbit coupling (SOC) or competing local interactions due to geometric frustration, and serve as the potential platform of next-generation memory devices and topological transport phenomena [10,11]. More generally, the electron/hole doping often induces charge-density-wave (CDW) and spin-densitywave (SDW) orders, which are considered as the key to understanding the high-temperature superconductivity [12,13] and colossal magnetoresistance [14,15]. Moiré graphene has a long-wavelength electronic state engineered by a misfit of two layers, and there is an urgent demand to understand the origin of its unconventional superconductivity [16].
The theoretical difficulty in dealing with incommensurate orders in quantum many-body systems lies in that they have exceptionally large unit cells whose periods depend sensitively on the parameters of the Hamiltonian. In methods that rely on the reciprocal lattice representation such as a mean-field calculation, the candidate period needs to be known apriori. Cluster-based methods such as the cluster dynamical mean-field theory [17][18][19][20] or the cluster perturbation theory [21,22] have crucial problems, since the orders with periods that match the cluster size are energetically favored by accident. In calculating quantum many-body states on finite clusters, such as exact diagonalization, density matrix renormalization group (DMRG) [23][24][25], and tensor network methods [26][27][28], the results often depend sensitively on the boundary conditions [29], size and shape of the clusters [30,31]. This becomes serious in dealing with incommensurate orders since the mismatch of the period between an order and size of the clusters ex-cludes the true orders from the lowest-energy-state candidates. This can also be a problem for variational Monte Carlo methods [32][33][34][35][36] and quantum Monte Carlo (QMC) methods [37,38] which can deal with large system size, when several orderings compete with each other. Notice that such a problem does not occur for classical systems, since the thermal equilibrium state does not require a coherence throughout the system unlike the wave functions in quantum phases. In fact, the open or free-boundary conditions can offer a proper evaluation of large-scale skyrmions in classical Monte Carlo simulations [39].
In this Letter, we develop a sine-square deformed mean-field theory (SSDMF), a protocol that can accurately evaluate the lowest-energy quantum state having an incommensurate order without knowing it apriori. The sine-square deformation (SSD) modifies a Hamiltonian H by an envelope function f SSD (r) as H SSD = drf SSD (r)H(r), where f SSD (r) smoothly scales down from unity at the system center (r = 0) towards zero at the edges [40][41][42]. By now it is widely used for quantum many-body ground states, dynamics, finite temperature, and for field theories [43][44][45][46][47][48][49][50][51][52][53]. Deforming the whole terms of the Hamiltonian is known to keep the nature of the Hamiltonian unchanged [54][55][56][57], since it works as a realspace renormalization [31]. At the same time, it reduces the boundary-effect and finite-size effects since the edges have zero-energies [29][30][31]58]. When applying the SSD using the DMRG, even commensurate-incommensurate quantum phase transitions are very accurately evaluated [30]. Unfortunately, their wave vectors were inaccessible for two reasons; the many-body wave functions do not provide meaningful spatial profile for the one-body quantities as well as two-point correlations with the SSD. Even when the spatial structures were available, the standard Fourier transformation is not applied to SSD systems which break the translational symmetry. We resolve this issue by developing a dubbed SSDMF and by introducing a deformed Fourier transformation. The ordering wave vector is extracted accurately even when the system size is much smaller compared to the period of the order, and the result does not change in extrapolating it to infinite system size.
SSDMF. Let us consider a generalized Hubbard Hamiltonian,Ĥ =Ĥ 0 +Ĥ U witĥ where N is the number of sites (N = L×L for the square lattice), t σ i,j = (t σ j,i ) * is the directed transfer integral from site j to i with spin σ, µ is the chemical potential,ĉ † i,σ /ĉ i,σ is the creation/annihilation operator of a fermion with particle densityn i,σ =ĉ † i,σĉ i,σ andn i = σn i,σ . We take the origin of the positional vector r i at the center of the system. The envelop function for a periodic boundary condition (PBC) is f (r) = f PBC (r) = 1 and for SSD is [31], with a radius R = R 0 + 1/2 where R 0 is the distance of the farthest site from the origin. This function smoothly modifies the energy scale of the Hamiltonian from f SSD (r) ∼ 1 at the center to ∼ 0 at the edges.
Originally the eigenstates are characterized by their wavenumber k. However, the deformation-induced term ∝ (1/2 − f SSD (r)) introduces a moderate scattering between these k-indexed states and generates a basis of localized wave-packets [30,31]. Each wave packet has a major contribution to the eigenstate ofĤ with the eigenenergy that matches the energy scale of its location. Since wave packets are localized, the finite size effect is suppressed and the excitation energy scales with system size L as ∝ 1/L 2 [31]. The wave packets at the edges have zero energy and serve as a reservoir, while a set of states that contribute to the center reproduce intrinsic ground state properties of the original Hamiltonian. The number of particles is automatically adjusted between them for a given value of µ/t σ i,j and U/t σ i,j . Our mean-field Hamiltonian for the SSD is given bŷ Here, n i and Ŝ i are the order parameters representing the particle and spin densities, respectively, which are determined self-consistently. The spin operator is defined aŝ is the Pauli matrix. The mean fields serve as internal SSD fields and are self-consistently determined by minimizing the energy. Hole-doped 1D Hubbard model. We first consider the hole-doped 1D Hubbard model with the PBC and SSD ( Fig. 1(a)) by setting t σ i,j = −1 uniform and only between nearest-neighbor sites. The hole doping shifts the nesting vector Q = π as shown in Fig. 1(b), and yields an incommensurate SDW of wavenumber Q and CDW of 2Q at the mean-field level [59]. The coexistence is allowed since the latter is the higher harmonics of the former [59,60]. Although the Bethe-ansatz solutions indicate the absence of any long-range orders [61,62], the mean-field approximation can still capture the dominant wavenumber that should appear in its correlation functions. It also provides a good description for actual materials where a finite inter-chain coupling is relevant. We will show how the SSDMF eliminates the finite-size effects arising from the long-period CDW/SDW order.
We consider U = 2.0 and µ = 0.5 which is less than half-filling, and set the initial value of the mean fields to those with a U(1) spin-rotational symmetry about the z-axis, which can describe the CDW/SDW order. Figures 1(c) and 1(d) show the spatial distribution of particle and spin densities for N = 30 and N = 60, respectively. The PBC system does not exhibit any spatial structure for N = 30, while the SSD system shows a clear modulation of the particle density by optimizing their particle and spin densities using the edges as a reservoir. For N ≥ 60, both the SSD and PBC show features of CDW and SDW orderings. To quantitatively evaluate the results, we extract the particle density by the following formula [40]: which is reduced to the standard average particle density for the PBC. Indeed, the PBC still has a translationally symmetric wave function for periods larger than the lattice spacings in the mean-field solution. However, for the SSD such a symmetry largely deteriorates at the edges, and simply averaging the whole particle density gives a physically meaningless result. Equation (5) extracts the contribution near the center by smoothly suppressing the contribution closer to the edges. As shown in the inset of Fig. 1(d), the extracted n f for the SSD is almost Nindependent; we find the value at the thermodynamic limit n f 0.88 already at N = 20. Whereas the ones obtained by the PBC show a large oscillation even at around N ∼ 100.
Next, we show that the ordering wavenumbers can be accurately evaluated. We subtract the uniform component from the particle density as δn i =n i − n f . We first demonstrate that the standard Fourier transformation,ô q = (1/N ) N j=1ô j e −iqxj for local operatorô i and wavenumber q = 2nπ/N (n ∈ Z), does not give meaningful results. In Figs. 2(a) and 2(b) these Fourier components, | δn q | and | Ŝ z q |, are shown as the function of q. Those for the SSD have the significantly large contribution from the edge, and although we find spikes in the SSD data, their positions are not equal to the PBC ones nor to those at the thermodynamic limit, making it hard to understand their physical implication. To resolve this issue, we define a deformed Fourier transformation aŝ and we take q as a continuous variable. Equation (6) is regarded as a generalization of the Fourier transformation applied to systems with spatial variation, and has several advantages over the standard transformation. First, the average particle density corresponds to the q = 0 component as n (deform) q=0 = n f , and is reduced to the standard Fourier transformation when we use f PBC (x j ). Second, the finite size effect that usually appears due to the discretization of q is formally smeared out. Figures 2(c) and 2(d) show | δn (deform) q | and | (Ŝ z q ) (deform) | as the function of continuous q. For the PBC, although the peak position does not change from Figs. 2(a) and 2(b), some additional peaks appear. For the SSD, our transformation completely suppresses the background contribution and the clear peak structure appears, whose position is close to that of the thermodynamic limit. The insets in Figs. 2(c) and 2(d) show the size dependence of the wave number extracted from the peak. Again the ones from the SSD have no significant size dependence, and the | and | (Ŝ z q ) (deform) | obtained by Eq. (6). Although q is taken as continuous in equivalently adopting Eq.(6) to PBC and SSD, for PBC it is meaningful for q = 2π/N ×(integer) which are peak positions. Insets: 1/N -dependence of the ordering wave number Q extracted from the peak position in the main panels.
peak position already reaches the thermodynamic limit at N ≥ 60. Therefore, the SSDMF can efficiently describe incommensurate orders even when the system size is small, and more importantly, there is no artifact from the mismatch of the period between the order and lattice; by minimizing the total energy of the system the excess/deficient particles and spins are automatically removed from the system center using the almost zeroenergy edge sites as a reservoir.
Square-lattice Hubbard model with SOC. We consider a square-lattice Hubbard model at half-filling with SOC which manifests in a spin-dependent transfer integral λ that rotates the spin about the z-axis, as shown in Fig. 3(a); the nearest-neighbor transfer integrals are t ↑ i,j = −t − iλ and t ↓ i,j = −t + iλ when propagating in the +x or +y direction. This form is realized by the combination of Rashba and Dresselhaus SOC with equal amplitudes, and is known to host persistent spin helix states [63,64]. It is convenient to rewrite the first term in Eq. (1) using an SU(2) gauge field as whereĉ i = (ĉ i,↑ ,ĉ i,↓ ) T , t eff = √ t 2 + λ 2 , and θ = nesting vectors Q θ = (π − θ, π − θ) and Qπ = (π, π), and that of the doubly degenerate bands (black). Finite SOC λ = t eff sin(θ/2) induces the incommensurate spiral order characterized by | Ŝ (deform) q | as the function of q = (qx, qy). The results for U = 2.0, µ = 1.0 are shown along the symmetric line in the Brillouin zone (inset of (c)), for (c) L = 8, θ = π/ √ 10 and (d) L = 28, θ = π/30. The purple dashed line denotes the ordering wave number at the thermodynamic limit. The bottom panels show the size dependence of the peak position Q = (Q, Q). The (π, π) peak for PBC is a size effect that approaches the true Q at N → ∞.
The Fermi surfaces of the noninteracting bands are shown in Fig. 3(b). There are two vectors, Q θ = (π − θ, π − θ) and Q π = (π, π), that perfectly nest the Fermi-surfaces and drive the system to an insulator at infinitesimal U . This Hamiltonian provides a prototype example that the numerical calculations mislead the conclusions about the incommensurate orderings: a recent mean-field study [65] reports the existence of CDW and incommensurate SDW orders. However, such possibility is excluded by the following analytical argument.
Let us perform a local unitary transformation using the operatorÛ = N i=1 exp(−iθ(x i + y i )Ŝ z i ) that rotates the spin quantization axis. SOC in Eq. (7) disappears asÛĉ † i e i(θ/2)σ zĉ jÛ † = σĉ † i,σĉ j,σ , and the Hamiltonian is transformed to the SU(2) symmetric Hubbard model, which shows a standard antiferromagnetic ordering. In addition, there is a rigorous proof that the CDW order is absent [66]. By transforming back the antiferromagnetic ordering of the rotating frame to the original frame, we obtain a spiral ordering. Therefore, the reported CDW/SDW orderings [65] are determinably the artifact of the standard mean-field framework.
Our SSDMF correctly captures the spiral-ordered ground state and its periodicity. We take U = 2.0, µ = 1.0 which is half-filling, and an SOC with θ = π/ √ 10 and θ = π/30. We prepare the initial mean fields that allow the spiral order in the xy-plane, which is characterized by the wave vector Q θ . Figures 3(c) and 3(d) show | Ŝ (deform) q | and the 1/L-dependence of peak positions Q = (Q, Q). Although Q/π takes an irrational number for θ = π/ √ 10, the SSDMF gives that exact value for the system as small as L = 4. When θ = π/30, the expected spiral order has a long periodicity of size 60 × 60; the method even succeeds in accurately extracting this period, Q = (0.967π, 0.967π) (unchanged up to N → ∞), which is much longer than the size of the cluster L ∼ 24. We checked several different initial states, finding that they converge to the spirals.
Charge gap. Finally, we show that the SSDMF yields an accurate charge gap of the bulk Hubbard model beyond the mean-field level. Since SSD systems have a gapless spectrum by construction due to zero-energy edge modes, the standard formula, ∆ c = (E(N + 2) + E(N − 2)−2E(N ))/4, does not give the appropriate charge gap, where E(N f ) is the N f -particle ground-state energy. For the SSD, the intrinsic filling factor observed as a parti-cle density at the system center is determined not by N f but by µ, which is numerically proven in the 1D Hubbard model using DMRG by comparing with the Bethe-ansatz solution [30]. Varying N f does not change the center particle density if µ is unchanged. For µ in the gapless region, if N f is larger/smaller than the value expected, the excess particles/holes concentrate on edges that are regarded as baths. When µ is inside the gap, the center particle density remains constant even when we change µ or N f . At the point where µ ≥ µ c starts to go outside the gap, the particles start to be doped at around the center site (see Fig. 4(a)).
In our SSDMF by varying µ for an appropriate choice of N f (not strict), and measuring the critical value µ c , the charge gap is evaluated as ∆ c = µ c − U/2, where U/2 is the chemical potential at half-filling. The initial mean fields should be chosen carefully such that the final form of mean fields is symmetric about the system center. Figures 4(b) and 4(c) show the charge gap for the 1D chain and square lattice. Here the SOC is turned off while turning it on does not change the result as proven by the local unitary transformation. The PBC ones give the standard mean-field solutions, whereas our SSDMF is close to the exact Bethe ansatz solutions [61,62] and highly accurate QMC results [67,68]. The correlation effect beyond the mean field approximation is seemingly erased by the SSD.
Conclusion. We proposed the SSDMF and the deformed Fourier transformation that accurately characterize the incommensurate orderings in an unbiased manner without knowing their spatial periods apriori. The SS-DMF also extracts the charge gap in high quality very close to the values of the Bethe ansatz and QMC results in the bulk limit. Therefore, we conclude that the method provides decisive conclusions to many topics of modern condensed matter physics related to topology, where often the contradictory results from different numerical solvers previously made the issue controversial.