Phase Diagram and Dynamics of the SU(N) Symmetric Kondo Lattice Model

In heavy-fermion systems, the competition between the local Kondo physics and intersite magnetic fluctuations results in unconventional quantum critical phenomena which are frequently addressed within the Kondo lattice model (KLM). Here we study this interplay in the SU(N) symmetric generalization of the two-dimensional half-filled KLM by quantum Monte Carlo simulations with N up to 8. While the long-range antiferromagnetic (AF) order in SU(N) quantum spin systems typically gives way to spin-singlet ground states with spontaneously broken lattice symmetry, we find that the SU(N) KLM is unique in that for each finite N its ground-state phase diagram hosts only two phases – AF order and the Kondo-screened phase. The absence of any intermediate phase between the N = 2 and large-N cases establishes adiabatic correspondence between both limits and confirms that the large-N theory is a correct saddle point of the KLM fermionic path integral and a good starting point to include quantum fluctuations. In addition, we determine the evolution of the single-particle gap, quasiparticle residue of the doped hole at momentum (π, π), and spin gap across the magnetic order-disorder transition. Our results indicate that increasing N modifies the behavior of the coherence temperature: while it evolves smoothly across the magnetic transition at N = 2 it develops an abrupt jump – of up to an order of magnitude – at larger but finite N . We discuss the magnetic order-disorder transition from a quantum-field-theoretic perspective and comment on implications of our findings for the interpretation of experiments on quantum critical heavy-fermion compounds.


Fakher F. Assaad
Institut für Theoretische Physik und Astrophysik and Würzburg-Dresden Cluster of Excellence ct.qmat, Universität Würzburg, Am Hubland, D-97074 Würzburg, Germany (Dated: October 18, 2019) In heavy-fermion systems, the competition between the local Kondo physics and intersite magnetic fluctuations results in unconventional quantum critical phenomena which are frequently addressed within the Kondo lattice model (KLM). Here we study this interplay in the SU(N ) symmetric generalization of the two-dimensional half-filled KLM by quantum Monte Carlo simulations with N up to 8. While the long-range antiferromagnetic (AF) order in SU(N ) quantum spin systems typically gives way to spin-singlet ground states with spontaneously broken lattice symmetry, we find that the SU(N ) KLM is unique in that for each finite N its ground-state phase diagram hosts only two phases -AF order and the Kondo-screened phase. The absence of any intermediate phase between the N = 2 and large-N cases establishes adiabatic correspondence between both limits and confirms that the large-N theory is a correct saddle point of the KLM fermionic path integral and a good starting point to include quantum fluctuations. In addition, we determine the evolution of the single-particle gap, quasiparticle residue of the doped hole at momentum (π, π), and spin gap across the magnetic order-disorder transition. Our results indicate that increasing N modifies the behavior of the coherence temperature: while it evolves smoothly across the magnetic transition at N = 2 it develops an abrupt jump -of up to an order of magnitude -at larger but finite N . We discuss the magnetic order-disorder transition from a quantum-field-theoretic perspective and comment on implications of our findings for the interpretation of experiments on quantum critical heavy-fermion compounds.

