Exactly Solvable Kondo Lattice Model in Anisotropic Limit

TL;DR

This paper introduces an exactly solvable Kondo lattice model that captures magnetic and electronic phases, including antiferromagnetic insulators and metallic states, with potential relevance to heavy fermion and hidden order materials.

Contribution

The model is exactly solvable without fine-tuning, describing itinerant electrons coupled to localized moments via longitudinal exchange, valid for arbitrary lattice geometry and filling.

Findings

Ground state is a Néel antiferromagnetic insulator at half-filling.

Finite temperature phases include a Mott insulator and a correlated metal.

Transition to antiferromagnetic insulator is a continuous 2D Ising transition.

Abstract

In this paper we introduce an exactly solvable Kondo lattice model without any fine-tuning local gauge symmetry. This model describes itinerant electrons interplaying with a localized magnetic moment via only longitudinal Kondo exchange. Its solvability results from conservation of the localized moment at each site, and is valid for arbitrary lattice geometry and electron filling. A case study on square lattice shows that the ground state is a N\'{e}el antiferromagnetic insulator at half-filling. At finite temperature, paramagnetic phases including a Mott insulator and correlated metal are found. The former is a melting antiferromagnetic insulator with a strong short-range magnetic fluctuation, while the latter corresponds to a Fermi liquid-like metal. Monte Carlo simulation and theoretical analysis demonstrate that the transition from paramagnetic phases into the antiferromagnetic insulator is a continuous $2D$ Ising transition. Away from half-filling, patterns of spin stripes (inhomogeneous magnetic order) at weak coupling, and phase separation at strong coupling are predicted. With established Ising antiferromagnetism and spin stripe orders, our model may be relevant to a heavy fermion compound CeCo(In$_{1-x}$Hg$_{x}$)$_{5}$ and novel quantum liquid-crystal order in a hidden order compound URu$_{2}$Si$_{2}$.

Equations

$$\hat{H}$$

