Determinant Quantum Monte Carlo Study of Exhaustion in the Periodic Anderson Model

TL;DR

This study uses determinant Quantum Monte Carlo to investigate how dilution affects the transition point between singlet formation and antiferromagnetic order in the Periodic Anderson model, providing insights into exhaustion physics.

Contribution

It introduces a sign-problem-free variant of the Periodic Anderson model to explore the impact of conduction electron dilution on magnetic phase transitions.

Findings

Dilution shifts the critical hybridization V for the singlet-antiferromagnetic transition.

The study demonstrates the feasibility of low-temperature QMC simulations for diluted models.

Results offer new understanding of exhaustion effects in heavy fermion systems.

Abstract

The Kondo and Periodic Anderson models describe many of the qualitative features of local moments coupled to a conduction band, and thereby the physics of materials such as the heavy fermions. In particular, when the exchange coupling $J$ or hybridization $V$ between the moments and the electrons of the metallic band is large, singlets form, quenching the magnetism. In the opposite, small $J$ or $V$, limit, the moments survive, and the conduction electrons mediate an effective interaction which can trigger long range, often antiferromagnetic, order. In the case of the Kondo model, where the moments are described by local spins, Nozi\`eres considered the possibility that the available conduction electrons within the Kondo temperature of the Fermi surface would be insufficient in number to accomplish the screening. Much effort in the literature has been devoted to the study of the temperature scales in the resulting `exhaustion' problem, and how the `coherence temperature' where a heavy Fermi liquid forms is related to the Kondo temperature. In this paper, we study a version of the Periodic Anderson model in which some of the conduction electrons are removed in a way which avoids the fermion sign problem and hence allows low temperature Quantum Monte Carlo simulations which can access both singlet formation and magnetic ordering temperature scales. We are then able to focus on a somewhat different aspect of exhaustion physics than previously considered: the effect of dilution on the critical $V$ for the singlet-antiferromagnetic transition.

Figures

Equations

$$\hat{H}=-$$

$$+U_f$$

$$P_m=\frac{4!}{m!(4-m)!}{p}^m(1-p{)}^{4-m};$$

$${\mathcal{S}}^{\mathrm{ff}}(\pi ,\pi )\equiv \frac{1}{N_{\mathrm{f}}}\sum _{\mathbf{i},\mathbf{j}}⟨S_{\mathbf{i}}^{\mathrm{f}}\cdot S_{\mathbf{j}}^{\mathrm{f}}⟩(-1{)}^{\mathbf{i}+\mathbf{j}},$$

$$C_{\mathbf{i}}^{fc}\equiv {\vec{S}}_{\mathbf{i}}^{\mathrm{f}}\cdot \text{sideset}{}^{\prime }{\sum }_{\mathbf{j}\in \mathcal{N}(\mathbf{i})}{\vec{S}}_{\mathbf{j}}^{\mathrm{c}},$$

$$\frac{1}{N}{S}^{\mathrm{ff}}(\pi ,\pi )=\frac{1}{3}{m}^2+\frac{a}{L},$$

Full text

I Introduction

A fundamental property of the description of a local magnetic moment embedded in a sea of conduction electrons provided by the Kondo model (KM) and the single impurity Anderson model (SIAM)Anderson (1961) is the screening of the moment through the formation of a Kondo singlet, a phenomenon which occurs below a characteristic Kondo temperature $T_{\rm K}$. This singlet formation is accompanied by the appearance of a narrow resonant state at the Fermi energy and a large electronic effective mass, enabling these models to provide a qualitative picture of heavy fermion physics –the enhancement of specific heat and magnetic susceptibility.Stewart (1984); Nozieres (1985); *Lee *et al. (1986) Certain features of this problem are amenable to exact analytic solution, e.g., via the Bethe ansatz.*Andrei *et al. (1983); Schlottmann (1987)