I. INTRODUCTION
Nowadays, we are witnessing a remarkable progress in experimental techniques and the emergence of promising platforms for exploring novel aspects of quantum manybody phenomena. One notable example of many-body physics is the Kondo effect which arises from entanglement of the impurity spin with surrounding conduction electrons and the formation, below the Kondo temperature T K , of a spin singlet ground state [1]. In fact, the role of the electron spin can be replaced by any other quantum degree of freedom with symmetry protected two-fold degeneracy, e.g., orbital momentum [2] while the simultaneous presence of both a spin and an orbital degeneracy might lead to an SU(4) symmetric Kondo physics [3][4][5]. The SU(4) Kondo effect was already observed in carbon nanotubes, quantum dots with orbitally degenerate states, double quantum dot systems, and in a nanoscale silicon transistor [6][7][8][9][10][11].
Furthermore, the advent of scanning tunneling microscopy has made it possible to fabricate artificial Kondo nanostructures with atomic precision [12][13][14][15][16][17][18] whose properties, in particular the onset of lattice effects, have recently been a subject of increasing theoretical attention [19][20][21]. Magnetic atoms or organic molecules with orbital degeneracy deposited onto a metallic surface provide the opportunity to realize an SU(4) symmetric Kondo effect [22]. In addition, it is possible to study its evolution upon increasing the number of periodically arranged magnetic centers as in Ref. [23] where, start-ing from a single iron(II) phthalocyanine molecule deposited on top of Au (111) surface, in consecutive steps a two-dimensional superlattice was created followed by the theoretical analysis [24].
Yet another very active field of research with promises to provide new insights into the SU(N ) symmetric generalization of the Kondo effect [25][26][27] are quantum simulations with alkali-earth-like atoms in optical lattices [28]. Thus far, building on theoretical proposals [29][30][31], subsequent experimental studies utilizing ytterbium and strontium isotopes reported the observation of SU(N ) symmetric spin-exchange interactions between different orbitals with N as large as 10 [32][33][34][35]. From a practical point of view, there are three crucial features required for the realization of the SU(N ) Kondo physics in such setups: (i) the existence of a metastable excited state playing together with the atomic ground state the role of orbital degrees of freedom loaded into an orbital state-dependent optical lattice [35], (ii) a large nuclear spin I > 1/2 of fermionic isotopes must decouple from the electronic degrees of freedom to guarantee the SU(N = 2I + 1) spin-rotation symmetry of the interactions, and (iii) an antiferromagnetic (AF) character of spin-exchange interactions in the limit with one fully localized orbital. Although the currently used isotopes with I > 1/2 realize ferromagnetic interorbital interactions [32][33][34][35], ongoing theoretical proposals [36][37][38][39][40], numerical simulations [41][42][43], and utilizing other isotopes of alkali-earth-like atoms [44] allow one to envisage a controllable implementation of a Kondo-singlet state with SU(N ) symmetric interactions in single-impurity and lattice situations in the near future.
Given all these experimental developments, it is timely to consider an SU(N ) generalization of the conventional SU(2) Kondo lattice model (KLM) and to elucidate what kind of correlated phenomena occur under novel conditions with N > 2 which is the purpose of this paper. The importance of the KLM stems from its capability to account for the essential aspects of 4f -orbital-based heavyfermion systems summarized in the seminal Doniach phase diagram [45,46]. On the one hand, the Kondo exchange interaction J between the local moments and conduction electrons promotes a Kondo-screened paramagnetic phase in which the local moments are quenched by spins of the conduction electrons. On the other hand, the conduction electron-mediated RKKY exchange interaction between the local moments drives them into a magnetically ordered state thus leading to a quantum phase transition [47][48][49]. In some cases, the latter corresponds to a spin-density-wave transition as described in the Hertz-Millis scenario [50,51].
However, accumulating experimental evidence suggests that a realistic description of various types of quantum criticality and non-Fermi-liquid effects in heavyfermion systems requires, in addition to the "Kondo axis" K tuning the ratio between the Kondo temperature T K and the strength of RKKY interactions, a second "quantum axis" Q ∼ 1/S tuning the strength of quantum zeropoint fluctuations of the impurity spin S [52,53]. The magnitude of Q can be tuned by either increasing geometric frustration or reduction of dimensionality; large Q paves the way for the realization of exotic proposals such as local quantum criticality [54,55], fractionalized Fermi liquids [56][57][58], and partial Kondo screening [59][60][61].
Alternatively, when the physical SU(2) spin symmetry of the quantum model is generalized to SU(N ), a large number of degrees of freedom makes the long-range magnetic order less likely. An exciting aspect of studying SU(N ) quantum antiferromagnets in various geometries is that they allow one to pin down the role of Berry phases on the emergence of quantum-disordered ground state. As pointed out by Haldane [62], the relevance of the Berry phase term implies that point defects (hedgehogs) of the Néel field in spacetime acquire a net geometric phase. On the square lattice, large-N calculations predicted that the proliferation of topological defects in the presence of nontrivial Berry phases naturally leads to columnar q q q = (π, 0) valence bond solid (VBS) order in the paramagnet [63][64][65]. Later on, the onset of VBS order at sufficiently large N was confirmed by variational Monte Carlo study [66] and quantum Monte Carlo (QMC) simulations on the square [67][68][69][70][71][72] and honeycomb [73][74][75] lattices. Extensive QMC simulations of extended SU(N ) Heisenberg models have also led to insight into the nature of the quantum phase transition separating Néel and VBS phases [76][77][78][79] lending strong support to the theory of continuous "deconfined" quantum criticality [80,81]. Moreover, by extending QMC studies to the bilayer geometry [82], it has been confirmed that finite interlayer coupling renders Berry phases irrelevant at the quantum critical point [83][84][85][86]. As a result, the continuous Néel-VBS transition turns first order [82].
In contrast to SU(N ) Hamiltonians with direct effective spin-exchange interactions, very little is known about the phase diagram of the SU(N ) KLM with RKKY interactions between the impurity spins mediated by conduction electrons. On the one hand, coherent Kondo screening [87][88][89][90][91][92][93][94][95][96] and the formation of the Kondo insulating (KI) phase at half-filling can be accounted for within the large-N approach [97]. Strictly speaking, this mean-field approximation is formally justified only in a limit where the spin symmetry of the original model is extended from SU(2) to SU(N ) with N → ∞. Nevertheless, the method is often applied to heavy-fermion models with only SU(2) symmetry [19] and is considered as a good starting point for studying dynamical properties of heavy-fermion metals using the 1/N expansion technique [98][99][100][101]. However, there is no way of assessing a priori the validity of a large-N approach at any finite N .
On the other hand, despite a few attempts to develop a controlled treatment of both magnetism and the Kondo effect within a single large-N expansion [102,103], its applicability remains restricted to quantum disordered phases and thus the large-N approach cannot be used to explore the full phase diagram of the model. Another caveat of large-N approximation is that finite hybridization order parameter breaks the local gauge symmetry and implies that the constraint of single occupancy on the f -sites is fulfilled only on average which motivated the development of alternative approaches [104][105][106]. This yields a further motivation for systematic studies of the SU(N ) KLM using an unbiased method which handles the constraint numerically exactly so as to assess the validity of large-N approximate treatments [97].
Here, by performing auxiliary-field QMC simulations [107] we shall map out the phase diagram of the SU(N ) KLM as a function of the coupling J/t and the number of flavors N . Given diverse phenomena found upon loss of the AF order in SU(N ) Hubbard and Heisenberg models of magnetic insulators [67][68][69][70][71][72][73][74][75][76][77][78][79], one could equally expect the emergence of novel phases in the SU(N ) KLM. Furthermore, previous QMC simulations of the SU(2) KLM predict that below the magnetic energy scale T RKKY , the single-particle gap scales as J [108,109]. This contrasts with an exponentially small gap found in the dynamical mean-field theory [89], large-N limit [97], and Gutzwiller approximation [110], all of them omitting spatial fluctuations. Hence, we shall elucidate necessary conditions for recovering the large-N limit in the singleparticle dynamics thus providing a valuable benchmark of the large-N approach. one will show that the imaginary part of the action takes the value inπ with n an integer. Hence for even values of N the negative sign problem is absent.
In the large-N limit, the saddle-point approximation: z,λ=z * ,λ * = 0 (16) becomes exact. Assuming space and time independent fields produces the large-N mean-field theory discussed in Appendix A. The Monte Carlo method can hence be seen as not only accounting for all fluctuations around the large-N saddle point, but also for assessing if the saddle point is stable or not.
As mentioned above, our calculations were carried within the projective formulation. To be able to pull out N in front of the action, we use an SU(N ) symmetric trial wave function corresponding to the large-N saddle-point Hamiltonian: Unless stated otherwise, our simulations were carried out at finite imaginary time step ∆τ t = 0.1 and to generate the trial wave function, we have used V = 0.1t.

III. OVERVIEW OF THE RESULTS
The hybridization of conduction electron states with the f -electron states in a lattice situation leads to the large Fermi surface of the heavy-fermion metal and to the hybridization gap in the KI phase at half-band filling. The factors controlling the large Fermi surface topology continue to attract considerable attention [111][112][113][114][115][116][117][118][119][120][121].
As discussed in Appendix A, the number of flavors N is a control parameter which tunes the relative importance of the RKKY interaction and the Kondo energy scale. Here, we are interested in the following questions: (i) does the order-disorder phase transition exist for any N > 2 at all or just the opposite -one immediately reaches the large-N limit with only the KI phase? (ii) assuming that the phase transition continues to exist, is the continuous nature of the transition specific to the N = 2 case? (iii) does a larger N stabilize any intervening phase in the vicinity of the magnetic phase transition, e.g., VBS order? (iv) given that at the mean-field level with a frozen f -spin Ansatz, magnetic ordering and Kondo screening are not compatible [109], what are the single-particle spectral properties of the AF phase at finite N > 2? (v) Does the quasiparticle (QP) dispersion continue to feature a flat band extending up to k k k = (π, π) point signaling remnant Kondo screening of the impurity spins? If so, how does the QP residue evolve as a function of N and across the magnetic order-disorder transition?
A partial answer to these questions is given in Fig. 1 showing the ground-state phase diagram of the SU(N ) KLM as a function of the inverse of the number of fermion components N and Kondo coupling J. In addition, colorcoded circles correspond to the QP residue Z k k k of the doped hole at momentum (π, π). Our main result in Secs. IV A and IV E is that the SU(N ) KLM in the fully antisymmetric self-adjoint representation supports magnetic ordering for each considered value of N , and that no other phases aside from the Kondo insulator and Néel state intervene. Intuitively, we expect the J = 0 and N → ∞ point to be singular.
For the ordering of limits lim N →∞ lim J→0 we expect an AF ground state whereas for lim J→0 lim N →∞ a paramagnetic one. Figure 1 confirms this point of view: increasing N reduces critical coupling J c and enhances the domain of stability of the KI phase at the expense of the AF state.
At N = 2, and in the magnetically ordered phase, we observe a remarkable coexistence of Kondo screening and antiferromagnetism that stands at odds with the meanfield result predicting only a very narrow coexistence region [109,122]. We show in Sec. IV B that this does not carry over over to the larger values of N where we observe an abrupt drop in the QP residue Z (π,π) across the magnetic transition see Fig. 1. Finally, in Secs. IV C and IV D we investigate the impact of an enhanced Hamiltonian symmetry on the spin excitation spectrum and single-particle spectral function, respectively. We proceed now to discuss numerical results which led us to the above phase diagram.

IV. NUMERICAL RESULTS
Numerical results were obtained with an N -dependent projection parameter ranging from 2Θt = 50 for N = 2 to 2Θt = 400 for N = 8, chosen sufficiently large to guarantee the convergence to the ground state |Ψ 0 , see Appendix C 1. Physical observables have been extrapolated to the thermodynamic limit based on the QMC data obtained on lattice sizes ranging from 6 × 6 to 14 × 14 with periodic boundary conditions. Finite-size scalings and representative raw QMC data are presented in Appendices C 2, C 4, and C 3.

A. Spin degrees of freedom
We define the Néel state for the SU(N ) quantum antiferromagnet as: For the square lattice, we can split the lattice into two sub-lattices A and B such that the nearest neighbors of one sublattice belong to the other. For this Néel state, one will show that for i = j We hence adopt the following definition of the spin-spin correlation function, with α = c, f . To pin down the nature of ground state of the SU(N ) KLM, we compute the staggered moment, The corresponding finite-size scaling is presented in Appendix C 2 and the extrapolated values for localized (conduction) electrons are plotted versus J/t in Fig. 2(a) [ Fig. 2(b)], respectively. On the one hand, the QMC data confirms that increasing N suppresses magnetism by shifting the magnetic order-disorder transition from J c /t 1.44(1) in the SU(2) symmetric case to progressively lower values of J/t: for N = 4 we find a critical point at 1.09(1); for N = 6 and N = 8 we find the critical points at 0.85(1) and 0.73(1), respectively. In this respect, the effect of finite N bears a similarity to that generated by geometric frustration, e.g., by next-nearestneighbor hopping t [112,117]. While the RKKY scale varies as 1 N , the Kondo scale is N independent such that comparing scales yields J c (N ) ∝ 1 log N in the large-N and small J limit (see Appendix A). We have used this form to plot a guide to the eye for the phase boundary in Fig. 1. On the other hand, the data in Fig. 2(a) suggests that in the magnetically ordered phase, the f -local moment m f remains large since up to N = 8 it exceeds 80% of the Néel value. Furthermore, while m f seems to grow continuously below J c at N = 2 and N = 4, one finds a rapid buildup of the f -local moment for larger N .
Once the magnetic order disappears at J c , the ground state is expected to develop a finite spin gap ∆ s (q q q), i.e., the energy difference between the singlet S = 0 ground state and the lowest excited spin triplet S = 1 state with momentum q q q. To compute ∆ s (q q q) without resorting to analytic continuation, we consider the imaginary-time displaced spin correlation functions, whereŜ(q q q, τ ) =Ŝ f,µ ν (q q q, τ ) +Ŝ c,µ ν (q q q, τ ) is the total spin. The spin gap ∆ s (q q q) can be extracted from the asymptotic behavior of S(q q q, τ ) at τ t 1 since S(q q q, τ ) ∝ exp (−τ ∆ s (q q q)). Here, we focus on the AF wavevector Q Q Q = (π, π), i.e., ∆ s ≡ min q q q ∆ s (q q q) = ∆ s (Q Q Q). A linear extrapolation of finite-size QMC estimates of ∆ s (N ) to the thermodynamic limit, see Appendix C 3, leads to the results plotted in Fig. 2(c). For each N we find that opening of the spin gap coincides with the vanishing of the magnetic moment. As shown in Fig. 2(c), upon increasing N the J-dependence of the spin gap approaches asymptotically the large-N behavior ∝ e −t/J . While the spin gap evolves smoothly across the transition, the magnetization shows an abrupt change, especially at larger values of N . This poses the question of the nature of the transition, first-order or continuous. To provide more insight, we plot in Fig. 3 the free-energy derivative, obtained on our largest 14 × 14 lattice. The smooth evolution of this quantity across the critical point suggests that the phase transition remains continuous for each N .

B. Charge degrees of freedom
An important quantity of direct experimental relevance is the QP residue Z k k k . Indeed, since the effective QP mass m * ∝ 1 Z k k k , the behavior of Z k k k along the Fermi surface reveals how electron interactions modify properties of a metal. Given that QMC simulations are restricted to the half-filled case, one possibility to get insight into properties of the metallic state at small doping is to consider the problem of a single-hole doped into the insulating phase. Then, assuming a rigid band scenario, one can , directly from the long-time behavior of the imaginary-time Green's function: Here, E n 0 is the ground-state energy at half-filling with n particles while E n−1 0 (k k k) corresponds to the energy eigenstate with momentum k k k in the n − 1 particle Hilbert space. Typical raw data of G(k k k, τ ) from QMC simulations with different system sizes L and the extrapolation to the thermodynamic limit of finite-size estimates of ∆ qp ≡ min k k k ∆ qp (k k k) = ∆ qp (k k k = (π, π)) and Z (π,π) is presented in Appendix C 4. Figure 4 plots ∆ qp and Z (π,π) as a function of J/t. We first discuss the evolution of these quantities in the paramagnetic phase. As shown in the insets of Fig. 4 both quantities evolve smoothly from N = 2 to N = ∞. The fact that we are able to recover the large-N results by extrapolating QMC data obtained by handling the constraint of no double occupancy on the localized felectron orbitals numerically exactly, validates the large-N approximate treatments of the constraint and confirms that the large-N theory is the correct saddle point of the SU(2) KLM. Furthermore, by comparing the Ndependence of the single-particle gap ∆ qp in Fig. 4(a) with that of the spin gap ∆ s in Fig. 2(c), we conclude that upon increasing N both quantities evolve in the paramagnetic phase (i.e. J/t = 2) to the asymptotic limit ∆ s = ∆ c with ∆ c = 2∆ qp being the charge gap, corresponding to the band insulator in the noninteracting periodic Anderson model. Our results hence provide a text-book numerical demonstration that the N = 2 Kondo lattice in the paramagnetic phase is adiabatically connected to the large-N saddle point.
Across the magnetic transition, ∆ qp and Z (π,π) show a very strong N dependence. In contrast to the N = 2 case with a smooth evolution of both quantities across the quantum critical point, a nonmonotonic behavior of the single-particle gap accompanied by an abrupt reduction of the QP weight on the magnetically ordered side is apparent. Although Z (π,π) shows an abrupt change, it remains finite in the magnetic phase. Hence we find the continued existence of the heavy-fermion band for all the (a) Single-particle gap ∆qp at momentum k k k = (π, π) and (b) the corresponding QP residue Z (π,π) after extrapolation the QMC data to the thermodynamic limit, see Appendix C 4. Insets show second-order polynomial fits to the QMC data in order to extract ∆qp and Z (π,π) in the N → ∞ limit; the extrapolated values ∆ N →∞ qp = 0.109(3) and Z N →∞ (π,π) = 0.0256(3) at J/t = 2 and ∆ N →∞ qp = 0.051(2) and Z N →∞ (π,π) = 0.0125(2) at J/t = 1.6 match well those obtained using the large-N approximation (triangles).
considered values of N down to the smallest J. Assuming a rigid-band scenario this implies that, in contrast to an Ansatz with frozen f -spins, the emergent heavyfermion metal at small coupling J is characterized by a large Fermi surface containing both conduction and localized electrons (see Sec. V ). As a consequence, the coherence temperature is expected to drop abruptly across the transition. On the magnetically ordered side, the QP gap tracks J/N in the small J limit. Finally, a notable feature in Fig. 4(a) is a broad plateau in the J-dependence of ∆ qp at N = 4. It is interesting to point out that a similar plateau was found in quantum cluster theories allowing for SU(2) symmetry-breaking AF order [117,121] as well as in the bond fermion theory [123].
To provide a theoretical framework for the above, we introduce in Appendix B a mean-field theory. Here we describe how this mean-field theory and fluctuations around the corresponding saddle point can account for the QMC results. Our numerics shows that for any fixed value of N the paramagnetic state is unstable to an RKKY driven magnetic instability and that deep in the magnetic phase the f -local moment is next to saturated. The Néel state of Eq. (18) is hence a good starting point to formulate a mean-field theory. This wave function breaks the U(N ) symmetry down to U(N /2) × U(N /2). The mean-field Hamiltonian derived in Appendix B possesses a U(N /2) × U(N /2) symmetry and is a generalization of the mean-field theory of Ref. [122] that captures both Kondo screening and magnetism to the SU(N ) group. In the mean-field Hamiltonian the RKKY interaction scales as 1/N . As a consequence, and owing to the nesting properties of the conduction-electron band, the magnetically induced QP gap scales as J/N . Our QMC results support this.
Following Ref. [124] one can define a model Hamiltonian that reduces to the KLM model at N = 2, that has the U(N /2) × U(N /2) symmetry beyond N = 2, and that reproduces the saddle point of Eq. (B11) in the large-N limit. It is very tempting to interpret our QMC results in terms of fluctuations -that are suppressed as a function of N -around this magnetic saddle point. In the limit N → ∞ [109], and deep in the magnetically ordered phase, the f -spins are frozen and Z (π,π) vanishes.