$$\hat{H}(q)=\sum _{i,j,\sigma }{t}_{ij}{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+\sum _{j\sigma }\frac{J\sigma }{4}{q}_j{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{j\sigma }$$

$$\hat{H}$$

$$G_{\sigma }(k,\omega )=\frac{1}{\omega -\frac{J^2}{16\omega }-{\epsilon }_k}=\frac{{\alpha }_k^2}{\omega -{\tilde{E}}_k^{+}}+\frac{1-{\alpha }_k^2}{\omega -{\tilde{E}}_k^{-}}.$$

$$S$$

$$\mathcal{Z}=\mathrm{Tr}{e}^{-\beta \hat{H}}={\mathrm{Tr}}_c{\mathrm{Tr}}_S{e}^{-\beta \hat{H}}=\sum _{\{q_j\}}{\mathrm{Tr}}_c{e}^{-\beta \hat{H}(q)}.$$

$$\hat{H}(q)=\sum _n{E}_n{\hat{d}}_n^{†}{\hat{d}}_n$$

$${\hat{c}}_j=|0⟩⟨j|=\sum _n|0⟩⟨j|n⟩⟨n|=\sum _n|0⟩⟨n|⟨j|n⟩=\sum _n{\hat{d}}_n{\varphi }_n^j.$$

$${\hat{c}}_j^{†}=|j⟩⟨0|=\sum _n|n⟩⟨n|j⟩⟨0|=\sum _n|n⟩⟨0|⟨n|j⟩=\sum _n{\hat{d}}_n^{†}({\varphi }_n^j{)}^{\ast }.$$

$${\mathrm{Tr}}_c{e}^{-\beta \hat{H}(q)}$$

$$\mathcal{Z}=\sum _{\{q_j\}}\prod _n(1+e^{-\beta E_n(q)})=\sum _{\{q_j\}}{e}^{-\beta F(q)}$$

$$F(q)=-T\sum _n\ln (1+e^{-\beta E_n(q)}).$$

$$r=\frac{e^{-\beta F_{new}(q)}}{e^{-\beta F_{old}(q)}}=e^{-\beta (F_{new}(q)-F_{old}(q))}$$

$$\tilde{r}=\frac{r}{1+r}$$

$$\hat{O}={\hat{O}}^c+{\hat{O}}^q.$$

$$⟨\hat{O}⟩=⟨{\hat{O}}^c⟩+⟨{\hat{O}}^q⟩=\frac{\mathrm{Tr}{\hat{O}}^c{e}^{-\beta \hat{H}}}{\mathrm{Tr}{e}^{-\beta \hat{H}}}+\frac{\mathrm{Tr}{\hat{O}}^q{e}^{-\beta \hat{H}}}{\mathrm{Tr}{e}^{-\beta \hat{H}}}$$

$$⟨{\hat{O}}^q⟩$$

$$⟨{\hat{O}}^q⟩\simeq \frac{1}{N_m}\sum _{\{q_j\}}{\hat{O}}^q(q)$$

$$⟨{\hat{O}}^c⟩=\frac{{\sum }_{\{q_j\}}{e}^{\beta h{\sum }_j{q}_j}{\mathrm{Tr}}_c{\hat{O}}^c(q)e^{-\beta \hat{H}(q)}}{{\sum }_{\{q_j\}}{e}^{-\beta F(q)}}$$

$$⟨{\hat{O}}^c⟩$$

$$⟨{\hat{O}}^c⟩\simeq \frac{1}{N_m}\sum _{\{q_j\}}⟨⟨{\hat{O}}^c(q)⟩⟩,$$

$$⟨\hat{H}⟩$$

$$C_{\mathrm{v}}$$

$$⟨⟨{\hat{H}}^2(q)⟩⟩$$

$$S_{qq}(Q)=\frac{1}{N_s^2}\sum _{ij}{e}^{iQ(R_i-R_j)}⟨q_i{q}_j⟩.$$

$$S_{qq}(Q)\simeq \frac{1}{N_m}\sum _{\{q_j\}}\left[\frac{1}{N_s^2}\sum _{ij}{e}^{iQ(R_i-R_j)}{q}_q{q}_j\right].$$

$${\chi }_{qq}(Q)$$

$$n_c=\frac{1}{N_s}\sum _j⟨{\hat{c}}_j^{†}{\hat{c}}_j⟩\simeq \frac{1}{N_m}\sum _{\{q_j\}}\left[\frac{1}{N_s}\sum _n{f}_F(E_n(q))\right].$$

$$N(\omega ,q)=\frac{1}{N_s}\sum _n\delta (\omega -E_n(q)),$$

Full text

I Introduction

Exactly solvable quantum many-body models play a fundamental role in understanding exotic quantum states in condensed matter physicsWen . In spatial dimension $d>1$, Kitaev’s toric-code, honeycomb model and certain $Z_{2}$ lattice gauge field models are prototypical examplesKitaev1 ; Kitaev2 ; Kogut ; Prosko , which shows novel quantum orders like Majorana quantum spin liquid, orthogonal metal, fractionalized Chern insulator, fracton order, Majorana superconductor and many-body localized stateKitaev2 ; Nandkishore ; Zhong1 ; Zhong2 ; Vijay ; Parameswaran ; Ng ; Smith . Their solvability is due to the existence of an infinite number of conserved quantities (for infinite lattice), resulting from intrinsic or emergent local $Z_{2}$ gauge symmetry. However, because of the designed local gauge symmetry, these models are far from the standard ones in condensed matter physics like Hubbard and Kondo lattice modelsHubbard ; Tsunetsugu . Thus, it is highly desirable to find solvable models without any fine-tuning local gauge structure.

In this paper, we introduce an exactly solvable model without any local gauge structure. It is a Kondo lattice model in its anisotropic limit, (also called the Ising-Kondo lattice model) which describes an itinerant electron interplaying with a localized $f$-electron moment via only longitudinal Kondo exchange Sikkema .

[TABLE]

$\hat{c}_{j\sigma}$ is the creation operator of conduction electron while $\hat{S}_{j}^{z}$ denotes localized moment of the $f$-electron at site $j$. $t_{ij}$ is hopping integral between $i,j$ sites and $J$ is the longitudinal Kondo coupling, which is usually antiferromagnetic, i.e. $J>0$. Chemical potential $\mu$ can be added to fix conduction electron’s density.

In literature, this model (with $x$-axis anisotropy) is proposed to account for the anomalously small staggered magnetization and large specific heat jump at hidden order transition in URu2Si2Sikkema ; Mydosh . (Due to lack of transverse Kondo coupling, enhancement of effective mass and related Kondo screening are irrelevant.) It can explain the easy-axis magnetic order and paramagnetic metal or bad metal behaviors in the global phase diagram of heavy fermion compoundsVojta ; Coleman2015 ; Si ; Coleman .

Importantly, in Eq. 1, we observe that $f$-electron’s spin/localized moment at each site is conservative since $[\hat{S}_{j}^{z},\hat{H}]=0$. Therefore, taking the eigenstates of spin $\hat{S}_{j}^{z}$ as bases, the Hamiltonian is automatically reduced to an effective free fermion model

[TABLE]

with $q$ emphasizing its $q$ dependence and $\hat{S}_{j}^{z}|q_{j}\rangle=\frac{q_{j}}{2}|q_{j}\rangle,q_{j}=\pm 1$. Now, the many-body eigenstate of the original model Eq. 1 can be constructed via the single-particle state of effective Hamiltonian Eq. 2 under given configuration of effective Ising spin $\{q_{j}\}$. So, Eq. 1 is exactly solvable and can be considered as an effective spinfull Falicov-Kimball (FK) modelFalicov . Because the procedure of reduction to the free fermion model involves only local conservation, the above model is solvable for arbitrary lattice geometry, spatial dimension and electron filling, in contrast to isotropic Kondo lattice modelAssaad , where notorious fermion minus-sign problem prevents exact solution/simulation. Moreover, including the external magnetic field along $z$-axis and Ising term $\sum_{ij}J_{ij}\hat{S}_{i}^{z}\hat{S}_{j}^{z}$ does not change the solvability.

To illustrate this, we consider the model (1) on the square lattice with nearest-neighbor-hopping $-t$. (See also Fig. 1(a)) In terms of analytical arguments and numerically exact lattice Monte Carlo (LMC) simulationCzajka , we have determined its ground state and its finite temperature phase diagram at half-filling in Fig. 1(b). There exist an antiferromagnetic insulator (AI), a Mott insulator (MI) and a correlated metal (CM). Both AI and MI have gapful single-particle excitations but CM are gapless metallic states. The transition from AI to CM or MI is a continuous Ising transition, and a smooth crossover appears between CM and MI with the opening of a gap at Fermi energy. The nature of these phases and transitions will be explored in the main text. For other bipartite lattice like honeycomb lattice, its phase diagram is similar to Fig. 1(b) with the Dirac semimetal replacing CM. When doping away from half-filling, we have established that various patterns of spin stripes occur at weak coupling and phase separation ultimately dominates at strong coupling. Since Ising antiferromagnetic order has been established in the heavy fermion compound CeCo(In1-xHgx)5Stock , our model could serve as a minimal model for this compound. Moreover, the spin stripe found here may be relevant to the electronic liquid crystal phenomena in a hidden order compound URu2Si2Okazaki . Finally, our model provides a benchmark for sophisticated numerical techniques like dynamic mean-field theory (DMFT)Kotliar and variational cluster approximation (VCA)Potthoff .

II Ground state

At zero temperature, when our system is half-filled ($\mu=0$) and on a bipartite lattice (here a square lattice), the ground-state configuration of $q_{j}$ has the two-fold degenerated checkerboard order $q_{j}=\pm(-1)^{j}$(also confirmed by our LMC simulation), according to the theorem proved by Kennedy and LiebKennedy . (For the generic lattice geometry, the ground-state configuration has to be searched by numerical techniques.) This can be seen as the Ising antiferromagnetic long-ranged order for a localized $f$-electron moment. For conduction electrons, the single-particle Hamiltonian in the ground state is thus

[TABLE]

which is just the familiar antiferromagnetic spin-density-wave (SDW) mean-field Hamiltonian with characteristic wave vector $\vec{Q}=(\pi,\pi)$ in Hubbard-like modelsKusko . But, we have to emphasize that this one is the exact ground state fermion Hamiltonian for $J>0$. Its quasi-particle dispersion is found to be $E_{k\sigma}=\pm\sqrt{\varepsilon_{k}^{2}+\frac{J^{2}}{16}}$, which splits the single free conduction electron band into two Hubbard-like bands with direct band gap $\Delta=\frac{J}{2}$. [$\varepsilon_{k}=-2t(\cos k_{x}+\cos k_{y})$ is a free conduction electron dispersion.] We conclude that the ground state of a half-filled model on the square lattice is an insulating antiferromagnetic state with Néel-like spin order.

III Finite-$T$ phase diagram

At finite temperature, one has to sum over all configurations of the effective Ising spin $\{q_{j}\}$, which can only be performed via Monte Carlo simulation. (See the Appendix for details.) We consider periodic $N_{s}=L\times L$ lattices with $L$ up to $16$. The resulting phase diagram is shown in Fig. 1(b). Here, AI is stable at low-$T$ since discrete Ising symmetry is able to break spontaneously into two dimensions $2D$ at finite $T$, and has a thermodynamic transition into CM in weak coupling or MI in strong coupling. The existence of AI is inferred by its checkerboard order parameter $\phi_{c}=\frac{1}{N_{s}}\sum_{j}(-1)^{j}\langle q_{j}\rangle$ and SDW structure factor $C_{\rm{SDW}}=\frac{1}{N_{s}^{2}}\sum_{j,k}(-1)^{j+k}4\langle\hat{s}_{j}^{z}\hat{s}_{k}^{z}\rangle$. ($\hat{s}_{j}^{z}=\frac{1}{2}(\hat{c}_{j\uparrow}^{{\dagger}}\hat{c}_{j\uparrow}-\hat{c}_{j\downarrow}^{{\dagger}}\hat{c}_{j\downarrow})$) In Fig. 2,

both $\phi_{c}$ and $C_{\rm{SDW}}$ saturate at low-$T$, which suggests the existence of AI. At high-$T$, $C_{\rm{SDW}}$ approaches zero and signals a transition to paramagnetic phases. This agrees with the divergence of specific heat $C_{\mathrm{v}}=\frac{\langle\hat{H}^{2}\rangle-\langle\hat{H}\rangle^{2}}{N_{s}T^{2}}$ and susceptibility $\chi_{q}=\frac{\langle S_{q}^{2}\rangle-\langle S_{q}\rangle^{2}}{T}$ ($S_{q}=\frac{1}{N_{s}^{2}}\sum_{j,k}(-1)^{j+k}q_{i}q_{j}$) at critical temperature $T_{c}$. We note that the maximal $T_{c}$ appears when Kondo coupling is comparable to the band-width of conduction electrons. ($J/t\sim 8$) The qualitative physics in AI can be understood by SDW mean-field theory, but such mean-field approximation underestimates thermal fluctuation and leads to unrealistic $T_{c}$.

III.1 Mott insulator and correlated metal

In the high-$T$ paramagnetic regime, based on the behavior of density of state (DOS) $N(\omega)$ for conduction electrons, magnetization under external magnetic field and conduction electrons’ distribution function $n_{c}(k)$, there exist MI and CM. (Fig. LABEL:fig:3 and 4)

MI appears at strong coupling and is a melting AI/SDW without long-ranged order but with (fluctuated) short-ranged order. Such short-ranged order gives rise to a single-particle gap observed in DOS (Fig. LABEL:fig:3(b)). Intuitively, conduction electrons mediate an antiferromagnetic Ising coupling $\sim\frac{J^{2}}{8t}$ between localized moments. At strong coupling, only temperature itself prevents the formation of long-ranged magnetic order, but short-ranged order survives. Thus, in short-time regime, conduction electrons feel the short-ranged order just as a true long-ranged order, and an excitation gap appears due to the effective molecular field applied to conduction electrons. It is important to note that the short-ranged order has classical Ising feature, thus it does not involve the resonance-valence-bond physics. So, MI here is irrelevant to quantum spin liquids but more like the featureless MI in Bose-Hubbard or FK modelsZhou ; Fisher ; Antipov . (We have checked that MI has very weak $T$-dependence on DOS, contrasting with radical reconstruction around coherent temperature in the usual Kondo insulatorJarrell .)

The spectral/single-particle feature of MI can be qualitatively captured by Hubbard-I approximationHubbard , where two Hubbard-like bands are observed and its DOS is similar to LMC’s results. In Hubbard-I approximation, the single-electron Green’s function is

[TABLE]

Here, the coherent factor $\alpha_{k}^{2}=\frac{1}{2}(1+\frac{\varepsilon_{k}}{\sqrt{\varepsilon_{k}^{2}+J^{2}/4}})$ and the quasiparticle dispersion is $\tilde{E}_{k}^{\pm}=\frac{1}{2}\left[\varepsilon_{k}\pm\sqrt{\varepsilon_{k}+J^{2}/4}\right]$. Although no antiferromagnetic order exists in MI, the strong local electron correlation (due to Kondo coupling) splits the band and drives the system into correlated MI. Since MI has a single-particle gap, its thermodynamics and transport properties are insulating and show exponential-$T$ behavior. ($C_{\mathrm{v}}$ in MI is extrapolated to vanishing at low-$T$, thus MI are not orthogonal metal-like exotic metalsNandkishore .)

At weak coupling, localized $f$-electron moment acts like uncorrelated random potential for conduction electrons, and the average over those disorder potentials leads to correlation correction for the latter one, which is CM. In contrast to gapped MI, CM has finite DOS around zero energy (Fermi energy) without SDW order.(Fig. LABEL:fig:3(a)) Thus, above $T_{c}$, the gapless CM has linear-$T$ behavior in specific heat ($C_{\mathrm{v}}\sim T$) and it shows Fermi liquid-like behavior. (Fig. 2(c))

When the external magnetic field $h$ is applied, in Fig. 4, we show the magnetization curvature of conduction electrons ($m_{c}=\sum_{j}\langle\hat{c}_{j\uparrow}^{{\dagger}}\hat{c}_{j\uparrow}-\hat{c}_{j\downarrow}^{{\dagger}}\hat{c}_{j\downarrow}\rangle$), localized $f$-electron moment ($m_{f}=2\sum_{j}\langle S_{j}^{z}\rangle$) and the total one ($m_{z}=m_{c}+m_{f}$). Under large $h$, all states evolve into a fully polarized state. At small $h$, CM shows characteristic metallic linear-$h$ behavior, [Fig. 4(a) Pauli-like susceptibility] and both MI and AI have a spin gap. [Figs. 4(b) and 4(c)] A strong first-order transition from AI to a fully polarized state is observed, indicating the absence of a field-induced magnetic quantum phase transition. The quasiparticle behavior in CM and MI can also be inspected in the conduction electrons’ distribution function $n_{c}(k)=\langle\hat{c}_{k\sigma}^{{\dagger}}\hat{c}_{k\sigma}\rangle$. [Fig. 4(d)] In CM, there is a clear jump in $n_{c}(k)$ around (underlying) Fermi surface [$k=(\pi/2,\pi/2)$] while only smooth evolution exists in MI.

There is a smooth crossover between CM and MI as inspecting the evolution of energy density $E=\frac{1}{N_{s}}\langle\hat{H}\rangle$ and double occupation $d_{c}=\frac{1}{N_{s}}\sum_{j}\langle c_{j\uparrow}^{{\dagger}}c_{j\uparrow}c_{j\downarrow}^{{\dagger}}c_{j\downarrow}\rangle$ versus $J$, where $E$ and $d_{c}$ have nearly linear dependence on $J$ and no singularity has been observed(Fig. LABEL:fig:4).

The transition from MI and CM to AI is continuous as seen from order parameter $\phi_{c}$ and $C_{\mathrm{SDW}}$ around $T_{c}$ [Fig. 2(a) and 2(b)]. The corresponding histogram for energy density $E$ only shows one peak structure around $T_{c}$ (not shown here), which excludes the strong first-order transition, however the possibility of a very weak discontinuous transition may be still possible. Assuming a continuous phase transition, the finite-size scaling analysis of checkerboard order parameter $\phi_{c}$ in Fig. LABEL:fig:6 suggests that critical behavior of our model belongs to a $2D$ Ising universality class with an order parameter critical exponent $\beta=1/8$ and correlation length critical exponent $\nu=1$. (The standard scaling form $\phi_{c}L^{\beta/\nu}=g[(T-T_{c})L^{1/\nu}]$ is used to extract $\beta$ and $\nu$Binder .)

III.2 Effective field theory analysis

We have provided data from LMC simulation, which explores the whole phase diagram and the magnetic transition. In this subsection, we use an effective field theory analysis to understand these phases and related phase transition.

In terms of path integral formalism, we have the following action:

[TABLE]

where $\bar{c}_{j\sigma},c_{j\sigma}$ are the Grassmann field for conduction electrons, $\phi_{j}$ denotes the localized moment, and $\lambda_{j}$ is the (dynamic) Lagrangian multiplier. Since above $T_{c}$, the localized moment is disordered, it is reasonable to assume $i\lambda_{j}=m^{2}$ with $m$ being the mass of $\phi_{j}$. This is equivalent to the saddle point approximation for constraint $\phi_{j}^{2}=1$. Then, integrating over $\phi_{j}$, the action only including conduction electrons reads

[TABLE]

In the return to Hamiltonian formalism, such action gives rise to a symmetric Hubbard model with effective Hubbard interaction $U_{eff}=\frac{J^{2}}{8m^{2}}$. Therefore, it seems that the high-$T$ paramagnetic states could be related to solutions of half-filled symmetric Hubbard model. However, such relation is not strict due to the saddle point approximation. This explains why at finite-$T$ there is first-order transition between metallic states and MI in Hubbard-like models while in our case, only a smooth crossover appearsKotliar . Moreover, when approaching $T_{c}$, the fluctuation effect cannot be neglected, thus it is not reliable to treat our system as an antiferromagnetic Heisenberg model as if $U_{eff}\rightarrow\infty$ for $m^{2}\rightarrow 0$. Such feature is consistent with different (spin) symmetry in our model ($Z_{2}$) and Heisenberg model ($SU(2)$).

Alternatively, if we integrate out conduction electrons and set $\phi_{j}=(-1)^{j}\phi_{j}$ to emphasize the dominating antiferromagnetic correlation, the resulting effective theory is the celebrated Hertz-Millis-Moriya theory with $\phi$ being antiferromagnetic Ising order parameter fieldHertz ; Millis ; Moriya . Because of the nesting of the Fermi surface on square lattice, antiferromagnetic quantum criticality is avoided and only thermal critical behaviors preserve, (Landau damping due to particle-hole excitation of conduction electrons is subleading and can be neglected away from quantum critical regime) thus the effective theory is replaced by the classic $\phi^{4}$ theory and it corresponds to $2D$ Ising universality class. This agrees with our LMC simulation.

IV Doping away from half-filling

Doping the half-filled system leads to inhomogeneous magnetic order or spin stripe order, as seen in Fig. 7. Here, the ground state configuration of the localized $f$-electron moment ($\langle\hat{S}_{j}^{z}\rangle$) versus chemical potential $\mu$ is shown, and various patterns of spin stripes occur to minimize the free energy. Conduction electrons have similar spin stripe structure but with opposite direction and small amplitude. When Kondo coupling $J$ is too large compared with the the band-width of conduction electrons, phase separation ultimately dominatesCzajka2 . These findings are qualitatively consistent with spiral/stripe magnetic phases in the isotropic Kondo lattice modelCosta ; Peters .

V Conclusion

In conclusion, we have provided a solvable lattice fermion Hamiltonian, whose solvability is due to the conserved localized moment at each lattice site. The case study on the square lattice has established the existence of Ising antiferromagnetic order, spin stripe, Fermi liquid-like metallic state and Mott insulator. Importantly, Ising antiferromagnetism and spin stripe orders found here may be relevant to the uniaxial magnetic order in CeCo(In1-xHgx)5 and the quantum liquid crystal order in URu2Si2. If frustration is introducedCzajka , an exotic metallic phase such as orthogonal metal can be exploredNandkishore ; Hohenadler . Topological states of matter like topological SDW is also an interesting issue when spin-orbit coupling is addedZhong2013 . These will be left for our future work.

Acknowledgments

We thank Hantao Lu, Yuehua Su and Yu Liu for helpful discussion on the phase diagram and nature of each phase. This research was supported in part by NSFC under Grant No. $11704166$, No. $11834005$, No. $11874188$ and the Fundamental Research Funds for the Central Universities.

Y.Z. proposed the project, Y.Z. and W.Y. carried out the calculation.

Appendix A Details of Lattice Monte Carlo simulation

When exploring the finite temperature properties, one has to sum over all configurations of effective Ising spin $\{q_{j}\}$, which can only be performed via Monte Carlo simulationCzajka .

A.1 Monte Carlo simulation

Firstly, we know that the equilibrium state thermodynamics is controlled by the partition function

[TABLE]

Here, the trace is split into $c$-fermion and $\hat{S}^{z}$, where the latter is transformed into the summation over all possible configuration $\{q_{j}\}$. For each single-particle Hamiltonian $\hat{H}(q)$, it can be easily diagonalized into

[TABLE]

where $E_{n}$ is the single-particle energy level and $\hat{d}_{n}$ is the quasi-particle. The fermion $\hat{d}_{n}$ is related into $\hat{c}_{j}$ via

[TABLE]

with $\phi_{n}^{j}\equiv\langle j|n\rangle$ while the creation operator $\hat{c}_{j}^{{\dagger}}$ has

[TABLE]

Now, the trace over $c$-fermion can be obtained as

[TABLE]

This is the familiar result for free fermion, however one should keep in mind that $E_{n}$ actually depends on the effective Ising spin configuration $\{q_{j}\}$, thus we write $E_{n}(q)$ to emphasize this fact.

So, the partition function reads

[TABLE]

where we have defined an effective free energy

[TABLE]

In this situation, we can explain $e^{-\beta F(q)}$ or $\rho(q)=\frac{1}{\mathcal{Z}}e^{-\beta F(q)}$ as an effective Boltzmann weight for each configuration of $\{q_{j}\}$ and this can be used to perform Monte Carlo simulation just like the classic Ising model. Specifically, we start with random or chosen configuration of $\{q_{j}\}$, then try to flip each $q_{j}$ to $-q_{j}$. The relative probability for such flip is determined by effective Boltzmann weight

[TABLE]

where $F_{old}(q),F_{new}(q)$ denotes the effective free energy before and after flip. This is the weight used in the classic Metropolic importance sampling. Alternatively, one can use the so-called bath algorithm, which means

[TABLE]

such that the probability is seemingly normalized. Then, we generate a random number $a$ from uniform distribution $[0,1]$ and compare this with $\tilde{r}$. If $\tilde{r}>a$, the flip of $q_{j}$ is accepted and $F_{old}$ is updated into $F_{new}$, otherwise we reset $q_{j}$ to its original value before the flip. By trying to flip $q_{j}$ over all sites, this is called a sweep and doing such sweep many times, the system can be equilibrium and the next sweeps are used to calculate physical observable.

A.2 On physical quantities

To calculate physical quantities, we consider generic operator $\hat{O}$, which can be split into part with only Ising spin $\{q_{j}\}$ and another part with fermions,

[TABLE]

Then, its expectation value in the equilibrium ensemble reads

[TABLE]

For $\hat{O}^{q}$, we have

[TABLE]

In the Metropolic importance sampling algorithm, the above equation means we can use the simple average to estimate the expectation value like

[TABLE]

where $N_{m}$ is the number of sampling and the sum is over each configuration. $\hat{O}^{q}(q)$ is a number since we always work on the basis of $\{q_{j}\}$.

For $\hat{O}^{c}$,

[TABLE]

and we can insert $\frac{e^{-\beta F(q)}}{e^{-\beta F(q)}}$ in the numerator, which leads to

[TABLE]

This means

[TABLE]

where $\langle\langle\hat{O}^{c}(q)\rangle\rangle=\frac{\mathrm{Tr}_{c}\hat{O}^{c}(q)e^{-\beta\hat{H}(q)}}{\mathrm{Tr}_{c}e^{-\beta\hat{H}(q)}}$ is calculated based on the Hamiltonian $\hat{H}(q)$. More practically, such statement means if fermions are involved, one can just calculate with $\hat{H}(q)$. Then, average over all sampled configuration gives rise to the wanted results.

A.3 Energy, correlation and others

The total energy of our system is an essential quantity and can be calculated as

[TABLE]

Here, the $\langle\langle\hat{H}(q)\rangle\rangle$ is just the summation over all quasi-particles for given configuration. Now, the specific heat $C_{\mathrm{v}}$ can be found by $C_{\mathrm{v}}=\frac{\langle\hat{H}^{2}\rangle-\langle\hat{H}\rangle^{2}}{T^{2}}$, which means

[TABLE]

where

[TABLE]

Next, we consider the spin-spin correlation function $S_{qq}(Q)$, which is defined as

[TABLE]

This object is designed to detect the static or dynamic (fluctuation) order with the characteristic wave-vector $Q$. If spin orders in certain $Q$, the value of the corresponding $S_{qq}(Q)$ should reach $\mathcal{O}(1)$. Now, using Eq. 12, we obtain

[TABLE]

Usually, we may also use its susceptibility

[TABLE]

to locate the position of long-ranged order. Obviously, if orders appear, $\chi_{qq}(Q)$ has to diverge for certain $Q$ at some critical temperature $T_{c}$. At $T_{c}$, the specific heat $C_{\mathrm{v}}$ diverges as well.

For total $c$-fermion density, in terms of Eq. 13, it is easy to show that

[TABLE]

Since we work on grand canonical ensemble, if chemical potential is setting to zero, the model is symmetric and $c$-fermion should be half-filled. So, in this situation, $n_{c}$ is $1$, irrespective of temperature.

If we are interested in the density of state $N(\omega)$ of $c$-fermion, we know that for $\hat{H}(q)$, it is given by

[TABLE]

thus,

[TABLE]

Finally, when we calculate fermion’s correlation function like $\langle\hat{c}_{i}^{{\dagger}}\hat{c}_{j}\hat{c}_{k}^{{\dagger}}\hat{c}_{l}\rangle$,

[TABLE]

Then, using the Wick theorem for these free fermions, we get

[TABLE]

Next, for each one-body correlation function like $\langle\langle\hat{c}_{i}^{{\dagger}}\hat{c}_{j}\rangle\rangle$, one can transform these objects into their quasiparticle basis,

[TABLE]

Similarly, we have

[TABLE]

Then, combining Eq. 15 and inserting this expression into the Eq. 14, we can calculate the fermion’s correlation function.

For instance, the spin-density-wave structure factor is calculated as

[TABLE]

Appendix B Effect of external magnetic field: Magnetization

If $z$-axis external magnetic field is included, our model is changed to

[TABLE]

With eigenstates of spin $S_{j}^{z}$ as bases, the corresponding effective free fermion model reads

[TABLE]

which is readily to be simulated by LMC with effective free energy $F(q)=-T\sum_{n}\ln(1+e^{-\beta E_{n}(q)})-h/2\sum_{j}q_{j}$. In the main text, by using above formalism, we have calculated the magnetization of conduction electron and localized moment under external field $h$.

Appendix C Phase separation

When doping away from half-filling, the microscopic phase separation is a leading effect for intermediate and large Kondo coupling. An example is given in Fig. 8, where the phase separation is clearly seen with different ferromagnetic droplets.

In addition, we note that if conduction electron’s density is very small ($n_{c}\sim 0.1$), ferromagnetism is stable both in weak and strong coupling regime.

Appendix D Mean-field approximation for SDW state

The SDW mean-field Hamiltonian can be obtained by decoupling the Kondo coupling term via

[TABLE]

where the magnetic order parameters are defined by

[TABLE]

thus, the mean-field Hamiltonian reads

[TABLE]

Here, the summation over momentum $k$ is restricted into the folded first Brillouin zone (for square lattice, the area is interior surrounded by $k_{y}=-k_{x}\pm\pi,k_{y}=k_{x}\pm\pi$) Therefore, its free energy can be found as

[TABLE]

with $E_{k\pm}=\pm\sqrt{\varepsilon_{k}^{2}+\frac{J^{2}m_{f}^{2}}{16}}$. So, the mean-field self-consistent equations can be derived by $\frac{\partial F}{\partial m_{c}}=0,\frac{\partial F}{\partial m_{f}}=0$,

[TABLE]

[TABLE]

From these two equations, one can obtain the equation for SDW critical temperature $T_{c}^{mf}$ in the mean-field approximation.

In Fig. LABEL:fig:5(a), we show the SDW critical temperature $T_{c}^{mf}$ versus $J$, which has linear dependence on $J$. However, since no thermal fluctuation is included, the critical temperature $T_{c}^{mf}$ shows a wrong tendency compared to the exact critical temperature $T_{c}$ from LMC.

Appendix E Hubbard-I approximation for MI

In the paramagnetic MI, one can follow the Hubbard-I approximation to give a rough solution for MI itselfHubbard . To this purpose, we use equation of motion (EOM) formalism and define the retarded Green’s function as

[TABLE]

By using the standard EOM relation,

[TABLE]

for our model, it follows that

[TABLE]

Here, $\Delta_{im}$ denotes $m$ is the nearest-neighbor site of $i$. For $\langle\langle\hat{S}_{i}^{z}\hat{c}_{i\sigma}|\hat{c}_{j\sigma}^{{\dagger}}\rangle\rangle_{\omega}$, we have

[TABLE]

If no further EOM are involved, to close EOM, we have to decouple $\langle\langle\hat{S}_{i}^{z}\hat{c}_{l\sigma}|\hat{c}_{j\sigma}^{{\dagger}}\rangle\rangle_{\omega}$ as

[TABLE]

In paramagnetic strong coupling regime, there is no magnetic order, thus $\langle\hat{S}_{i}^{z}\rangle=0$. Meanwhile, the contribution from $\langle\langle\hat{S}_{i}^{z}\hat{c}_{l\sigma}|\hat{c}_{j\sigma}^{{\dagger}}\rangle\rangle_{\omega}$ vanishes due to the decoupling and the above equation has a complete solution:

[TABLE]

which can be written as

[TABLE]

Now, performing the Fourier transformation

[TABLE]

we have

[TABLE]

Here,$-t\Delta_{im}e^{ik(R_{m}-R_{i})}=-t\sum_{\delta}e^{-ik\delta}=\varepsilon_{k}$ and we have

[TABLE]

Since $R_{i},R_{j}$ is arbitrary, the above equation is valid if $\left(\omega-\frac{J^{2}}{16\omega}-\varepsilon_{k}\right)G_{\sigma}(k,\omega)-1=0$, which gives the single-particle Green’s function as

[TABLE]

Here, the coherent factor $\alpha_{k}^{2}=\frac{1}{2}\left(1+\frac{\varepsilon_{k}}{\sqrt{\varepsilon_{k}^{2}+J^{2}/4}}\right)$ and

[TABLE]

Therefore, although no SDW state appears, the band of conduction electron has been split into two Hubbard-like bands. This is driven by the local Kondo interaction and the resulting state is MI if half-filling is satisfied. Moreover, the DOS in the present Hubbard-I approximation can be obtained by

[TABLE]

where $\Gamma$ is the damping factor and we give an example in Fig. LABEL:fig:5(b) where the DOS of conduction electron is the constant DOS.