The periodic Anderson model (PAM) extends the single impurity problem to the dense limit, i.e., to a lattice of magnetic moments, raising the possibility of the emergence of magnetic ordered states. This may occur due to an indirect coupling between local moments mediated by the conduction-electron polarization (oscillations of the spin density), which is known as the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction.Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957) The Fermi wavevector $k_{\rm F}$ of the conduction electrons determines the oscillation wavelength between moments separated by distance $R$, J_{\rm RKKY}(R)\sim\big{(}k_{\rm F}\,{\rm cos}(2k_{\rm F}R)\big{)}/R^{3}. Thus the density of conduction electrons $n_{c}$, via $k_{\rm F}$, plays a crucial role in the magnetic ordering pattern.*Xavier *et al. (2004); Peters and Kawakami (2015); *Costa *et al. (2017); *Igoshev *et al. (2017); Zhong et al. (2019) Due to its importance to heavy fermion physics, the competition between singlet formation and the magnetic ordering has been investigated through many different methods, from analyticalVidhyadhiraja and Logan (2004); Hewson (1993); *Tsunetsugu *et al. (1997); Rice and Ueda (1985, 1986) to numerical.*Vekić *et al. (1995); Jarrell (1995); Rozenberg (1995); *Georges *et al. (1996); *Tahvildar-Zadeh *et al. (1997); *Vidhyadhiraja *et al. (2000); *Pruschke *et al. (2000); Capponi and Assaad (2001); *Benlagra *et al. (2011); Wu and Tremblay (2015); *Aulbach *et al. (2015a); *Hu *et al. (2017); *Schäfer et al. (2018) In particular, quantum Monte Carlo (QMC) simulationsVekić *et al. (1995); *Hu *et al. (2017) have provided evidence of the existence of a quantum phase transition from a staggered AF phase to a spin-liquid state in the two-dimensional PAM at half filling.

This competition is strongly affected by the electronic density: as $n_{c}$ is reduced, fewer conduction electrons are available to screen the local moments. Indeed, Nozières introduced the idea of “exhaustion” to describe the increased difficulty in singlet formation. Even when, naively, $n_{c}$ is large, only conduction electrons within $k_{\rm B}T$ of the Fermi surface are available for screening. Thus, Nozières suggested that the singlet formation would occur at an energy scale called the “coherence temperature” $T_{\rm coh}$, much lower than the Kondo temperature of the SIAM. In this picture, a particular functional form $T_{\rm coh}\sim N(E_{\rm F})\,T_{\rm K}^{2}/N_{\rm imp}$, where $N(E_{F})$ is the density of states at the Fermi energy and $N_{\rm imp}$ is the number of local moments, reflects the availability of only those conduction electrons within $T_{\rm K}$ of the Fermi surface.

Considerable numerical effort*Vidhyadhiraja *et al. (2000); Meyer and Nolting (2000a); Hewson (1993); *Ono et al. (1991) has gone into evaluating $T_{\rm coh}$ and its relation to $T_{\rm K}$, specifically on the validity of Nozières’ original suggestion $T_{\rm coh}\sim T_{\rm K}^{2}$. The situation is potentially complex for a number of reasons. First, the singlets formed at this scale could be rather different from those envisioned in the simpler SIAM where a single $f$ moment is screened by conduction electrons. Instead, below $T_{\rm coh}$, a much more complex tangle of spin correlations might emerge in which $f$ electrons also screen each other, i.e., singlets between $f$ electrons develop. Second, for the PAM, there are additional energy scales associated with $f$ electron charge fluctuations. In this case, it has been foundVidhyadhiraja et al. (2000) that the detailed relation $T_{\rm coh}\sim N(E_{\rm F})\,T_{\rm K}^{2}/(\,\alpha(U_{f},V)\,N_{\rm imp}\,)$ is also affected by the scales of the on-site $U_{f}$ and interband hopping $V$ energies, as opposed to a simple counting of the relative numbers of conduction and local electrons. Despite the great experimental and theoretical effortNakatsuji *et al. (2004); *Curro *et al. (2004); Yang and Pines (2008); *Yang *et al. (2008); *Shirer *et al. (2012); Wirth and Steglich (2016); *Yang *et al. (2017); *Jiang *et al. (2014); Jiang and Yang (2017); *Costa *et al. (2019), the exact relation between these two energy scales ($T_{\rm K}$ and $T_{\rm coh}$) is still an open question.

Finally, it has been suggestedNozieres (1998); Meyer and Nolting (2000a) that the physics of exhaustion might be fundamentally different depending on the strength of the coupling between the conduction electrons and local moments. For large couplings, the singlets are local and dilution of conduction electrons leaves behind well-defined local “bachelor spins” which must then find a way to form singlets. For small couplings, the screening is largely collective in the first place, even before conduction-electron dilution. This would suggest that the nature of exhaustion differs markedly at small and large $V$, in line with a more complex relation between $T_{\rm coh}$ and $T_{\rm K}$ described in Ref. Vidhyadhiraja **et al., 2000.

In this work, our main interest is to investigate how the dilution of conduction electrons affects the response of magnetic quantities in the PAM. In particular, and differently from previous work,Kaul and Vojta (2007); Watanabe and Ogata (2010); Burdin and Lacroix we are interested in determining the evolution of the quantum critical point (QCP) when the number of conduction electrons differs from the localized ones, that is, $n_{c}\neq n_{f}$. In addition, we do not rely on Pauli blocking, i.e. the restriction of conduction electron excitations to a temperature window around $E_{\rm F}$. Instead, we directly remove orbitals and their associated electrons in order to introduce depletion. In so doing, we can examine not just Nozière’s original exhaustion limit, when the ratio $p=n_{c}/n_{f}<1$, but the opposite case, $p=n_{c}/n_{f}>1$, as well. This problem is investigated with the aid of an exact numerical approach, namely the determinant quantum Monte Carlo (DQMC) method,*Blankenbecler *et al. (1981); Hirsch (1985); *White *et al. (1989); Assaad (2002); dos Santos (2003); *Gubernatis *et al. (2016) by introducing a model which allows us to control the ratio $p\equiv n_{c}/n_{f}$ without running into the sign problem.*Loh *et al. (1990); Troyer and Wiese (2005) The description of the model and the methodology are presented in the next section. Our results are shown in Sec. III, while our main conclusions are summarized in Sec. IV.

II Model and numerical method

Our work is focused on the PAM, whose Hamiltonian in real space reads

[TABLE]

Here,$c^{\dagger}_{i\sigma}\,(c^{\phantom{\dagger}}_{i\sigma})$ and $f^{\dagger}_{i\sigma}\,(f^{\phantom{\dagger}}_{i\sigma})$ are the creation (annihilation) operators of conduction and localized electrons, respectively, in the standard second quantization formalism. Similarly, $n^{c}_{{\bf i}\sigma}$ and $n^{f}_{{\bf i}\sigma}$ are site-number operators for $c$ and $f$ electrons, with $n_{\bf i}=\sum_{\sigma}(n^{f}_{{\bf i}\sigma}+n^{c}_{{\bf i}\sigma})$ being the total occupation on site ${\bf i}$. $t_{ij}$ denotes the nearest-neighbor (NN) hopping between $c$ electrons, and $V_{ij}$ represents the nonlocal hybridization between $f$-orbitals and its NN conduction sites. In the $f$-orbital dilution case, when the $f$ electron is removed from site ${\bf i}$, $V_{ij}^{f\rightarrow c}=0\ (j\in\langle i,j\rangle)$ as Fig. 1(b) shows. On the other hand, if the $c$ electron on site ${\bf i}$ gets diluted, we set $t_{ij}=0,V_{ij}^{c\rightarrow f}=0\ (j\in\langle i,j\rangle)$ as the exhaustion case plotted in Fig. 1(c). Figure 1(a) illustrates the nonlocal nature of the hybridization for the undiluted geometry; here, for simplicity, just the one-dimensional analog is shown. The local character of $f$ sites is denoted by the momentum independence of the $f$ level, $\epsilon_{f}$, and by the strong Coulomb repulsion $U_{f}$ for doubly occupied orbitals. We set $t=1$ as the energy scale and explore the “symmetric limit,” $\mu=\epsilon_{f}=0$, for which the density of the $c$ and $f$ electrons is half filled, $\langle n^{c}_{{i}\sigma}\rangle=\langle n^{f}_{{i}\sigma}\rangle=\frac{1}{2}$, a property which holds for all $t,U_{f},V,$ and temperatures $T$, due to the particle-hole symmetry (PHS) of a bipartite lattice with NN hopping terms.

At this point, we should mention that on-site (local) hybridizations are more commonly studied rather than those with NN couplings, such as the doping effect of the on-site hybridization PAM.Wei and Yang (2017) Then, it is worth emphasizing the differences between both cases beyond their dispersive character. While the former leads to a charge gap for the noninteracting limit ($U_{f}=0$) at half filling, the latter is a metal for any hybridization strength. As a consequence, the dispersive $V_{\bf k}$ is more appropriate to describe a metallic system.*Huscroft *et al. (1999); *Held *et al. (2000) Further, in the case of on-site hybridization and $U_{f}\neq 0$, there is evidence of conduction-electron localization when $f$ orbitals are removed, accompanied by an enhancement of spin-spin correlations around the unpaired noninteracting sites, breaking singlets and leading to a magnetic ground state even at large $V$.*Titvinidze *et al. (2014, 2015); *Aulbach *et al. (2015b); *Benali *et al. (2016); *Costa *et al. (2018) This effect may overestimate the magnetic response, in particular, the value of the critical hybridization for a given interaction strength, $V_{c}(U_{f})$. By contrast, such effects are strongly attenuated in the dispersive case since unpaired sites are less likely in a more connected lattice, as in the case of nonlocal hybridization. In view of this, the latter seems more relevant to study the evolution of the critical point in diluted systems.

The properties of the model are investigated using the determinant quantum Monte Carlo method,*Blankenbecler *et al. (1981); Hirsch (1985); *White *et al. (1989) which allows for an exact solution (to within statistical sampling errors) of the PAM Hamiltonian on finite-size lattices. Here we present highlights of the method, the details of which can be found in a number of reviews; see, e.g., Refs. Gubernatis **et al., 2016; Assaad, 2002; dos Santos, 2003. The underlying step is the construction of a path integral for the partition function ${\cal Z}$ by discretizing the inverse temperature $\beta=L\Delta\tau$, and breaking the full imaginary-time evolution operator $e^{-\beta\hat{H}}$ into incremental pieces $e^{-\Delta\tau\hat{H}}$. This allows for the use of the Trotter approximation,Trotter (1959); Suzuki (1976); Fye (1986) $e^{-\Delta\tau\hat{H}}\sim e^{-\Delta\tau\hat{K}}e^{-\Delta\tau\hat{\cal U}},$ which isolates the interaction term $\hat{\cal U}$, containing $U_{f}$, from the quadratic kinetic-energy pieces containing $t,\mu,$ and $V$.

The interacting $\hat{\cal U}$ term is decoupled in a quadratic form by performing a discrete Hubbard-Stratonovich (HS) transformation, with the inclusion of auxiliary fields $S({\bf i},\tau)$ in both real and imaginary coordinates, that are coupled to the spin of electrons. Therefore, the path integral consists entirely of quadratic forms and the fermionic trace can be evaluated, resulting in a trace over HS fields of a product of determinants, ${\rm det}{\cal M}_{\uparrow}(\{S({\bf i},\tau\})\,{\rm det}{\cal M}_{\downarrow}(\{S({\bf i},\tau)\})$, of matrices whose dimension is the number of spatial sites of the lattice. The trace over $S({\bf i},\tau)$ is carried out by sampling them through conventional Monte Carlo methods. Here, in addition to the usual single moves, we also perform global moves,*Scalettar *et al. (1991) which improves the ergodicity of the system.

Although the DQMC method is exact, it suffers from the infamous minus-sign problem,*Loh *et al. (1990); Troyer and Wiese (2005) which arises from the possibility of the product ${\rm det}{\cal M}_{\uparrow}(\{S({\bf i},\tau)\}\,{\rm det}{\cal M}_{\downarrow}(\{S({\bf i},\tau)\}$ being negative for certain field configurations, corresponding to a negative density matrix. The sign problem is worse at low temperatures, large lattice sizes, or strong interactions, and its dependence with the electron filling or geometries is quite nontrivial.*Mondaini *et al. (2012); *Iglovikov *et al. (2015) However, this problem is absent for systems with PHS since it implies constraints over the two determinants, leading to a positive total sign. Notice that our Hamiltonian of Eq. (II) is particle-hole symmetric at half filling, and hence the sign problem is absent throughout this work.

As mentioned earlier, the system is diluted through the direct removal of orbitals and their associated electrons. In fact, the more obvious way to reduce conducting electrons would be by simultaneously setting $\mu<0$ and $\epsilon^{f}<0$ since it preserves the number of $f$ electrons, but lowers the $c$ occupancy. However, this leads to a severe sign problem and the energy scales for singlet formation and antiferromagnetic (AF) order are no longer accessible. By contrast, our approach preserves PHS, which only requires that the hopping and hybridization should be between NN sites at half filling, hence avoiding any sign problem. Figures 1(b) and 1(c) illustrate our dilution procedure for localized ($p=n_{c}/n_{f}>1$) and conduction ($p=n_{c}/n_{f}<1$) orbitals, respectively. We should also note that depletion affects the respective neighborhoods in different ways. With the concentration of $c$ sites being $p=n_{c}/n_{f}<1$, the probability that an $f$ site is connected to $0\leq m\leq 4$ active $c$ sites is

[TABLE]

see the discussion of Figs. 8 and 9. In our measurements, for example of the structure factor below, spatial configurations with isolated $f$ electrons ($m=0$) are not included, since these sites contribute a “trivial” Curie-law free moment $\chi\sim 1/T$. Similarly, with the concentration of $f$ sites being $q=n_{f}/n_{c}<1$, the probability that a $c$ site is connected to $0\leq m\leq 4$ active $f$ sites is also given by Eq. (2), but with $p$ being replaced by $q$.

We investigate the magnetic properties by performing measurements of spin-spin correlation functions for $f$ orbitals, and their Fourier transform, the spin structure factor,

[TABLE]

where $N_{\mathrm{f}}$ is the number of $f$ sites connected to at least one $c$ site. Here, we define the fermionic spin operators as \vec{S}^{\rm f}_{\bf i}=\big{(}\,f^{\dagger}_{{\bf i}\uparrow}\,,f^{\dagger}_{{\bf i}\downarrow}\big{)}\,\,\vec{\sigma}\,\,\big{(}\,f^{\phantom{\dagger}}_{{\bf i}\uparrow}\,,f^{\phantom{\dagger}}_{{\bf i}\downarrow}\big{)}^{T}\,, with $\vec{\sigma}$ being the Pauli spin matrices; a similar expression applies to $\vec{S}^{c}_{\bf i}$. The phase factor $(-1)^{{\bf i}+{\bf j}}$ takes opposite signs on the two sublattices, corresponding to a staggered pattern.

For singlet formation, we examine a correlator function,

[TABLE]

with the sum being over sites ${\bf j}$ in the neighborhood ${\cal N}({\bf i})$ of ${\bf i}$. The prime on the sum emphasizes that in the $c$-diluted case, some $f$ orbitals have less than four $c$ neighbors; see Eq. (2). By contrast, in the case of $f$ dilution, every surviving $f$ site necessarily has four neighboring $c$ orbitals. The data reported here are obtained from lattice sizes up to $L=12$, with averaging over 20 different disorder realizations. In general the requisite number of realizations in simulations with disorder must be determined empirically and is a complex interplay between self-averaging on sufficiently large lattices, the strength of the disorder, and the location in the phase diagram. We show the results in Fig. 2. For any dilution case, the averaged $S(\pi,\pi)$ are consistent regardless of the number of realizations. These plots justify the use of 20 realizations in our work. The error bars shown in the following results reflect both statistical and disorder sampling fluctuations.

III Results and discussion

Let us first discuss the undiluted PAM [Fig. 1(a)], with which the diluted case should be compared. One should notice that due to the Mermin-Wagner theorem,Mermin and Wagner (1966) long-range order is expected to occur only at $T=0$. Therefore, as the temperature is lowered, the correlation length $\xi$ associated with spin correlations grows, but it is limited by the finite size of the system. Figures 3($a$)–3($d$) illustrate this for the AF structure factor, plotted as a function of the inverse temperature, $\beta=1/T$, for different lattice sizes $L$ and $U_{f}$, fixing $V=1$. At high temperatures (small $\beta$), $S^{\rm ff}(\pi,\pi)$ is independent of lattice size, due to the short-range character of the spin correlations. As $\beta$ increases, $\xi$ grows, ultimately reaching the linear lattice size $L$, so that the structure factor increases and stabilizes at a finite value. The growth of $S^{\rm ff}(\pi,\pi)$ with $L$ at low temperatures suggests the existence of long-range order, which should be verified through scaling arguments, as discussed below. Figures 3($e$)–3($f$) show that the behavior for the diluted cases is similar, irrespective of $n_{f}$ being larger or smaller than $n_{c}$.

There are several additional features of Fig. 3 which are worth noting: Most importantly, by comparing Figs.  3(b) and  3(e), we see that dilution of conduction orbitals enhances the AF structure factor, as it should. Second, the AF structure factor decreases as $U_{f}$ decreases, due to the magnitude of the local moments getting smaller as a result of increasing charge fluctuations. For instance, $U_{f}=2$ requires larger $\beta$ to reach the ground state than $U_{f}=5$. In fact, since the AF exchange in the Heisenberg limit is $J\sim t^{2}/U_{f}$, a large $\beta$ is also required for $U_{f}\gg W=8t$ (not shown).

We probe the existence of long-range ordering by performing a finite size scaling (FSS) analysis of the structure factor. According to spin-wave theory,Huse (1988) the AF structure factor scales with system size as

[TABLE]

with $m^{2}$ being the square of the AF order parameter, and $N=N_{f}$. Figure  4 exhibits this FSS for different values of $U_{f}$ and $V$, for the undiluted PAM. It is interesting to notice that despite the NN hybridization, our results are similar to those of the onsite case.*Vekić *et al. (1995); *Hu *et al. (2017) Our data suggest that the AF-singlet QCP is located at $V_{c}\simeq$ 1.3, 1.2, 1.0, 0.9, and 0.8 (each estimate carries a rough error bar of 0.05), for $U_{f}=$ 6, 5, 4, 3, and 2, respectively. Also, similarly to the onsite case, $V_{c}$ for the undiluted PAM is not too sensitive to the value of $U_{f}$; it only changes by approximately 50% over a range where $U_{f}$ is increased by a factor of three.

We now turn our attention to the FSS analysis for the diluted system, starting with the case $n_{f}<n_{c}$ (opposite to the exhaustion limit); see, e.g., Fig. 1(b). Here we take $q=n_{f}/n_{c}=7/8$, with the number of $c$ orbitals being $L^{2}$. Following the preceding analyses, Figs. 5($a$)–5($e$) display the data for different values of $U_{f}$ and $V$. A comparison with Fig. 4 reveals that the critical points $V_{c}(U_{f})$ are very close to those for the undiluted case. As expected, the absence of some local moments reduces the $V$ required to destroy AF order, but the effect is around 10%, i.e., of the order of $1-n_{f}/n_{c}$. This may be attributed to the longer-range character of the effective RKKY interaction between the local moments.

On the other hand, the exhaustion scenario $n_{f}/n_{c}>1$ is dramatically different. For instance, Fig.6 presents the scaling analysis for $n_{f}/n_{c}=4/3$, corresponding to one quarter of the conduction orbitals removed. It is quite evident that when the number of localized electrons is larger than the conduction ones, the AF-singlet quantum critical point is shifted to much larger values of $V$: long-range AF ordering is stabilized for $V\lesssim 2$ when $U_{f}=4$, and for $V\lesssim 3$ when $U_{f}=5$. This should be contrasted with the small changes [relative to the undiluted $V_{c}(U_{f})$] that occur when $n_{f}/n_{c}<1$ or which accompany altering $U_{f}$ at $n_{c}/n_{f}=1$.

Figures  4$(f)$, 5$(f)$, and 6$(f)$ share the feature that the AF structure factor $S^{\rm ff}(\pi,\pi)$ grows monotonically with $U_{f}$. We expect that at fixed inverse temperature $\beta$, these curves will eventually turn over and decrease, owing to the $1/U_{f}$ behavior of the exchange constant at strong coupling. In the single-band Hubbard model, maximal AF correlations occur at $U_{f}\sim 8$ for $\beta\sim 12$.*Scalettar *et al. (1991) In summary, from these preceding results, we observe that the critical hybridization for an AF-singlet transition grows dramatically when the conduction-electron count is smaller than the localized ones, while $V_{c}(U_{f})$ is only weakly changed in the opposite situation. All results for the QCP are summarized in the phase diagram of Fig. 7.

We have also investigated the singlet formation by calculating the local singlet correlator, given by Eq. (4). Figure 8 shows $C_{fc}$ as a function of $V$, that is, as the AF-singlet transition is traversed, for different values of $U_{f}$. Figures 8$(a)$–8$(c)$ display results for several values of $U_{f}$, and for the filling ratios $n_{f}/n_{c}=$1, 7/8, and 4/3, respectively. For these cases, $C_{fc}$ increases in magnitude, from small values to $|C_{fc}|\sim 0.4$, as $V$ changes from $V\sim 0.5$ to $V\sim 1.0$. The curves for the three filling ratios also all exhibit a crossing pattern: At weak hybridization $V$, $C_{fc}$ is largest in magnitude at weak coupling $U_{f}=2$. However, as $V$ increases, $C_{fc}$ becomes largest in magnitude at $U_{f}=6$; we interpret this as occurring because large $U_{f}$ yields the most well-formed moments on the $f$ sites.

We have seen that $V_{c}$ for the destruction of AF order increases dramatically for the exhaustion value, $n_{f}/n_{c}=4/3$, when $n_{f}$ exceeds $n_{c}$. It is intriguing that, in Fig. 8 there is not as great a reflection of this in the values of $V$ at which singlet correlators develop. That is, $|C_{fc}|$ grows from small values to $|C_{\rm fc}|\sim 0.4$ in the same range $0.5\lesssim V\lesssim 1.0$ for $n_{f}/n_{c}=4/3$ as for $n_{f}/n_{c}=7/8$ and 1. It appears, therefore, that the interval from $V\sim 1$ to $V_{c}\sim 3$ (for $U_{f}=5,6$) is characterized by relatively large values of the singlet correlator, even though AF order remains. Whether it corresponds to a partially screened region with coexistence between singlet and AF is a challenge to resolve conclusively with the DQMC methodology used here. While not evident in the onset of $C_{fc}$, the enhancement of $V_{c}$ by exhaustion appears to be reflected in the nonmonotonic evolution of $C_{fc}$, which is unique to the $n_{f}/n_{c}=4/3$ case; see, e.g., Fig. 8($c$). The similarity of evolution of the onset for the three filling ratios is emphasized by replotting the data of Figs. 8($a$)–8($c$) for all three ratios (at a single value $U_{f}=4$) in Fig. 8(d).

The data shown in Fig. 8 result from averaging $C_{fc}$ over all $f$ orbitals. However, as discussed in relation to Eq. (2), diluting conduction electrons leads to nonequivalent $f$ orbitals, depending on the number of active connected $c$ sites. Data for $C_{fc}$ can therefore be decomposed according to whether this number is 1, 2, 3, or 4 (the largest value for a square lattice with only near-neighbor $f$-$c$ hopping). These are, respectively, denoted by $1n,2n,3n,4n$ in Fig. 9. We show this only for the case of exhaustion, $n_{f}/n_{c}=4/3$. As it might have been intuitively expected, $C_{fc}$ is largest in magnitude for $4n$, and smallest for $1n$. It appears that $C_{fc}$ also begins to grow in magnitude at smaller $V$ for $4n$ than for $1n$.

Fig.10 shows the behavior of the double occupancy, $D_{\alpha}=\langle n^{\alpha}_{\bf i\uparrow}n^{\alpha}_{\bf i\downarrow}\rangle$, $\alpha=c$ or $f$, with $V$ for $U_{f}=3$ and $U_{f}=5$; the left panels show $D_{f}$, while and the right panels $D_{c}$. The onsite repulsion $U_{f}$ makes $D_{f}$ small, especially at small $V$ where quantum fluctuations are suppressed. By contrast, $D_{c}$ changes very little over the range of $V$ examined, taking on values close to the uncorrelated limit, $D_{c}\sim\langle n^{d}_{\bf i\uparrow}\rangle\langle n^{d}_{\bf i\downarrow}\rangle\sim 1/4$.

IV Conclusions

One of the quantitative conclusions of past QMC studies of the PAM (in its different variants) is that the position of the QCP’s, which signal the AF-singlet transition in the ground state, is rather weakly dependent on the parameters of the model, especially on the value of the on-site repulsion $U_{f}$. This was noted in Refs. Vekić **et al., 1995; Hu **et al., 2017, where the phase boundaries were found to be rather vertical in the $U_{f}-V$ plane, in contrast with mean field theory predictions of much larger $dV_{c}/dU_{f}$. Neither was $V_{c}$ found to vary much in going from 2D to 3D,*Huscroft *et al. (1999) or with changes in the momentum dependence of the $f$-$c$ hybridization.*Huscroft *et al. (1999); *Held *et al. (2000) In the former case, the noninteracting density of states is divergent at half-filling for 2D, and finite for 3D. In the latter case, the on-site and intersite forms for $V_{\bf k}$ give rise to very different band structures: insulating for $V$ independent of ${\bf k}$, and metallic for nearest neighbor hybridization. Despite these seemingly important differences, $V_{c}$ was found to be not only rather immune to changes in $U_{f}$, but also insensitive to the underlying band structure.

By contrast, one of our key results here, then, is that $V_{c}$ can vary dramatically with the ratio $n_{c}/n_{f}$. This effect is summarized in the phase diagram in Fig. 6. The physics of exhaustion is starkly evident. While dilution of local electrons barely shifts the phase boundary (singlet formation occurs slightly earlier), dilution of conduction electrons delays singlet formation to values of $V$ as much as a factor of three greater than in the balanced case, the conventional PAM. This result has the potential to lend qualitative insight into heavy fermion materials since their doping can proceed both by the elimination of moments, e.g. Ce1-xLaxCoIn5, and also by changes to the conduction electrons e.g. CeCo1-xCdxIn5. Theoretical descriptions of the latter situation have focused on the impurity-induced changes to the hybridization $V$ rather than changes to $n_{c}$.

A second key, and rather unexpected, conclusion is that in the interval $1\lesssim V\lesssim 3$ in which exhaustion induces AF order (for $U_{f}=5$) the singlet correlator is large. In the undiluted case the same $V_{c}$ pinpoints where the AF structure factor becomes small and the singlet correlator becomes large. In the presence of exhaustion, however, these two events no longer share a common $V_{c}$, and there is an extended region where both $|C_{fc}|$ is large and AF order is still present.

The large increase in $V_{c}$ found for $n_{f}/n_{c}>1$ reflects a significant increase in the stability of the AF phase, the mechanism of which can be attributed to the increased difficulty of forming singlets when the number of conduction electrons available for screening is reduced. Although we have not evaluated $T_{K}$ and $T_{\rm coh}$ here (it is possible to do so with QMC in impurity limit using, e.g., the Hirsch-Fye method,Hirsch and Fye (1986) but much harder for the lattice) it is reasonable to suppose that these temperatures will be reduced to reflect the decreased tendency towards forming singlets. We have also emphasized that the ground state structure factor $S^{\rm ff}(\pi,\pi)$ is larger for $n_{f}/n_{c}=4/3$ than for $n_{f}/n_{c}=1$, consistent with the general trend towards stronger AF with (a moderate degree of) exhaustion.

The focus of this paper has been on the competition of AF order and singlet formation (for both the cases $n_{f}/n_{c}>1$ and $n_{f}/n_{c}<1$). The physics of exhaustion appears as a large increase in AF stability. It is worth noting that the PAM has also been extensively studied as a model of ferromagnetism (FM).Meyer and Nolting (2000b); *Batista *et al. (2003) The Nagaoka theorem notwithstanding, FM appears to be very difficult to be achieved in single band systems such as those described by the Hubbard Hamiltonian. Exploration of $\mathbf{q}=0$ order would be an interesting avenue to pursue, e.g. in a more extreme limit $n_{f}/n_{c}\gg 1$ than that considered here.

Acknowledgement* — L.Z. and T.M. were supported by NSFC (Grant No. 11774033) and Beijing Natural Science Foundation (Grant No. 1192011). The numerical simulations in this work were performed on the HSCC of Beijing Normal University and Tianhe in Beijing Computational Science Research Center. N.C.C and R.R.dS. acknowledge support by the Brazilian agencies CAPES, CNPq, and FAPERJ. The work of R.T.S was supported by Grant No. DOE-DE-SC0014671.*