C. Spin excitation spectrum
To get further insight into the nature of AF and KI phases in the SU(N ) symmetric situation, we consider the dynamical spin structure factor S(q q q, ω). We have extracted this quantity from the QMC imaginary-time displaced spin correlation functions defined in Eq. (22), by using the stochastic analytic continuation method [125]. In the above, we consider the total spin. The spin-density-wave approximation presented in Appendix B breaks explicitly the SU(N ) symmetry and hence it fails to capture Goldstone modes. Since the Néel state has the U(N/2) × U(N/2) symmetry but the Hamiltonian a U(N ) one, we expect a total of dim [126,127] that should become apparent in the dynamical spin structure factor. One expects that this large number of Goldstone modes will destabilize the ordering and in this respect, it is interesting to see that in the KLM an AF state can be stabilized for each N . Concerning the energetics, and as argued in Appendix A, the effective RKKY coupling J RKKY scales as J 2 N ( f ) N and the single-particle gap as J N . Hence in the small J limit, the Goldstone modes are expected to be located well below the particle-hole continuum.
in the AF phase with J/t = 0.8 is clearly seen for N > 2. For each N considered in Figs. 5(a-c), the particle-hole continuum lies above the Goldstone modes such that we should interpret the data solely in terms of an SU(N ) quantum spin model. Adopting this point of view, the relevant energy scale is the spin-wave velocity that at fixed J scales as 1/N . In terms of this energy scale, it becomes apparent that the width of the dynamical spin structure factor at say wavevector q = (0, π) grows as a function of N . We interpret this as a consequence of scattering between a growing number of Goldstone modes. One should also mention that as a function of growing values of N , the distance to the magnetic critical point is suppressed. Although for all considered values of parameters the magnetic moment is well developed, this could certainly play a role in the interpretation of the N -dependence of the spectrum. Figures 5(d-f) plot the dynamical spin structure factor at J/t = 1.6 as a function of N . At N = 2, we are close to the quantum critical point, and the triplon mode shows a minimal gap at the AF wavevector, Q Q Q = (π, π). Triplons will condense at the transition to generate the magnetic order. In this regime triplons are bound electronhole pairs and the binding originates from vertex corrections. Enhancing N from this point, damps vertex corrections such that the bound triplon mode will progressively merge in the particle-hole continuum. Precisely this effect is seen in Figs. 5(d-f).

D. Single-particle spectral function
We move on to discuss in more detail the single-particle dynamics. To this end, we have computed the singleparticle spectral function A(k k k, ω) of the conduction electrons. It is related to the imaginary-time Green's function defined in Eq. (24) via: Again, we use the stochastic analytic continuation method to extract A(k k k, ω). The evolution of A(k k k, ω) upon increasing N in the AF phase at J/t = 0.8 is shown in Figs. 6(a-c). As expected for the half-filled case, all the spectra display a clear hybridization gap which, in agreement with findings in Sec. IV B, becomes gradually smaller at larger N . Furthermore, the spectral function features a flat heavyfermion band extending to the k k k = (π, π) point with relatively low spectral weight. The continued presence of this band around k k k = (π, π) even in the AF phase shows that the heavy fermions undergo a magnetic instability such that Kondo screening is still present in the ordered phase.
A direct consequence of the magnetic ordering is a back-folding of the Brillouin zone and the emergence of additional low-energy spectral feature around the k k k = (0, 0) momentum. It arises due to the scattering of the heavy QP off spin fluctuations with the AF wavevector Q Q Q = (π, π) and thus it corresponds to the shadow of the band in the vicinity of the k k k = (π, π) point. As apparent, the shadow band becomes less pronounced at larger N . This can be traced back to the combined effects that the Z-factor at k k k = (π, π) drops as a function of N and that the magnetic moment is reduced as N grows from 2 to 6 at J/t = 0.8. In Figs. 6(d-f) we show the evolution of A(k k k, ω) upon increasing N in the KI phase at J/t = 1.6. In the disordered phase with N = 2, only a precursive feature of the shadow band is visible, see Fig. 6(d): despite a much larger spectral weight of the heavy QP band at k k k = (π, π) point with respect to J/t = 0.8, the precur- sive feature has relatively low intensity and it is shifted by the energy corresponding roughly to the spin gap ∆ s . As shown in Figs. 6(e,f), this feature becomes broad and consequently more difficult to resolve at larger N . This can be traced back to the fact that as a function of N , the triplon mode approaches the particle-hole continuum broadens, and ultimately disappears.
Finally, in Fig. 7 we plot A(k k k, ω) in the KI phase at J/t = 1.6 for our largest N = 8 together with that obtained in the large-N approach. As apparent, the large-N approximation produces a single-particle spectrum which compares favorably with the QMC spectral function. One of the key properties of the large-N selfenergy, is its locality: Thereby, despite all the caveats of the large-N approximation -finite hybridization order parameter which breaks the local gauge symmetry -it can be considered to be well suited to account for the essence of Kondo screening deep in the KI phase.

E. VBS correlation function
Generically, enhancing the symmetry group from SU(2) to SU(N ) leads to VBS orders. To confirm the absence of this instability in the SU(N ) KLM, we have  computed the VBS correlation function for the f -spins:  Figs. 8(b,d)]. As expected, the VBS correlation function S f VBS (q q q) for N = 2 is featureless and lattice-size independent throughout the transition confirming that the SU(2) KLM is dominated by magnetic fluctuations. Given numerical evidence for enhanced columnar q q q = (π, 0) dimer correlations in SU(N ) Hubbard and Heisenberg models [66][67][68][69][70][71][72][73][74][75], one could expect that the same physics shows up in the SU(N ) KLM. In contrast, even though the lineshape of S f VBS (q q q) gets sharper at N = 4, a dominant cusp feature is found at the AF wavevector Q Q Q = (π, π), see Figs. 8(b,d). Same behavior is observed across the magnetic order-disorder transi-tion for larger N = 6 and N = 8 shown in Figs. 9(a,c) and 9(b,d), respectively. Thus, we conclude that there are no significant columnar dimer fluctuations in the phase diagram and we interpret the cusp feature at Q Q Q = (π, π) as a fingerprint of the perfectly nested conduction-electron Fermi surface in the noninteracting limit.

