{'fraction_non_alphanumeric': 0.03276924375422012, 'fraction_numerical': 0.011352126941255909, 'mean_word_length': 4.369703149784727, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ferromagnetism in the Two-Dimensional Periodic Anderson Model

Using the constrained-path Monte Carlo method, we studied the magnetic properties of the two-dimensional periodic Anderson model for electron fillings between 1/4 and 1/2. We also derived two effective low energy theories to assist in interpreting the numerical results. For 1/4 filling we found that the system can be a Mott or a charge transfer insulator, depending on the relative values of the Coulomb interaction and the charge transfer gap between the two non-interacting bands. The insulator may be a paramagnet or antiferromagnet. We concentrated on the effect of electron doping on these insulating phases. Upon doping we obtained a partially saturated ferromagnetic phase for low concentrations of conduction electrons. If the system were a charge transfer insulator, we would find that the ferromagnetism is induced by the well-known RKKY interaction. However, we found a novel correlated hopping mechanism inducing the ferromagnetism in the region where the non-doped system is a Mott insulator. Our regions of ferromagnetism spanned a much smaller doping range than suggested by recent slave boson and dynamical mean field theory calculations, but they were consistent with that obtained by density matrix renormalization group calculations of the one-dimensional periodic Anderson model.


I. INTRODUCTION
Identifying the origin of itinerant ferromagnetism in metals and specifying simple models exhibiting it are two of the most intriguing and long-standing problems in condensed matter physics. Here we report the results of lowenergy perturbation theory calculations and supporting zero temperature quantum Monte Carlo (QMC) simulations that suggest the existence and mechanisms for ferromagnetic (FM) ground states in the two-dimensional periodic Anderson model. ¿From a historical point of view, the one-band, nearestneighbor hopping Hubbard model was one of the first models proposed to describe itinerant ferromagnetism; however, the ferromagnetic phase has never been found at physical parameter values. The numerical calculations, for example, have narrowed the extent of this phase to a small region around the Nagaoka point [1], that is, the strong coupling limit for one hole doped away from half-filling. Paradoxically, in two-dimensions (2D), this model exhibits anti-ferromagnetism at half-filling and anti-ferromagnetic correlations around half-filling at weak and intermediate couplings.
Recently, Guerrero and Noack [2] listed several possible extensions of the Hubbard model that should enhance ferromagnetism: (i) the addition of frustrating hopping terms [3][4][5], (ii) the inclusion of more than one orbital per unit cell, and (iii) the addition of more general nearestneighbor interactions [3,6]. In fact, a number of frustrated models with more general interactions, such as the t − t ′ Hubbard model, and multiband models, such as the periodic Anderson model (PAM), have ferromagnetic ground states [2,6,7]. In this paper, we focus on the properties of the two-dimensional PAM.
The PAM is often used to describe the essential physics of many transition metals and rare-earth and actinide metallic compounds including the so-called heavyfermion systems [8]. The model includes a band of "light" uncorrelated electrons hybridized with a band of heavy strongly correlated electrons. Despite intense efforts to determine its phases, only a few controlled analytical approximations and numerical calculations exist for D ≥ 2. Instead, previous work often studied the single-impurity Anderson model and focused on the competition between Kondo screening and the direct RKKY coupling between the localized spins [9,10]. This competition is present when the number of conduction electrons is at least similar to the number of singly-occupied low-lying f -levels. In D ≥ 2, much less attention has been given to the region of the phase diagram where the density of electrons in the conduction band is small; however, it is for this case that ferromagnetism has been established in one dimensional systems [7,[11][12][13][14][15].
Using the density matrix renormalization (DMRG) method in the one-dimension, Noack and Guerrero [7], for example, found partially and completely saturated ferromagnetism in the PAM. They considered a parameter regime where the energy of orbital state ǫ f and the strength of the Coulomb repulsion U were adjusted so each orbital had just one electron. The position of the orbital energy was below the lower band of the noninteracting problem, and one electron per orbital cor-responds to a 1/4-filled non-interacting problem. For a sufficiently large value of U , the model exhibited a ferromagnetic ground state. Beyond an interaction-dependent value of the doping and a doping-dependent value of U , this state disappeared. The ferromagnetic phase was a peninsula in a phase diagram that was otherwise a sea of paramagnetism except at 1/4 and 1/2 filling where the ground state of the PAM was antiferromagnetic.
Ferromagnetism seems to be readily found by meanfiled approximations in any dimensions [16][17][18][19][20][21][22]. Using a slave-boson mean-field theory (SBMFT) for the symmetric PAM, Möller and Wölfe [16] found results similar to those of Noack and Guerrero. At 1/2 filling, they found a paramagnetic (PM) or antiferromagnetic (AF) phase depending on the value of the Coulomb repulsion U . By lowering the density of electrons from 1/2 filling, they also found a smooth crossover from AF to FM order via a spiral phase. Just before 1/4 filling, they got a firstorder transition from FM to AF order. More recently, the SBMFT calcualtions of Doradziński and Spalek [17,18] found wide regions of ferromagnetism in the intermediate valence regime that surprisingly extended well below 1/4 filling.
In the low temperature dynamical mean-field theory (DMFT) calculations, Tahvildar-Zadeh et al. [19] also found a region of ferromagnetism and studied its temperature dependence. At very low temperatures, their ferromagnetic region extended over a wide range of electron density and in many cases embraced the electron density of 3/8. At 3/8 filling they proposed a specific Kondoinduced mechanism for ferromagnetism that has the conduction electrons in a spin-polarized charge density-wave anti-aligned with the ferromagnetically aligned local moments on the valence orbitals. More recently Meyer and Nolting [20][21][22] appended perturbation theory to DMFT and also predicted ferromagnetism over a broad ramge of electron filling extending below 1/4 filling.
Our previous [23] and new QMC results qualitatively agree with the DMRG work; however, the phases we find quantitatively and qualitatively disagree with those derived from the mean-field approximations. Quantitatively, we find ferromagnetism in a narrower doping range than the one predicted by the DMFT and SBMFT calculations. For densities between 3/8 and 1/2, QMC predicts a paramagnetic region. whereas mean-field theory predicts a ferromagnetic states in part of that region. In fact, at a filling of 3/8 where these calculations predict ferromagnetism, we find a novel ground state of an entirely different symmetry. Instead of ferromagnetism, QMC finds a resonating spin density-wave (RSDW) state; that is, the ground state was a linear combination of two degenerate spin-density waves characterized by the (π, 0) and (0, π) wave vectors.
We remark that the quantitative differences between the DMRG and QMC calculations and the DMFT and SBMFT calculations most likely result from the expected breakdown of mean-field theory in one and two dimensions. Probably the RSDW state was not found because it was not sought. On the other hand, trying to understand the mechanism for ferromagnetism is more fundamental. It points to the long-standing difficulty of building an understanding of the PAM upon the better understood single impurity Anderson model (SIAM) or the analogous problem of building an understanding of the Kondo lattice model upon the better understood single Kondo impurity problem. Noziéres addressed this later problem and proposed a reconciliation in terms of what he calls "the exhaustion picture" [24,25].
We note that the electron densities near 1/4 filling place the work of Noack and Guerrero [7] in the exhaustion regime. When the Coulomb repulsion U associated with the double occupancy of an orbital is large (strong coupling), it is often argued that around 1/4 filling the PAM behaves like a 1/2-filled one-band Hubbard model. In two-dimensions this Hubbard model has an antiferromagnetic ground-state generated by a superexchange interaction J sx ∼ V 2 t 4 /U 5 (for the symmetric case) where V is the strength of the hybridization in the PAM. It has also been established that the Kondo lattice model and the Hubbard model for strong coupling are isomorphic with particles in one becoming holes in the other [27]. The relevance of the Hubbard model gains additional support from Noziéres's long-standing argument that the dynamics of the screening clouds are described by an effective Hubbard model.
Because of these mappings, it is seems consistent to suggest the strong-coupling physics found by Noack and Guerrero would map onto the domain of a onedimensional 1/2-filled Hubbard model. The venerable theorem of Lieb and Mattis [26], however, excludes the possibility of the Hubbard model (with nearest neighbor hopping) from showing ferromagnetism in onedimension. Their proof relied on the obvious ability to order the electrons along the chain. The two-band nature of the PAM, however, prevents this ordering by allowing processes not possible in the Hubbard model. A similar situation would occur for the two-legged Hubbard model, if it in fact shows ferromagnetism. We will argue that in the two-dimensional PAM these same processes are responsible for the ferromagnetism. They are RKKY and Nagaoka-like and are excluded in Noziéres's [24,25] picture and Tahvildar-Zadeh et al.'s [19] interpretation of it. It is important to remark that ferromagnetism in the PAM is obtained for a large region of parameters which include realistic values.
We will base our arguments on the predictions of effective Hamiltonians generated from the PAM by perturbation theory and the results of our QMC simulations. For the PAM parameters studied, these effective Hamiltonians suggest a paramagnetic or anti-ferromagnetic state at 1/4 filling, a ferromagnetic region between 1/4 and 3/8 filling, a RSDW at 3/8 filling, a paramagnetic region between 3/8 and 1/2 filling, and an antiferromagnetic state at 1/2 filling. We see all these features in the QMC simulations.
In the Section II we will define the PAM and sketch our derivation of the effective Hamiltonians. In Section III we will summarize our numerical method, noting finite size limitations. Our results will be presented in Section IV. In the Section V, the Concluding Section, we will give a detailed contrast between our picture and select other works.

II. MODELS
The PAM is described by the Hamiltonian where d † iσ and f † iσ create an electron with spin σ in d and f orbitals at site i in a square lattice, and n f iσ = f † iσ f iσ is the number operator for the f -electrons of spin σ at site i. Elsewhere we will use a similar notation to denote quantities like n d iσ = d † iσ d iσ , the number operator of d-electrons. The lattice has N sites and the hopping amplitude t between d-orbitals is only to nearest-neighbor (n. n.) sites. The hopping amplitude V hybridizes different orbitals on the same site. We used periodic boundary conditions. ¿From (1) we define H 0 , the resulting Hamiltonian when U = 0. H 0 has two dispersive bands separated by a gap For a square lattice, the energy e k = −2t(cos k x +cos k y ). This band structure for H 0 is illustrated in Fig. 1. We note if ǫ f becomes very negative (doping way below the bottom of the lower band), ∆ approaches |ǫ f |. We also note that the widths of the upper and lower band are As ǫ f becomes very negative, W − approaches zero and W + approaches W = 8t, the band width when U = V = 0. The operators which create quasi-particles in the lower and upper bands are of the form The symmetric PAM, which has the electron filling ρ = 1/2 and U = −ǫ/2, has particle-hole symmetry. This symmetry is sufficient to prevent the fermion sign problem that plagues QMC simulations. Such simulations, performed by Vekić et al. [28], suggest the existence of a charge and spin gap for small values of U with the spin gap disappearing when U is increased to some U c ∼ 2. Above U c the system exhibits long-range antiferromagnetic order.
In a previous work we presented QMC results for the asymmetric model. QMC simulations for the asymmetric model experience a sign problem which is the reason why we used the constrained-path Monte Carlo method [29].
For fixed values of ǫ f , we varied U and hole-doped away from half-filling down to a filling of 3/8. For a large enough value of U we also found anti-ferromagnetism at 1/2 filling. This state was rapidly suppressed upon doping. At 3/8 filling we saw a sharp peak appear in the spin-spin correlation at the wavenumbers k = (0, π) and (0, π). We interpreted this peak as consequence of a state resonating between two degenerate spin-density waves characterized by reciprocal wavevectors (0, π) and (0, π).
In the present work we explore the doping range from 1/4 to 3/8 filling, arguing for a region of ferromagnetism. Part of our arguments will be based on the properties of effective Hamiltonians for two different regions of parameters used in the simulations. These Hamiltonians will be derived in Section IV. We now summarize the the constrained-path Monte Carlo method.

III. NUMERICAL METHOD
The constrained path Monte Carlo (CPMC) method is extensively described and benchmarked elsewhere [29].
Here we only discuss its basic strategy and approximation. In the CPMC method, the ground-state wave function |ψ 0 is projected from a known initial wave function |ψ T by a branching random walk in an over-complete space of Slater determinants |φ . In such a space, we can write |ψ 0 = φ χ(φ)|φ . The random walk produces an ensemble of |φ , called random walkers, which represent |ψ 0 in the sense that their distribution is a Monte Carlo sampling of χ(φ), that is, a sampling of the ground-state wave function.
More specifically, starting with some trial state |ψ T , we project out the ground state by iterating where E T is some guess of the ground-state energy. Purposely ∆τ is a small parameter so for H = T + V we can write where T and V are the kinetic and potential energies.
For the study at hand, the initial state |ψ T is the direct product of two spin Slater determinants, i.e., Because the kinetic energy is a quadratic form in the creation and destruction operators for each spin, the action of its exponential on the trial state is simply to transform one direct product of Slater determinants into another. While the potential energy is not a quadratic form in the creation and destruction operators, its exponential is replaced by sum of exponentials of such forms via the discrete Hubbard-Stratonovich transformation provided U ≥ 0 and cosh ∆τ J = e −∆τ U/2 . Accordingly we re-express the iteration step as is the set of Hubbard-Stratonovich fields (one for each lattice site), N is the number of lattice sites, P ( x) = ( 1 2 ) N is the probability distribution for these fields, and B σ ( x) is an operator function of these fields formed from the product of the exponentials of the kinetic and potential energies.
The Monte Carlo method is used to perform the multidimensional integration over the Hubbard-Stratonovich fields. It does so by generating a set of random walkers initialized by replicating |ψ T many times. Each walker is then propagated independently by sampling a x from P ( x) and propagating it with B( x). After the propagation has "equilibrated," the sum over the walkers provides an estimate of the ground-state wave function |ψ 0 .
In practice we performed an importance-sampled random walk by using the transformed iterative equation In this equatioñ Thus importance sampling changes the probability distribution of the Hubbard-Stratonovich fields, biasing it towards the generation of states with large overlap with the initial state. The factor Z is the normalization constant for the new distribution. It is associated with the weight assigned to each walker and, the weight is used in a branching process to control the variance of the results. We will not discuss this process here. We used two different estimators for the expectation values of some observable O. One is the mixed estimator and the other is the back-propagated estimator where |ψ 0 is the QMC estimate of the ground state and ℓ is typically in the range of 20 to 40. For observables that commute with the Hamiltonian, the mixed estimator is a very accurate one and converges to the exact answer as |ψ 0 converges to exact ground state. For observables that do not commute with the Hamiltonian, like correlation functions, the back-propagated estimator has been found to give very accurate estimates of ground-state properties. Significant differences between the predictions of these two estimators often exist.
To completely specify the ground-state wave function for a system of interacting electrons, only determinants satisfying ψ 0 |φ σ > 0 are needed because |ψ 0 resides in either of two degenerate halves of the Slater determinant space, separated by a nodal surface N that is defined by ψ 0 |φ σ = 0. The degeneracy is a consequence of both |ψ 0 and −|ψ 0 satisfying Schrödinger's equation. The sign problem occurs because walkers can cross N as their orbitals evolve continuously in the random walk. Asymptotically they populate the two halves equally, leading to an ensemble that has zero overlap with |ψ 0 . If N were known, we would simply constrain the random walk to one half of the space and obtain an exact solution of Schrödinger's equation. In the constrained-path QMC method, without a priori knowledge of N, we use a trial wave function |ψ T and require ψ T |φ σ > 0. This is what is called the constrained-path approximation.
The quality of the calculation clearly depends on the trial wave function |ψ T . Since the constraint only involves the overall sign of its overlap with any determinant |φ , it seems reasonable to expect the results to show some insensitivity to |ψ T . Through extensive benchmarking on the Hubbard model, it has been found that simple choices of this function can give very good results [29].
Besides as a starting point and as a condition constraining a random walker, we also use |ψ T as an importance function. To reduce variance, we use ψ T |φ σ to bias the random walk into those parts of Slater determinant space that have a large overlap with the trial state. For all three uses of |ψ T , it clearly is advantageous to have |ψ T approximate |ψ 0 as closely as possible. Only in the constraining of the path does |ψ T = |ψ 0 generate an approximation.
We constructed |ψ T = σ |φ σ T from the eigenstates of the non-interacting problem. Because the z-component of the total spin angular momentum is a good quantum number, we could choose unequal numbers of up and down electrons to produce trial states and hence ground states with S z = 1 2 (N ↑ − N ↓ ). Whenever possible, we would simulate closed shells of up and down electrons, as such cases usually provided energy estimates with the least statistical error, but because we wanted to study the ground state energy as a function of S z , we frequently had to settle for just the up or down shell being closed. In some cases, the desired value of S z could not be generated from either shell being closed. Also we would select the non-interacting states so |ψ T would be translationally invariant, even if these states used did not all come from the Fermi sea. The use of unrestricted Hartree-Fock eigenstates to generate |φ σ T instead of the non-interacting eigenstates generally produced no significant improvement in the results.

A. Effective Hamiltonian: Wannier Orbital Approach
Our first effective Hamiltonian explicitly targets cases where the lower band of the non-interacting model is very flat. Such cases exist for −ǫ f W/2 > V > 0. In this regime of parameters and around 1/4 filling, single electron occupancy of the f-states can occur because of the depth of the orbital state as opposed to the double occupancy penalty of the Coulomb repulsion. We will begin by building Wannier orbital operators [30] for each band from the quasi-particle operators defined in (5) where Rewriting H 0 in the Wannier basis, we find that H 0 is simply the sum of two hopping terms corresponding to the lower and the upper bands. By construction, no hybridization exists between these two bands. The cost for this simplification is non-zero hoppings τ α ij and τ β ij between any pair of Wannier orbitals i and j in the same band. Next we rewrite the interaction term (19) in the Wannier basis This expression appears more complex than the one in the original basis; however, from it we can more conveniently derive a low energy effective Hamiltonian for electron fillings less than 1/2 filling. To do this we first require that ∆ > U , i.e. the system is a Mott insulator for 1/4 filling, so we can initially consider a H with g ii ′ ll ′ = j a ij a i ′ j a lj a l ′ j . To identify more easily the physically different contributions, we rewrite the U term to produce where ′ means that there are no repeated indices. We see that H 1 band is an extended Hubbard model with long range hoppings, a ferromagnetic exchange interaction, correlated hoppings, and a term which destroys a spinanti-aligned pair of electrons in sites i ′ and l ′ and creates an anti-aligned pair at i and l.
Again for a large range of parameters, the lower band of the PAM is quite flat. If we regard |ǫ f | as very large, we can Taylor series expand Eq. 6, substitute the result into Eq. 16, and obtain Matrix elements for i and j beyond nearest neighbors (n. n.) are smaller by higher powers of V /|ǫ f |. Thus the Wannier operator α † j is predominately f † j as the amplitudes a ij and b ij strongly decrease with the distance between i and j.
With these results we see that J e il and ω ill ′ are poportional to t 2 V 4 /|ǫ f | 6 , while g ii ′ ll ′ is proportional to t 3 V 6 /|ǫ f | 9 . This means that we can neglect the last term of H 1 band , In this one band Hamiltonian the Coulomb repulsion is no longer just on-site. Its spatial extension depends on the spatial extension of the Wannier orbitals. Because of original form of this interaction, it still only affects singlet states, but these states can now be non-local. We note that even if U/t ∼ 1 this one-band model is in the strong-coupling regime:Ū ∼ U and τ α ij ≈ t(V 2 /|ǫ f | 2 ) if i and j are nearest neighbors soŪ /τ α ij ≫ 1. In the new basis the narrow band appears narrower, the heavy fermions appear heavier, but the interaction experiences little renormalization.
The non-locality is the origin of the ferromagnetic Heisenberg term. In this term, provided U is not arbitrarily large, the ferromagnetic exchange interaction J e il ∼ U t 2 V 4 /|ǫ f | 6 , however, is smaller than the antiferromagnetic super-exchange interaction which is of order we remark is that the magnitude of the AF interaction is very small. For this reason, the lowest order terms involving the upper band are crucial to determine the magnetic phase of the system doped above 1/4 filling.
These lowest order upper band processes come from terms in the Eq. 20 with one β operator and can be written as Here the terms with l = l ′ can be neglected when the lower band is flat so we can rewrite Eq. 27 aŝ Thus the lowest order inter-band process are correlated hoppings between α i and β j orbitals and are proportional to the spin polarization (opposite to the spin of the electron which hops) surrounding the j site: the hoping occurs out of regions of ferromagnetically aligned electron spins.
It is interesting to note that the origin of this ferromagnetic alignment is not an exchange mechanism but processes involving charge-transfer. To see this more clearly we show in the Appendix that by using a Schrieffer-Wolff transformation we can reduce the multi-band term to an effective one band term |ǫ f | 2 s 2 k , for i and j second and third n. n.
with s k = cos k x + cos k y . The hopping between two α orbitals is only possible if both sites are occupied with electrons having the same polarization. Then it is clear that H 1 band induces a ferromagnetic interaction between the localized f -state which comes from the itinerancy of the carriers trough the upper band. In addition, to maximize the energy gain, the added carriers must occupy k ∼ 0 states. In this way, the magnitude of the hopping t ij (k) is maximized at the same time the denominator E + (k) −ǫ f − U is minimized (tends to ∆ − U ). The electrons can be added to k ∼ 0 states only if the background is ferromagnetic and, of course, opposite to the spin of the added electron. These charge-transfer processes enhancing ferromagnetism involve the states in the lower part of the upper band.
This inter-band process is illustrated in Fig. 2a. The representation emphasizes the collective nature of the mechanism. The virtually hopping electron has reached a k = 0 band state. It is anti-aligned with the f -orbitals which are ferromagnetically aligned among themselves. The moment of this band state partially compensates the fully saturated ferromagnetic alignment of these orbitals. This compensation picture differs from "the exhaustion picture"of Noziéres [24,25] evoked by Tahvildar-Zadeh et al. [19] The compensation is on a collective state to collective state basis and not the collective state to single moment basis argued by Noziéres. This difference highlights the difficulty building the physics of the periodic Anderson model from the physics of its single impurity version.
The process in Fig. 2b contrasts that in Fig. 2a. This process in Fig. 2b compensates on a site-to-site basis and is the one present for a paramagnetic or an antiferromagnetic alignment of the f -orbitals. The energy cost for virtual hopping is higher in than that in Fig. 2a. This leads to larger energy denominators in perturbation theory and in turn leads to a smaller lowering of the ground state energy. This process more effectively lowers conduction electron kinetic energy than mechanism (b). In the latter the conpensation is a one-to-one on-site process and is present for an anti-ferromagneitc or paramagnetic aligment of the f -orbitals. A paramagnetic alignment is shown

B. Effective Hamiltonian: Canonical Transformation
Here we will present an effective low energy model valid in a different region of parameter space: U ≫ W/2, −ǫ f ≫ V > 0. In this regime of parameters and around 1/4 filling, single electron occupancy of the fstates occurs mainly because of the double occupancy penalty of the Coulomb repulsion. As we will see, the interaction between the moments in these f -states is dominated by the RKKY interaction.
To derive this effective Hamiltonian, we will make a fourth order Schrieffer-Wolff transformation [35] as in Ref. [34]:H where the transformation operators W 1 and W 3 are of order V and V 3 . With this transformation we get a new HamiltonianH without terms of order V and V 3 . By means of another canonical transformation we eliminate the term of order V 2 . In this way we get the low energy effective Hamiltonian H spin correct through V 4 . The details of the derivation are given in the Appendix.

The final expression for the effective Hamiltonian is a Heisenberg Hamiltonian
where The RKKY contribution J (0) ij is given by where ij is associated with virtual processes where one felectron at sites i and j go through the conduction band, doubly occupy the d-sites j and i, and one of these two electrons comes back to an empty f -orbital. The expression for J The other contributions to J ij are given by where m kk ′ , n kk ′ , x kk ′ and y kk ′ are defined in the Appendix. These four contributions to J ij are about one order of magnitude smaller than J (0) ij . We expressed the magnitude of the effective spin-spin interactions J ij (33) in terms of k-space states because in the absence of hybridization they are the eigenstates of the problem. From this point of view, the ferromagnetism for low density of carriers comes from the RKKY term J (0) ij and the small volume enclosed by the Fermi surface: When the Fermi volume is small, the transferred k-vector between different points on the surface is small, and the relevant phase factors of Eq. 34 are positive for small distances between sites i and j. This positivity gives an effective ferromagnetic interaction between near neighbors, which is proportional to the density of carriers and competes with the usual anti-ferromagnetic superexchange interaction included in the other components of J ij . Because the ferromagnetic interaction increases with the concentration of carriers, we expect the appearance of a ferromagnetic phase above some critical concentration ρ c . Due to the non-interacting case (U = 0) being paramagnetic, we also expect a ferromagnetic phase above some critical value U c . To derive a simple expression for the phase boundary U c ≡ U (ρ c ), we have to compare the effective ferromagnetic interaction J (0) ij with the the antiferromagnetic super-exchange interaction, which is proportional to V 4 t 2 /U 5 for the symmetric PAM. From Eq. 34 it is clear that, for the dilute case, J (0) ij is proportional to the concentration of carriers ρ c and V 4 /tU 2 . Consequently, the phase boundary of the ferromagnetic region is given by which implies that This expression qualitatively agrees with the numerical results of Ref. [2]. This ferromagnetic phase must disappear above some upper critical concentration of carriers where the volume enclosed by the Fermi surface is no longer small. For this concentration we do not expect an ordered phase unless nesting is present. In presence of nesting, E − (k) = E − (k + Q) = E Fermi , the Fermi surface may be unstable towards the development of a spin density wave. This instability is manifested in Eq. 34 from the divergence of the Q Fourier component of J ij . For a commensurate state, Q is a high symmetry point of the Brillouin zone and equals one-half of a reciprocal lattice vector. For a square lattice two such points exist: (0, π) for ρ = 1/4 and (π, π) for ρ = 1/2. In the absence of nesting no well defined Q vector connects different points of the Fermi surface, and preferred ordering for the spins is absent. For this reason, if we were to increase ρ over the upper ferromagnetic phase boundary, we would expect a paramagnetic phase for any concentration different from 1/4 and 1/2. This picture is similar to the one already obtained in 1D [7]. The only difference between 1D and 2D is in 1D nesting exists at the Fermi level for any concentration of carriers. This nesting explains why spiral critical order at 2k F is obtained for any value of ρ in 1D [7].

C. Quantum Monte Carlo Results
All our simulations were performed for lattices of 4 × 4 and 6 × 6 unit cells. The cost of performing these simulations is approximately the same as simulating an 8 × 8 and 12 × 12 one-band Hubbard model. As mentioned previously, we used ground states derived from the noninteracting problem as |ψ T . When we simulated the PAM with ǫ f = −5, we used the non-interacting states for ǫ f = −2 because this choice consistently produced a lower estimate of the ground state energy with a smaller statistical error then we would obtain if we had used the non-interacting states for ǫ = −5.
Most of our simulations were performed for electron fillings of 1/4 through 3/8. In our previous work [23], where we studied the region from 3/8 to 1/2 and found the RSDW ground state at 3/8 filling for a 6 × 6 system. At 1/2 filling, we found the expected AF ground for U 2 in agreement with other QMC simulations. In between, we found a paramagnetic (PM) ground state. For these previous studies, we had ǫ f = −2. As part of the present study, we repeated some simulations over this 3/8 to 1/2 range with ǫ f = −5 but did not find any indications of a FM state. We note again that the predictions of DMFT and SBMFT are inconsistent with these PM and RSDW states. The effective Hamiltonians presented in this paper are expected to be inappropriate for this filling range. We will not discuss this range further.
In the present work, we also did a series of rough simulations at fillings less than 1/4 for various values of U and ǫ f . We only found PM ground states. DMFT and SBMFT find FM ground states for identical parameter choices. We note that we found a PM ground state for one electron removed from 1/4 filling, even for large values of U . As we will discuss below, when one electron was added to 1/4 filling, we found a FM state for ǫ f = −5. Clearly the properties of the PAM are asymmetric about 1/4 filling in contrast to generally accepted suggestions that it behaves as a half-filled, nearest neighbor hopping Hubbard model which displays particle-hole symmetry. Figures 4 and 3 show the main results of our simulations. They plot the energy difference between a polarized and a singlet ground state as a function of electron filling. Figure 4 is for the 4 × 4 lattice, and Fig. 4 is for the 6 × 6 lattice. In both figures, t = 1, V = 0.5, U = 4, and ǫ = −5. For the non-interacting problem, these parameters are the "flat band" case illustrated in Fig. 1b.
At 1/4 filling, the two energies are equal to within statistical error, indicating the ground state is likely PM. This is not what is expected. Our effective one-band Hamiltonian has long range hopping. So for a sufficiently large value of U , we expect an AF ground-state as was observed in the one-dimensional DMRG calculations. Calculating J RKKY for the effective Heisenberg model, we found that when i = j the magnitude of the exchange interaction is about two orders of magnitude smaller than the magnitude of our statistical error. Thus we should not expect to see this state easily. If we used an antiferromagnetic state for |ψ T , we would see a AF ground state, but the energy of this biased result and the one from a PM |ψ T were typically equal within statistical error. By doping 1/4 filling with one extra election, we found that the FM state has a lower energy. Adding one electron more produced a PM state. This behavior with the doping of one and two electrons is reminiscent of Nagaoka ferromagnetism and its instability in the U = ∞ nearestneighbor-hopping Hubbard model. While our model is not the same, the proposed mechanism for ferromag-netism is very similar: An electron lowers its kinetic energy by moving through a ferromagnetic background. In our mechanism, the electron lowers its kinetic energy by inter-band processes enabling its hopping between two occupied f -states anti-aligned with its spin. The need for inter-band processes differentiates our mechanism from Nagaoka's. We believe the disappearance of the FM with the second electron is not analogous to the known instability of the Nagaoka state but rather is a finite-size effect. It is well documented that the shell structure of a finite-sized non-interacting problem is reflected in the behavior of the filling dependence of the energy in the interacting problem [31][32][33]. At least to a first approximation, the chemical potential is constant in a shell and is discontinuous between such shells. According to our effective one-band model, if we were to add one electron to the 1/4-filled case, it would most effectively lower its kinetic energy by moving between oppositely aligned f -states by virtually hopping through a k = 0 conduction (upper) band state. The Pauli principle blocks, or frustrates, this process for a second electron of the same spin if finite-size effects make the k = (±π/2, ±π/2) states of the next shell energetically inappropriate: The second electron would have to enter the k = 0 shell oppositely aligned with the first, and the ferromagnetic background of f -states cannot accommodate the two different alignments.
The influence of finite-size and inter-band process is also seen in Fig. 4 where we show the energy ratio E/|E(S z = 0)| as the number of electrons increases from 18 to 26 electrons, from one closed conduction band shell, N e = 18, through another, N e = 24, and ending at an open conduction electron shell. At first the polarized ground state becomes much lower than the PM state and then two become approximately equal at N e = 24. This variation is highlighted in Figs. 5a-e which show the ground state energy as a function of S z for different electron fillings. For 18 and 26 electrons, the energy per site either does not vary within statistical error or increases. For 20,22, and 24 electrons, it has a clear minimum but the minimum for 24 electrons is just barely below unity. We note that the minimum in energy and the maximum in S z occurs for a half-filled conduction electron shell. It is as if a Hund's rule coupling is operative in a nonmulti-orbital situation. In Fig. 4, S max is the value of S z for which E/|E(S z = 0)| is a minimum. If this ratio had no clear minimum, then we used S max = 1 2 (N e − N s ) where N s is the number of lattice sites (the number of f -orbitals) and N e is the number of electrons.
N e = 24 corresponds to 3/8 filling. Instead of being a possibly a weak FM state, we argue that it is actually an unpolarized state RSDW state; however, while we do not clearly see this state for the 4 × 4 lattice, we clearly saw it for the 6 × 6 lattice when N e = 56 [23]. There is an important and subtle difference between the two lattices: Our performance gain in constructing |ψ T is based on using closed shell states of the non-interacting problem, and these states mix information about both the valent and conduction bands. For these states, the N e = 24 singlet for 4 × 4 is a half-filled shell, but the N e = 56 singlet for 6 × 6 is a closed shell. This difference makes it much more difficult to see the RSDW state for the 4 × 4 lattice. We emphasize that the shell in Fig. 4 is different. It is a remnant of a shell in the conduction band.
We remark that our effective Heisenberg models admits a RSDW state. In particular, within the effective Heisenberg model, the ferromagnetic nature of the RKKY exchange interaction depends on the volume enclosed by the Fermi surface being small. Increasing the number of electrons increases the volume, lengthens the k-vector for transfer across the Fermi surface, and decreases the strength of the interaction until the Fermi surface becomes unstable because of nesting. Direct calculation of the Fourier components of J ij being zero if sites i and j are on different sublattices and behaving antiferromagnetically if the sites are on the same sublattice. This is precisely the pattern observed in the spin-spin correlation function reported in our previous work, but we do not see it for the 4 × 4 lattice size. Finite sizes effects and the limited range of system sizes we can afford to simulate have prevented us from performing the scaling analysis need to establish the RSDW state as one of long-range order.
Both the ferromagnetism and conduction band shell structure is evident for the 6 × 6 case shown in Fig. 3.
Adding one electron to the 1/4 filled case (N e = 36) produces a FM state. Adding another closes the k = (0, 0) shell which frustrates the ferromagnetism. Adding more electrons to reach N e = 54 successively populates the k = (±π/3, ±π/3) and (±2π/3, ±2π/3) conduction band shells. At the middle of each shell, the energy difference is a maximum and the shell is half filled.

V. CONCLUSIONS
Our numerical results indicate that several important features of the phase diagram of the one dimensional PAM [2,7] are preserved in two dimensions. In both dimensions and for half-filling, the Coulomb interaction induces an insulating gap, and the system can have AF order (AF insulator), or it can remain in a paramagnetic state (Kondo insulator) if there is a strong enough hybridization V between both bands. In this latter case there is also a spin gap associated with the energy necessary to break a Kondo singlet. The AF order originates in the Fermi surface nesting at Q = (π, π).
When doped away from half-filling, the system in two dimensions becomes paramagnetic. In one dimension, however, critical inconmensurate correlations peaked at Q = 2k f develop. This can be understood with the effective Heisenberg theory derived in Section IV.B: The RKKY interaction J (0) ij has a divergent Fourier component at Q = 2k f because there is nesting for any concentration of electrons. This results in a Luttinger liquid with spin-spin correlations which are critical and peaked at Q = 2k f . The situation is different in two dimensions where there is nesting only for fillings 1/2, 3/8, and close to 1/4 (small Fermi surface). Our numerical evidence suggests that the system is paramagnetic between 3/8 and 1/2 filling. For 3/8 filling, where the nesting appears at Q = (0, π) and (π, 0), our CPMC results indicate the presence of a RSDW phase. In this phase, the two interpenetrating sublattices are decoupled, and there is AF order in each. This spatial ordering can be also understood by considering the effective Heisenberg theory for the PAM: It is clear that in going from 1/2 to 1/4 filling the nearest neighbor RKKY interaction changes its sign (AF close to 1/2 and FM close to 1/4). A cancellation then must occur at some intermediate concentration. We can easily see that this intermediate filling is 3/8 where there is nesting for two different wave vectors: Q = (0, π) and (π, 0). These two transferred wave vectors give canceling contributions to the nearest neighbor RKKY interactions. For instance, the nearest neigbhor in the x-direction will feel an antiferromagnetic interaction coming from Q = (π, 0) and a ferromagnetic one of the same magnitude coming from Q = (0, π). On the other hand, the next nearest neighbors have an overall phase π for both wave vectors (constructive interference) and the effective interaction is therefore AF. This is the origin of the RSDW state obtained for 3/8 filling.
By decreasing the filling below 3/8, the CPMC results indicate that the system becomes paramagnetic again down to some critical filling beyond which non-saturated ferromagnetism appears. Again this behavior is related to the nesting of the Fermi surface. Below 3/8 filling there is no nesting down to some small concentration of conduction electrons near 1/4 filling where the Fermi surface of the conduction electrons can be very well approximated by a small sphere and the wave vector 2k f connecting two different points is very close to zero. Under these conditions, the RKKY interaction J (0) ij diverges with a negative value indicating the presence of a ferromagnetic instabilty.
This mechanism is not the only one giving rise ferromagnetism in the PAM for small concentrations of carriers. In the Mott insulating regime (U < ∆) where the one band effective model H 1 band is a valid effective low energy theory, the ferromagnetism is related to the long range hopping processes involving charge virtually tranferring from the lower to the upper band. The relevance of these processes is related to the flatness of the lower band and the dispersive character of the upper one. Without considering processes connecting the two different bands (H 1 band = H (0) 1 band ), the effective one band Hamiltonian is a an extended Hubbard model (the hopping is not restricted to nearest neighbors). Due to the flatness of the lower band, the effective ratioŨ /τ α ij is very large, and the model is in a Nagaoka-like region even for small values of the bare interaction U . The relevant low energy scale which determines the magnetic ordering comes from the comparison between ρτ α ij (energy per link of the Nagaoka state) and τ α ij 2 /Ũ (energy per link of the AF state). As this difference is extremely small, any other term added to the Hamiltonian, which favors one of the two competitive phases, will be relevant. As we have explained in Section IV.A, H 1 band stabilizes the ferromagnetic phase. If the localized electrons are polarized in the same direction, the added carriers can gain energy from virtual processes transferring charge between k ∼ 0 states of the lower and the upper band.
It is also important to make some additional comments about the finite size effects in our numerical results. In an infinite system it is necessary to have a finite critical concentration of electrons in order to induce the ferromagnetic phase. We argue that This must be the case because the gain in kinetic energy of the added particles must overcompensate the loss of magnetic energy of the localized electrons (which is proportional to the system size). The numerical results show that the system becomes ferromagnetic with the addition of only one conduction electron. This is because one electron added to a 4×4 or 6×6 system corresponds to a finite concentration. Besides, as it is explicitly shown in Figs. 4 and 3, there are closed shell effects that give a non-monotonic behav-ior for the energy difference between the ferromagnetic and the paramagnetic state as function of concentration. It is clear that for small systems and closed shell conditions, the non-interacting (paramagnetic) solution will be more stable under the introduction of correlations. The small systems will generally have larger energy gaps between two shells.
Due to size effects, we must assume that we are dealing with states of long-range order, and we cannot say much about the order of the paramagnetic-ferromagnetic transitions. To estrablish long-range order and to determine the order of the transitions, it is necessary to scale the magnetization with the system size for a fixed filling. To do this properly requires larger systems than the ones considered in this paper.
The following expression forŴ 1 satisfies Eq. A4 and With thisŴ 1 , the second order term inH becomes where m kk ′ = γ k + γ k ′ , the components of the spin operator S i are and the components of s σσ ′ are the Pauli matrices divided by two. t (1) ij and t (2) ij are defined as Finally H spin is obtained by means of a second canonical transformation which eliminates H (2) fromH By expanding the exponentials of Eq. A13, we get an elimination condition that defines S: A S that satisfies this equation is ij )(n iσ − n jσ )f † jσ f iσ (A15) The second term in Eq. A16 reduces to where n d k 0 = 1 2 σ d † kσ d kσ 0 and · · · 0 means the expectation value relative to non-interacting conduction delectrons (V = 0). The first term of Eq. A17 gives the RKKY interaction.
The third term of Eq. A16, H (4) = 1 4 [Ŵ 1 , [Ŵ 1 , H (2) ]], can be expressed as APPENDIX: B In this appendix we show the derivation of H 1 band in the following way: The canonical transformation is given by the following equation:H = eT 1 He −T1 (B3) The terms of order H 1 band are eliminated under the following condition: