Sign problem free quantum Monte-Carlo study on thermodynamic properties and magnetic phase transitions in orbital-active itinerant ferromagnets

The microscopic mechanism of itinerant ferromagnetism is a long-standing problem due to the lack of non-perturbative methods to handle strong magnetic fluctuations of itinerant electrons. We have non-pertubatively studied thermodynamic properties and magnetic phase transitions of a two-dimensional multi-orbital Hubbard model exhibiting ferromagnetic ground states. Quantum Monte-Carlo simulations are employed, which are proved in a wide density region free of the sign problem usually suffered by simulations for fermions. Both Hund's coupling and electron itinerancy are essential for establishing the ferromagnetic coherence. No local magnetic moments exist in the system as a priori, nevertheless, the spin channel remains incoherent showing the Curie-Weiss type spin magnetic susceptibility down to very low temperatures at which the charge channel is already coherent exhibiting a weakly temperature-dependent compressibility. For the SU(2) invariant systems, the spin susceptibility further grows exponentially as approaching zero temperature in two dimensions. In the paramagnetic phase close to the Curie temperature, the momentum space Fermi distributions exhibit strong resemblance to those in the fully polarized state. The long-range ferromagnetic ordering appears when the symmetry is reduced to the Ising class, and the Curie temperature is accurately determined. These simulations provide helpful guidance to searching for novel ferromagnetic materials in both strongly correlated $d$-orbital transition metal oxide layers and the $p$-orbital ultra-cold atom optical lattice systems.