A. Experimental relevance
Heavy fermion systems are prototype materials to study quantum criticality of the magnetic order-disorder transition. Given the complexity of the problem, theoretical and experimental studies on the quantum criticality in heavy-fermion systems explore various routes to approach the quantum critical point (QCP) [47][48][49]. One possibility is to modify the strength of the Kondo coupling, e.g., by varying chemical or external pressure. Another route is to tune intersite interactions between the f -moments by considering systems with different dimensionality or with geometrical frustration.
In some materials, e.g., Ce 1−x La x Ru 2 Si 2 [128], the data are consistent with predictions of the conventional spin-density-wave theory [50,51], which considers the f -electrons as itinerant on both sides of the QCP. In this case, the dominant critical AF fluctuations modify neither the shape nor size of the large Fermi surface which incorporates both conduction electrons and the f -electron states. In other compounds such as CeCu 6−x Au x [129] and YbRh 2 Si 2 [130][131][132], there are indications for the breakup of composite heavy-fermion QPs and the concomitant collapse of the large Fermi surface driven by local critical magnetic fluctuations [54,55]. Moreover, the reconstruction of the Fermi surface may also occur away from the QCP -within the magnetically ordered phase [133] or even more exotically -outside [134], paving the way for an intervening phase where the local f -moments are neither Kondo screened nor antiferromagnetically ordered.
Here, in order to gain novel insight into the quantum criticality in heavy-fermion systems, we have considered the SU(N ) generalization of the KLM. Given that increasing N changes the degree of quantum fluctuations of the local f -moments, it allows one to investigate the impact of magnetic fluctuations on the coherent Kondolattice formation in a single setup. For each N , we find that quantum critical phenomena are driven by fluctuations of the AF order parameter. Importantly, we do not observe breakdown of Kondo screening which continues to exist on the magnetically ordered side of the phase diagram.
However, our findings show that increasing N strongly modifies the behavior of the QP residue Z (π,π) across the QCP. As such they have important implications for the interpretation of experimental data. Considering that in reality experiments are performed at small but finite temperatures, a rapid decrease of the QP residue resolved for N = 8, see Fig. 4(b), could be easily mistaken in the isothermal measurement of Hall coefficient as that in Ref. [131], as evidence for a collapse of the large Fermi surface at the QCP.
B. Quantum-field-theoretic perspective Throughout the J-N plane, the charge degrees of freedom are gapped. Hence charge fluctuations around halffilling -that mix SU(N ) spin representations -will not contribute in the low-energy effective field theory and can be safely omitted. The remaining degree of freedom is an SU(N ) spin in the totally antisymmetric representation corresponding to a Young tableau consisting of a single column and N /2 rows. Since we observe AF phases, the low-energy effective model is that of an SU(N ) quantum antiferromagnet: in the aforementioned representation. The generalization of Haldane's SU(2) spin coherent state path integral formulation [62] to the SU(N ) group has been carried by Read and Sachdev [63]. It is beyond the scope of this paper to review the derivation, and we will only cite the final result. As for the SU(2) case, the SU(N ) spin coherent state |q is obtained by an SU(N ) rotation of the Néel state [63,135]. It satisfies the relation: and the ± sign refers to the A and B sublattices. Q µ ν (i) hence corresponds to the Néel order parameter, that owing to the sign convention is uniform in space, and whose low-energy fluctuations are governed by the action: with Berry phase S B [63] and non-linear σ (NLσ) model, (32) In the above ρ s corresponds to the spin stiffness and c to the velocity. For the SU(2) case, we can write Q = n · σ with n a unit vector and σ the vector of Pauli spin matrices. With this parameterization, the above reduces to the well known O(3) NLσ model with Berry phase. In contrast to the one-dimensions, smooth space-time variations of the Néel order parameter have a vanishing Berry phase [62]. For the above U(N ) model, the order parameter manifold corresponds to U (N ) U (N/2)×U (N/2) . Since the second homotopy group of this space is given by Z, skyrmions are well defined, and one can carry over the arguments put forward by Haldane for the SU(2) case. In particular, skyrmion number changing field configurations (hedgehogs) carry a non-trivial Berry phase and quadruple hedgehogs insertions carry no Berry phase and hence do not interfere destructively. This suggest that the Hilbert space splits into four distinct classes corresponding to the skyrmion number modulo 4. Proliferation of quadruple hedgehog configurations has been argued to correspond to the VBS state [81] and is the essence of the notion of deconfined quantum criticality.
With this background that links the Berry phase to a four-fold degenerate VBS state and the lack of any (0, ±π) and (±π, 0) singularities in the VBS order parameter across the order-disorder QCP in the QMC data, see Sec. IV E, we conclude that the Berry phase can be omitted in the effective field theory. A similar result holds for the SU(2) bilayer Heisenberg model [83][84][85][86]. The SU(N ) KLM hence provides a lattice realization of NLσ model of Eq. (32). To the best of our knowledge, the critical exponents as well as the very nature of the transition as a function of N are unknown. A 1/N expansion study of the critical exponents has been carried out in Ref. [136] for the general representation corresponding to Young tableau of m (N − m) rows and one column on sublattice A (B). As pointed out in the paper [136], the results require m/N to be a small number and cannot be carried over to the self-conjugate representation where m = N/2.

VI. CONCLUSIONS
Our major findings can be summarized by the following points.
(i) A serious caveat of the large-N approximation is that it introduces a finite hybridization order parameter. It breaks the local gauge symmetry and implies that the constraint of single occupancy on the f -sites is fulfilled only on average [97]. Here, we have handled the constraint of no double occupancy numerically exactly with QMC simulations. By extrapolating finite-N QMC data to the N → ∞ limit, we were able to recover the large-N results in the KI phase. This validates large-N approximate treatments of the constraint and confirms that the large-N theory is the correct saddle point of the SU(2) KLM.
(ii) Up to N = 8 we observe a magnetically ordered phase. The RKKY interaction scales as 1/N and the Kondo energy is N -independent such that matching the two energy scales gives J c (N ) ∝ 1 log(N )N ( F ) . This form is consistent with our data. Since the charge degrees of freedom are gapped throughout the phase diagram, the magnetically ordered state should be understood in terms of an SU(N ) quantum antiferromagnet on a bilayer square lattice. Let us consider the representations discussed in Ref. [65], consisting of a Young tableau of m (N − m) rows and one column on sub-lattice A (B). The Néel bro-ken symmetry phase of the model has Lorentz symmetry and accordingly 2N m − 2m 2 Goldstone modes [126,127]. This count matches the dimension of the manifold on which the NLσ model of Sec. V B is defined. The number of Goldstone modes is a measure of the fluctuations around the Néel state and is maximal for the representation m = N/2 considered here. It is hence interesting to compare our result to that of Ref. [82] for the SU(N ) bilayer Heisenberg model with nearest-neighbor couplings at m = 1: at N = 8 no magnetic ordering is present. To reconcile this apparent contradiction we have to take into account the range of the RKKY interaction. To a first approximation, it is given by the inverse single-particle gap in the KI phase at a value of J just above J c (N ). With the above form for J c (N ) and large-N form for the single-particle gap, ∆ qp ∝ e − 1 JN ( F ) , we find that the range of the RKKY interaction grows as a power of N . We believe that this enhanced range of the interactionvery specific to the KLM -is the key to stabilize antiferromagnetism at large N .
(iii) We have argued in Sec. V B that the Berry phase could be omitted, such that the KLM provides a unique possibility to study the critical phenomena associated with the NLσ model of Eq. (32) at m = N/2. To the best of our knowledge, and as discussed in Sec. V B, this universality class has never been studied.
(iv) Since the range of the RKKY interaction grows as a function of N we expect that the interplay between charge and spin degrees of freedom will become more mean-field-like. In fact, at large N we observe an abrupt reduction of the QP residue Z (π,π) upon entering the AF phase. This behavior is very reminiscent of that observed in mean-field calculations of the SU(2) KLM that take into account both antiferromagnetism and Kondo screening [109]. Within a rigid band shift assumption to describe the heavy-fermion metallic state at small doping, this means that the coherence temperature drops by up to an order of magnitude across the magnetic transition. Isothermal measurement of the Hall coefficient that are below the coherence temperature in the paramagnetic heavy-fermion phase and above it in the magnetically ordered state, would be interpreted as a breakdown of Kondo screening [131]. Nevertheless, we find the signature of the heavy-fermion band for all the considered values of N down to the smallest J. With the aforementioned rigid band shift, the emergent heavyfermion metal at small coupling J is characterized by a large Fermi surface containing both conduction and localized electrons. In the magnetically ordered phase, back-folding of the Fermi surface accounts for the reduced translation symmetry. This abrupt reduction of Z (π,π) upon entering the AF phase at large N could be interpreted as a sign of a first-order transition. We however did not observe the typical signs -level crossingof a first-order transition in our simulations.
and the N 2 − 1 generators of SU(N ) satisfy the normalization condition: To estimate the energy scale of the RKKY interaction, we will first consider a single impurity at the origin with a frozen f -spin: Within first order perturbation theory in J the frozen fspin at the origin produces ripples in the spin texture that follow the spin susceptibility of the conduction electrons, Here, with (k) = −2t (cos(k x ) + cos(k y )), f (k) = 1 1+e β (k) and χ c (r, iΩ m = 0) = 1 L 2 q e −iq·r χ c (q, iΩ m = 0). We now consider a second impurity at position r that Kondo couples to the conduction electrons according to Eq. (A1). At the mean-field level, the interaction energy between the two spins, reads: T a,f r r r . (A7) Comparing the above expression to the RKKY Hamiltonian: that describes the effective SU(N ) Heisenberg interaction between the impurity spins gives: Hence, the RKKY interaction measured relative to the kinetic energy scales as 1 N .

The Kondo scale
In contrast, we now argue that the Kondo scale is Nindependent in the large-N limit. To formulate the large-N mean-field saddle-point, we use the completeness relation, a T a α,β T a α ,β = to show that up to a constant: Using the mean-field Ansatz V = 2 D † i /N and imposing the constraint N σ=1f † i,σfi,σ = N 2 on average yields the gap equation: where ∆ = JV /2, E = 2 2 + ∆ 2 , N ( ) = 1 L 2 k δ( − (k)), and f the Fermi function. For the particle-hole symmetric case considered here, the f -electron half-filling constraint is satisfied by symmetry such that no Lagrange multiplier has to be introduced. At the Kondo temperature, T K , ∆ vanishes, such that T K is given by: In the above, we have used the particle-hole symmetric condition N ( ) = N (− ). As apparent, the above equation is independent on N such that at the mean-field level, the Kondo temperature does not scale with N . Finally we note that for a density of states of width W and in the small J limit, 3. Functional form of the critical coupling Jc(N ) We can now compare scales to estimate the the critical value of J c (N ) where the Kondo effect gives way to magnetic ordering: In the above, we have used N ( f ) = 1 W for the aforementioned flat density of states, and for instance, considered the value of the spin-susceptibility at a distance given by the lattice spacing. In the large-N limit where we expect J c to be small, we obtain: Appendix B: Spin-density-wave approach for the SU(N ) KLM The data presented in the main text suggests that for each N magnetism and Kondo screening coexist, and that in the magnetically ordered phase, the f -local moment is large since up to N = 8 it exceeds 80% of the Néel value. In this appendix, we generalized the mean-field theory of Ref. [122] that captures both Kondo screening and magnetism to the SU(N ) group. To do so, we consider the following explicit form of the SU(N ) generators. For α > β included in the set of [1, N ] we consider the N 2 −N off-diagonal generators: alongside the N − 1 diagonal operators: (B2) In the above, n runs from 1, . . . , N − 1. This definition of the SU(N ) generators satisfies the normalization condition of Eq. (A3) and similar forms hold for the felectrons. As in Ref. [122], the off-diagonal operators will account for Kondo screening whereas the diagonal ones for magnetism. With the above, the SU(N ) Kondo interaction reads: To account for the Kondo effect, we adopt the mean-field Ansatz: For the magnetism, it is convenient to carry out an orthogonal transformation of the diagonal operators: such that: Identical forms hold for the f -electrons. In the Néel state, we have, with the AF wave vector Q = (π, π). This motivates the Ansatz: The above Ansätze break the U(N ) symmetry down to a U(N /2) × U(N /2) one such that it becomes convenient to introduced the notation: with µ = 1, . . . , N/2 and σ = ±1. The mean-field Hamiltonian is then given by: Here, (k) = −2t (cos(k x ) + cos(k y )) such that (k + Q) = − (k) and particle-hole symmetry pins the average f -occupation to N/2. The saddle-point equations then read: with F = − 1 β log Tre −βĤ M F . Several comments are in order.
• The underlying particle hole-symmetry pins the foccupation to half-filling such that no Lagrange multiplier is required to enforce this constraint on average.
• Since the µ index does not appear in the Hamiltonian matrix, the above has a U(N /2) × U(N /2) symmetry that generalizes of the U(1) × U(1) symmetry presented in Refs. [109,122]. One will notice that at N = 2 we recover precisely Eq. (45) of Ref. [109]. In this case, and assuming m f = m c = 0 but V = 0 as appropriate for the KI phase, one finds the single-particle dispersion relation: QP gap ∆ qp = −E − k=(π,π) and residue: Z (π,π) = 1 2 1 − (k = (π, π)) 2 (k = (π, π)) + (JV ) 2 , (B14) for a doped hole away from half-filling. Solving selfconsistently the saddle-point equation for the hybridization order parameter V and using the above relations for ∆ qp and Z (π,π) lead us to the large-N results shown in Figs. 4 and 7(b) in the main text. On the other hand, assuming that the AF order is present m f = 0 and m c = 0 and the spin degrees of freedom are frozen such that V = 0, the corresponding dispersion relation reads: and the QP gap tracks J/N as does the QMC data in Fig. 4(a) in the main text.
• While the magnetic energy scales as order N 0 the kinetic and hybridization energies scale as oder N . This can be seen explicitly in the last constant term of Eq. (B11) and is consistent with the above discussion of the Kondo and RKKY energy scales. As a consequence, we expect the J = 0 and N → ∞ point to be singular. For the ordering of limits lim N →∞ lim J→0 we expect an AF ground state whereas for lim J→0 lim N →∞ a paramagnetic one.
• It is very tempting to follow ideas presented in Ref. [124] and to formulate a U(N /2) × U(N /2) field theory that possesses the above mean-field Hamiltonian as a saddle point in the large-N limit and that reproduces the U(2) invariant KLM model at N = 2. At N = 2, and in the magnetically ordered phase, the we observe a remarkable coexistence of Kondo screening and antiferromagnetism that stands at odds with the mean-field results predicting only a very narrow coexistence region [109,122]. As a function of N , fluctuations around the magnetically ordered saddle point are reduced and we expect a stronger mean-field-like competition between magnetic ordering and Kondo screening. It is very interesting to see that the QMC data supports this line of thought.