We investigate thermodynamic properties of itinerant ferromagnetism by using the nonperturbative method of quantum Monte-Carlo simulation, which is shown free of the sign problem in a multi-orbital Hubbard model in the square lattice in a large region of fermion density. The spin magnetic susceptibility χ is local-moment-like exhibiting the Curie-Weiss law in the off-critical temperature region, while the compressibility κ typically exhibits the itinerant nature, which is finite and weakly temperature-dependent. χ further grows exponentially as approaching zero temperature for the SU(2) invariant models. The long-range ferromagnetic ordering appears when the symmetry is reduced to the Ising class, and the Curie temperature can be accurately determined.   , which has also become a recent focus of ultra-cold atom physics [24][25][26][27][28][29][30][31] . Stoner proposed the exchange interaction among electrons with parallel spins as the driving force for FM 1 . Along this direction, the local density approximation (LDA) of the density functional theory extended with spin polarization has achieved great success 32,33 . For example, the ground state magnetic moments of FM metals are calculated accurately 34 . The implementation of correlation effects in LDA is greatly improved by the methods of LDA+U 35 , LDA+dynamical mean-field theory 36,37 , and LDA+Guzwiller projection [38][39][40] .
Another approach to FM is based on model Hamiltonians, which is important to illustrate the non-perturbative aspect of FM. The Stoner criterion overlooks correlations among electrons with opposite spins 23 : electrons can delicately organize their wavefunctions to reduce repulsion and remain unpolarized even with strong interactions. The Lieb-Mattis theorem shows that the ground state of a rigorously 1D system is a spin singlet no matter how strong interaction is 5 . Therefore, exact theorems establishing FM are important to provide reference points for further investigations. Well-known examples include the Nagaoka theorem 6, [41][42][43][44][45] , which applies to the infinite U Hubbard models in two and above dimensions with doping a single hole on the half-filled background, and the "flat-band" FM in certain lattices with dispersionless band structures 13,46 . Significant progress has been made on the study of the critical behavior of the itinerant FM transitions 8,14,[16][17][18][19] .
Thermodynamic properties are among the major challenges of understanding itinerant FM 10,47-49 . The magnetic susceptibilities of FM metals typically exhibit the Curie-Weiss (CW) law of where C is the Curie constant 50 and T 0 is the Curie temperature at the mean-field level. The CW-law often applies to a wide range of temperature T 0 < T ≪ T f with T f the Fermi temperature, which implies spin incoherence in this temperature region. It is natural for local moment systems but not obvious for itinerant fermions with Fermi surfaces. In contrast, the random phase approximation yields χ(T ) ∝ 1/(T 2 − T 2 0 ), which is not observed in experiments 10,48,49 . The reason is that the paramagnetic phase close to T 0 cannot be viewed as a simple uncorrelated state but contains dynamic fluctuations of FM domains. A lot of efforts have appeared to explain the CW-law such as the self-consistent renormalization theory including spin mode coupling 9,10,47 , the direct exchange from the Coulomb integral 11,51 , and spin incoherence due to the Hund's coupling 52 .
Recently a multi-band Hubbard model is proved to be ferromagnetic in the ground state in the strong coupling limit 21 , whose Hamiltonian in the 2D square lattice is given in Eq. 2 plus Eq. 3 below. It can also be generalized to the 3D cubic lattice. This result establishes a stable itinerant FM phase over the filling region 0 < n < 2, where n is the occupation number per site. The band structure is quasi-1D consisting of orthogonal rows and columns, and electrons do not transit among different lines. When one electron in a row meets another one in a column at the crossing site, their spins are aligned by Hund's coupling. Usually Hund's rule can only polarize electrons on the same site. Remarkably, in this case it does polarize electrons in the entire system in the limit of infinite intra-orbital repulsion 20,21 . Although the electron band structure is quasi-1D, the FM correlation and ordering are genuinely 2D or 3D.
The above ground state FM phase provides a wellcontrolled starting point for studying thermodynamic properties of itinerant FM. For this purpose, we use the continuous time stochastic series expansion (SSE) quantum Monte-Carlo (QMC), which is shown free of the sign problem for the systems to be investigated below, and thus high numeric precision is feasible. T 0 is extracted based on the CW-law in the off-critical region, which is much lower than the temperature scale of charge coherence T ch . The filling dependence of T 0 (n) is calculated and the maximal T 0 reaches one tenth of the hopping integral. When entering the critical region, for the SU(2) symmetric models, χ(T ) grows exponentially. The true FM long range order is achieved by reducing the symmetry to the Ising class and the FM critical temperature T c is determined accurately.
We consider the following multi-orbital Hubbard model in the square lattice with the quasi-1D band structure 21 . This is not only of academic interest but also a good approximation to various physical systems. For example, in many transition metal oxides with the layered structure, the dispersions of d xz and d yz -orbital bands are highly anisotropic 20,53-56 , i.e., the longitudinal bonding parameter t is much larger than the transverse one t ⊥ , and the same feature also appears in the p x and p y -orbital systems 57 . For simplicity, below we use the 2D p-orbital system as an example. After neglecting the small t ⊥term, the following Hamiltonian is arrived in which particles in the p x(y) -orbital only move along the x(y)-direction, respectively; t is scaled to 1 below. The interaction part H int contains the standard multi-orbital Hubbard interaction [58][59][60][61] as where a = x, y; n a,σ = c † a,σ c a,σ and n a = n a,↑ + n a,↓ ; S a = c † a,α σ αβ c a,β is the spin operator of the a-orbital. The U and V -terms describe the intra-and inter-orbital Hubbard interactions, respectively; the J-term is Hund's coupling and J > 0 represents its FM nature; the ∆term describes the pairing hopping process between two orthogonal orbitals.
Remarkably, the above Hamiltonian is free of the QMC fermion sign problem in the limit of U → +∞ in which the ground state FM theorem applies 21 . In the coordinate representation, a convenient set of many-body bases are defined by ordering fermions according to their real space positions along one row by another and then along one column by another. The periodical and antiperiodical boundary conditions are employed for each chain if the particle number in that chain is odd and even, respectively, which is feasible because the particle number in each chain is separately conserved. Under the above bases and boundary conditions, the manybody Hamiltonian matrix of Eq. 2 plus Eq. 3 satisfies the prerequisite of the Perron-Frobenius theorem: all the off-diagonal matrix elements, which arise from the kinetic energy term and Hund's coupling, are seminegative-definite at U → +∞ (The pair hopping term vanishes in this limit.). This fact leads to a nodeless ground state wavefunction and the absence of the QMC sign problem 21 . In the thermodynamic limit, the effect of boundary conditions should not change the bulk phases. The QMC method we will use is the SSE with the directed loop update algorithm [62][63][64][65][66] , and the sign problem is also absent at finite temperatures within this algorithm. The grand canonical ensemble is employed to ensure the ergodicity of the particle number distribution in each chain. The spin magnetic susceptibility χ is represented as the equal-time correlation function The QMC results of χ −1 (T ) at V = 0 are presented in range from 0.01 to 0.1, which means that spin remains incoherent at temperatures well below t . Due to non-Gaussian fluctuations, the actual FM critical temperature T c significantly deviates from T 0 defined in Eq. 1, and can even be suppressed to zero for the SU(2) invariant systems 67 . Nevertheless, T 0 is an important energy scale which roughly equals the energy cost of flipping an individual fermion spin in the ground state.
The charge channel physics is reflected by the compressibility κ(T ) defined as whose results are presented in Fig. 1 (b) with V = 0. κ is proportional to 1/T in the high temperature incoherent regime, and saturates at low temperatures in the metallic phase. The crossover temperature scale T ch between these two regimes is typically the chemical potential at zero temperature. In our case at V = 0, due to the infinite U and the 1D band structure, T ch is roughly the Fermi temperature of spinless fermions at the same density. For most values of n presented, T ch is at the order of t , and thus κ saturates in the temperature region presented. As for the case of n = 1.8, T ch ≈ 0.1, and thus κ(T ) does not saturate in the simulated temperature region. Because of the strong FM tendency and V = 0, the inter-orbital interaction is suppressed, κ(T ) can be well fitted by that of spinless fermions. Comparing χ and κ, the spin coherence temperature T 0 is much lower than the charge coherence temperature T ch . These two distinct temperature scales are a common feature of itinerant FM, which also appears in the 1D spin incoherent Luttinger liquids 68 .
Because the particle number of each chain is separately conserved, the momentum space distribution function is essentially 1D-like. Nevertheless, each chain is not isolated but interacts with others through multiorbital interactions. Without loss of generality, we define n F (k) = σ p x,σ (k)p x,σ (k) for a horizontal x-chain. The case of n = 1 is studied below, which is equivalent to n x = 0.5 in this x-chain, and its T 0 ≈ 0.08 as shown in Fig. 2 (a) below. The simulated results of n F (k) are presented in Fig. 1 (c) with the periodical boundary condition and a discussion on the boundary condition is in Sec. II of the Supplementary Material. We define a reference wavevector as the Fermi wavevector k 0 f = π 2 of spinless fermions at the same density. At a low temperature T = 1/β = 0.1 close to T 0 , n F (k) is only slightly larger than 1 at k ≪ k 0 f and it smoothly decays to zero with a half-width approximately equal to k 0 f . n F (k) is rounded off compared to that of spinless fermions at the same temperature. The system remains unpolarized with a FM correlation length ξ at the order of 10 ∼ 20 as estimated in Sec. II of the Supplementary Material. Although n F (k) does not looks much different from that of spinless fermions, it is a consequence of strong interactions because its upper bound is 2 as k → 0. It implies that the phase space for thermal fluctuations is not restricted to a small region close to k 0 f , and thus its entropy capacity is enhanced. It is consistent with the real space picture of fluctuating FM domains as T approaches T 0 .
Next we consider the effect of a large inter-orbital repulsion V . The ground states remain fully spin polarized as proved in Ref. [21], and χ −1 (T ) still exhibits the CW law as shown in Fig. 7 in the Supplementary Material. The most prominent effect of V is at the commensurate filling of n = 1 as shown in Fig. 1 (d), in which the system is in the Mott-insulating state. As a result, κ(T ) is suppressed to nearly zero at 0 < T < 0.5, which is small compared to the charge gap at the order of V . The energy scale of orbital exchange is at the scale of J orb = 2t 2 /V , which results in the antiferro-orbital ordering, i.e., the staggered occupation of the p x and p y orbitals 69 . As n moves away from 1, electrons become itinerant again. Nevertheless, the values of κ(T ) at V = 8 are suppressed compared to those with the same values of n and T at V = 0.
The filling dependence of T 0 (n) is presented in Fig.  2 (a) for both V = 0 and V = 8. The FM coherence is built up due to the itinerancy of fermions 21 , thus T 0 approaches zero in both limits of n → 0 (the particle vacuum) and n → 2 (the hole vacuum). At V = 0, the maximal T 0 appears around n = 1 where electrons are most mobile. The density dependence of T 0 (n) at V = 0 is nearly symmetric with respect to n = 1 exhibiting an approximate particle-hole symmetry. In contrast, it is highly asymmetric at large V . In this case, T 0 is strongly suppressed at 0 < n < 1, in which both charge and spin carriers are electrons. A large V penalizes two electrons occupying the same site, thus the effectiveness of the Hund's rule is suppressed. After n passes 1, a quick increase of T 0 appears because the extra particles on the Mott background of n = 1 can move easily to build up FM coherence. T 0 reaches the maximum roughly at the middle point between n = 1 and 2. As n → 2, T 0 becomes insensitive to V . In this region, most sites are doubly oc- cupied as spin-1 moments, and holes are itinerant but do not carry spin. The hole's motion threads spin moments along its trajectory and aligns their orientations, which is not much affected by V . Next we present the results of the Curie constant C. Assuming the local moment picture, the simple molecule field method yields C per spin moment as 1 3 S(S + 1) 10 , where S is the spin magnitude. In our itinerant case, the values of S fluctuate: C = 0 for the empty site, 1 4 for the singly occupied site, and 2 3 for the doubly occupied site forming a spin-1 moment. respectively. We plot the normalized Curie constant C/n v.s. n in Fig. 2 (b). C/n approaches 1 4 as n → 0, and 1 3 as n → 2 where most sites are spin-1 moments. Generally, C/n lies between these two limits. At V = 0, as n increases, the number of onsite triplets smoothly increases and so does C/n. Nevertheless, at large V , the onsite triplet formation is strongly suppressed at 0 < n < 1, and thus C/n is stuck at 1 4 . After n passes 1, C/n starts to increase nearly linearly toward 1/3. As n → 2, V hardly affects the number of onsite triplets, and thus C/n also becomes insensitive to V as T 0 does. Now we move to the critical region for FM ordering. Actually the total spin of each chain is not separately conserved, thus FM correlations are intrinsically 2D in spite of the quasi-1D band structure. Although there is no quantum fluctuation for the FM order parameter as a conserved quantity, thermal fluctuations are so strong that long-range FM ordering cannot appear at finite temperatures for the SU(2) symmetric models 67 . In Fig. 3 (a), we present the crossover of χ −1 (T ) from the off-critical region to the critical region based on the finite size scaling in Sec II of the Supplementary Material. It exhibits clear deviation away from the CW-law as T ∼ T 0 = 0.08. In the critical region, the FM order parameter already develops a non-zero magnitude, and its directional fluctuations are described by the O(3) non-linear σ-model. The FM correlation length increases exponentially as approaching zero temperature. χ(T ) evolves to the exponential form fitted by χ = Ae b T 0 T 70,71 , and the result in Fig. 3 (a) shows b = 3.1 ± 0.3 at n = 1, V = 0, and J = 2.
In order to obtain the FM long-range order, we modify the Hund's rule term of Eq. 3 to reduce its symmetry to the Ising class: We introduce J and J ⊥ for the spin components in the xy-plane and along the z-direction, respectively, and choose J ⊥ > J . The z-component FM structure factor S ⊥ (T, L) is defined as T χ(T, L). For the case presented in Fig. 3 (b), the finite size scaling of S ⊥ (T, L)/L 2 yields the critical temperature T c ≈ 0.134. This result is also checked from the scaling in the critical region in Fig. 3 (c) and (d). S ⊥ L −2+η v.s. T is plotted with η = 1 4 from the anomalous dimension of the 2D Ising universal class. The crossings of curves yield the value of T c consistent with that of the previous scaling. Furthermore, a good data collapse is achieved by employing the scaling form S ⊥ L −2+η = f ((T − T c )L 1 ν ) with ν = 1 of the 2D Ising class. In Sec. III of Supplementary Material, the mean-field value T 0 ≈ 0.20 is extracted. The suppression of T c from T 0 is due to critical fluctuations.
In summary, we have performed the SSE QMC simulations to study thermodynamic properties of a 2D itinerant FM model over a large region of densities based on a multi-orbital Hubbard model. The simulations are free of the sign-problem, and thus non-perturbative results can be obtained with high numeric precision. The spin magnetic susceptibility exhibits the CW-law in the off-critical region, and crosses over to the exponential growth in the critical region. There is a large temperature regime T 0 < T < T ch , in which the spin channel is incoherent while the charge channel exhibits the metallic behavior. The compressibility is weakly temperature dependent and saturates to its zero temperature value. The true FM long-range transition appears when the symmetry class is reduced from SU(2) to Ising. The finite size scaling in the critical region gives rise to an accurate determination of the FM transition temperature. In order to maintain ergodicity, the directed loop algorithm is carried out in the grand canonical ensemble in which the chemical potential µ is the characteristic variable. Nevertheless, in realistic systems, it is more natural to fix the average fermion number per site n rather than to fix µ. For example, in a system with a fixed average value of n, µ(T ) changes as varying the temperature T . Therefore, in presenting the simulation results for a fixed value of n, we have carefully adjusted µ to maintain n invariant. The obtained relations of µ(T ) at various values of n with the system size L = 8 are plotted in Fig.  4. The dependence of µ on L is weak, and thus we use the same set values of µ(T ) for even larger sample sizes except the case of V = 0 and n = 0.40. In this case, the finite size effect of µ is relatively strong, and we calculate µ(T ) at L = 10 and use it for larger sample sizes. By this method, we can maintain the values of n invariant within the error of ∆n = 0.006 for all the simulations. In Fig. 5 and Fig. 6, the finite size scaling of χ(T, L) are plotted with J = 2 at V = 0 and V = 8, respectively. In both cases, curves are fitted with the scaling ansatz χ(L, T ) = χ 0 (T ) + aLe −L/b to extrapolate the spin susceptibility χ(T ) for the infinite system size. For the data sets with β < 6, we have simulated lattice sizes L × L up to L = 30.
In Fig. 7, we present the spin susceptibility χ(T ) at a large value of inter-orbital repulsion V = 8 based on the finite size scaling presented in Fig. 6. For all the values of n, the curves of χ −1 (T ) exhibit the CW-law. The mean-field Curie temperature T 0 and Curie constant C extracted at V = 8 and J = 2 are extracted and presented in Fig. 2 (a) and (b) in the main text, respectively.
In Fig. 8, we present the finite size scaling of χ(T, L) in the critical region from with β from 7 to 12.5 with parameters V = 0, J = 2 and n = 1. The extrapolated values of χ(T ) are used in Fig. 3 (a) in the main text. Let us look at the curve with β = 10, the dependence of χ(T, L) on L converges at large values of L. The starting of convergence takes places at values of L at the order from 10 to 20, and thus we estimate the FM correlation length ξ at β = 10 also in this range. Next we discuss the momentum space distribution n F (k) which is 1D-like because the particle number of each chain is separately conserved. Due to this feature, we make a special choice of boundary conditions to remove the sign problem: the periodical (anti-periodical) boundary condition for a chain if its particle number is odd (even). Our simulation uses the grand canonical assemble, and thus both configurations with even and odd particle numbers are sampled which can also be easily distinguished. For configurations with odd particle numbers, the values of k take 2pπ/L with p = 0, ±1, ... ± (L − 1), and for configurations with even particle numbers, the values of k take 2pπ/L with p = ± 1 2 , ..., ±(L − 1 2 ). Because of this reason, we simulate n F (k) by separately sampling configurations of even and odd particle numbers and present both of them in Fig. 9 (a) and (b), respectively. Results of the sample sizes with L = 30 and 50 are presented, which shows that the finite size dependence is very weak. At both sample sizes L = 30 and 50, the differences caused by using periodical or anti-periodical boundary conditions are rather small. And the result under the periodical boundary condition is presented Fig.  1c in the main text. For the Ising class Hamiltonian with the modified Hund's rule coupling with J ⊥ = 2J = 4, we extract its mean-field value of the Curie temperature T 0 following the method presented in the main text. The spin susceptibility χ(T ) is obtained after the finite size scaling for χ(T, L) in the off-critical region as shown in Fig. 10 (a). By the linear extrapolation of χ −1 (T ) in the off-critical region, we obtain T 0 from the interception of χ −1 (T ) on the temperature axis as shown in Fig. 10 (b). The scaling is performed in the region T > 0.5 below which the deviation from the CW behavior appears. The linear extrapolation of χ −1 (T ) gives rise to T 0 = 0.20 ± 0.01.