Appendix C: Supplemental data
Here we provide further details about the QMC simulation results discussed in the main text.

Convergence to the ground state
In this appendix we check the dependence of the QMC results on the projection parameter Θ. In order to ensure that a given result corresponds to the ground state we have performed test simulations on the 12 × 12 system at a variety of projection parameters Θ. The energy scales of the KLM, the single-ion Kondo temperature, coherence temperature, and the RKKY scale, they all become smaller on decreasing J/t. The calculations become more expensive in the SU(N ) case since as shown in Appendix A, the RKKY scale ∝ 1 N . Consequently, increasingly large projection parameters are required to reach the AF ground state and the issue becomes partic-ularly severe for small values of J, see Figs. 10(a,d,g,j).

Spin structure factor
As discussed in Sec. IV A, we have estimated the onset of long-range magnetic order from the behavior of the   staggered magnetic moment: extracted separately for the f -and c-electrons. The corresponding finite-size scaling analysis of the spin structure factor S α (Q Q Q) is shown in Figs. 10(b,e,h,k) and 10(c,f,i,l), respectively. Long-range AF order is present when lim L→∞ S α (Q Q Q)/L 2 extrapolates to a finite value.

Spin gap
In Sec. IV A, the gap for spin excitations ∆ s (q q q) was obtained by considering the imaginary-time displaced spin correlation functions, S(q q q, τ ) = µ,ν Ŝ µ ν (q q q, τ )Ŝ ν µ (−q q q) ,  whereŜ(q q q, τ ) =Ŝ f,µ ν (q q q, τ ) +Ŝ c,µ ν (q q q, τ ) is the total spin. The spin gap ∆ s for a given linear system size L has been extracted from the asymptotic behavior of S(Q Q Q = (π, π), τ ) ∝ exp (−τ ∆ s ) at large imaginary time τ . Extrapolating to the thermodynamic limit, one finds for each N that the spin gap scales to a finite value in the KI phase and vanishes inside the AF state due to the emergence of Goldstone modes of the broken continuous SU(N ) symmetry group, see Figs. 11(a,d,g,j).

Single-particle dynamics
As described in Sec. IV B, to probe the single-particle dynamics we have measured the imaginary-time displaced Green's function for the conduction electrons, G(k k k, τ ) = 1 N σ Ψ n 0 |c † k k k,σ (τ )c k k k,σ |Ψ n 0 . (C3) The single-particle gap ∆ qp at momentum k k k = (π, π) and the corresponding QP weight Z (π,π) were extracted by fitting the tail of the Green's function to the form Z k k k e −∆qp(k k k)τ . As an example, Fig. 12 shows raw data of G[k k k = (π, π), τ ] for N = 4, 6, and 8 obtained from QMC simulations with different system sizes L at J/t = 0.6. The good quality of the data allowed us to determine finite-size estimates of ∆ qp and Z (π,π) directly on the imaginary-time axis. The corresponding extrapolation of both quantities to the thermodynamic limit is performed in Figs. 11(b,e,h,k) and 11(c,f,i,l), respectively. Note enhanced finite-size effects in vicinity of the quantum critical point.

Imaginary-time discretization ∆τ
In Sec. IV E, we have calculated the VBS correlation function S f VBS (q q q). We used the imaginary-time step ∆τ t = 0.1 in the discrete Trotter-Suzuki decomposition in Eq. (9) which yields an error O(∆τ 2 ). In order to exclude that the cusp feature at the AF wavevector Q Q Q = (π, π) is an artifact related to the Trotter-Suzuki decomposition, we have repeated QMC simulations for N = 4, 6, and 8 with a twice smaller imaginary-time step ∆τ t = 0.05. The corresponding dimer structure factors S f VBS (q q q) shown in Fig. 13 look qualitatively very similar to those in Figs. 8(b,d) and 9 in the main text.