Quasiparticle Band Structure and Spin Excitation Spectrum of the Kondo Lattice
R. Eder
 


Karlsruhe Institute of Technology,
Institut für Festkörperphysik, 76021 Karlsruhe, Germany


A formulation of the Kondo lattice Hamiltonian in terms of bond particles
is derived and solved in two different approximations. The bond particles 
correspond to the eigenstates of a single unit cell
and are bosons for states with even electron number and fermions
for states with an odd electron number. As a check various 
physical quantities
are calculated for the 1D Kondo insulator and good agreement with
numerical results is obtained for J/t>1.
 
71.27.+a,75.30.Mb,71.28.+d
Quasiparticle Band Structure and Spin Excitation Spectrum of the 
Kondo Lattice
    R. Eder
    December 30, 2023
===============================================================================




§ INTRODUCTION


Metallic compounds containing Cerium, Ytterbium or Uranium - the so-called 
Heavy  fermions - continue to be a much studied field of solid state physics. 
These materials show a number of remarkable phenomena which are widely believed
to be caused by the strong Coulomb repulsion between the electrons in the 
4f-shells of Cerium and Ytterbium or the 5f-shell of Uranium. A long known 
phenomenon is the crossover 
from a lattice of localized f electrons coexisting with weakly or moderately
correlated conduction electron bands at high temperature, to an exotic Fermi 
liquid with strongly correlation-enhanced effective masses at low temperature, 
whereby the f electrons now contribute to the Fermi-surface 
volume<cit.>.
The low temperature Fermi-liquid phase can undergo magnetic ordering 
transitions whereby the transition temperature often can be 
tuned to zero by external parameters
such as temperature, pressure, magnetic field or alloying, resulting in 
quantum critical points and the ensuing non-Fermi-liquid behavior
and superconducting domes<cit.>. 
An intriguing feature thereby is the fact that whereas the f electrons do 
contribute to the Fermi surface volume in the paramagnetic phase they seem
to `drop out' of the Fermi-surface volume at some of these transitions.

Heavy fermion compounds can be described by the
Kondo-lattice model (KLM). In the simplest case of no orbital degeneracy
each unit cell n contains one f orbital and one conduction band orbital,
and denoting the creation operators for electrons in these orbitals by
f_n,σ^† and c_n,σ^† the Kondo lattice Hamiltonian 
is H=H_t + H_J with

    H_t   =   ∑_m,n ∑_σ  t_m,n c_m,σ^† c_n,σ^,
    
    
    H_J   =   ∑_n h_n,
    
    h_n   =    J S_n,f·S_n,c,
    S_n,c   =   1/2 ∑_σ,σ'   c_n,σ^† τ_σ,σ' 
     c_n,σ'^.

Here τ denotes the vector of Pauli matrices 
(an analogous definition holds for S_n,f).
An important feature of the KLM is the constraint to have precisely one 
electron per f orbital:

    ∑_σ f_n,σ^† f_n,σ^   =    1,

which must hold separately for each n. In the following
we consider a lattice with N unit cells and N_c conduction electrons.
The KLM can be derived from the more realistic periodic
Anderson model (PAM) by means of the Schrieffer-Wolff 
transformation<cit.>.

The impurity versions of the Kondo and Anderson model are well 
understood. Approximate solutions can be obtained by
variational wave functions<cit.>,
mean-field (or saddle point) approximation to the
exchange term<cit.>,
or Green's function techniques where the hybridization or exchange between 
the f level and conduction band are treated as 
perturbation<cit.>.
Thereby the inverse of the degeneracy n_f of the f level, 1/n_f, 
often plays the role of a small parameter<cit.>. Exact solutions of the
impurity models can be obtained by renormalization group<cit.> 
and Bethe ansatz<cit.>. 
The ground state of the impurity  is a singlet formed from
the f electron on the impurity and an extended screening state
formed from the conduction band states, whereby for weak coupling (J≪ t) 
the binding energy - the so-called Kondo temperature - is 
k_BT_K∝ We^-1/ρ J. Here W and ρ are
bandwidth and density of states of the conduction band.

The lattice versions of the model are less-well understood.
A noteworthy result is the fact that even for the KLM,
where the f electrons are strictly localized, they do contribute
to the Fermi surface volume<cit.> provided the
system is a Fermi liquid. Approximate results for the KLM 
and PAM have been obtained in the mean-field (or saddle-point)
approximation. For the KLM the exchange term, which is quartic in electron 
operators is 
factorized<cit.>, whereas for the PAM a 
slave-boson representation is 
used<cit.>. The models also have been studied
using Gutzwiller-type trial wave functions<cit.>.
The resulting band structure is consistent
with a simple hybridization picture: a dispersionless effective f band
close to the Fermi energy of the decoupled conduction band
hybridizes with the conduction electron band via an effective
hybridization matrix element ∝ k_B T_K at weak coupling.
This results in a Fermi surface with a volume corresponding to
itinerant f electrons and the `heavy bands' characteristic of
heavy fermion compounds. 
In addition to mean-field theories, a large amount of quantitative
results has also been gathered by numerical methods such as 
Density Matrix Renormalization Group (DMRG)<cit.>,
Dynamical Mean-Field (DMFT) calculations<cit.>,
Quantum Monte-Carlo (QMC)<cit.>, 
series expansion (SE)<cit.>,
Variational Monte-Carlo (VMC)<cit.>,
exact diagonalization<cit.>,
Dynamical Cluster Approximation<cit.>
and Variational Cluster Approximation<cit.>.
For the paramagnetic phase the numerical techniques produce band 
structures which are consistent with the hybridization picture, 
whereby it has to be kept in mind that numerical methods often have
problems to access the limit of small J/t and thus to reproduce the
Kondo scale k_BT_K. However, the heavy quasiparticles and the fact that
the f electrons do participate to the Fermi surface in the KLM and PAM
are reproduced.

Considerable effort was devoted to a study of the magnetic phase 
transitions which are believed to be due to a competition<cit.> 
between the singlet formation in the impurity model and the
RKKY-interaction between the f spins<cit.> mediated by the conduction 
electrons. A controversial question is whether
the heavy quasiparticles persist at the magnetic transition, so that
this may be viewed as the heavy Fermi liquid undergoing a conventional
spin-density-wave transition, or whether the magnetic ordering
suppresses the heavy Fermi liquid alltogether. Numerous studies have addressed 
magnetic ordering<cit.>
but open questions remain.

It is the purpose of the present manuscript to present a theory for
the single-particle band structure and spin excitation spectrum of the
Kondo lattice which relies on the interpretation of the eigenstates of
a single cell as fermionic or bosonic particles, which we call bond particles.
Bond particle theory was proposed originally by 
Sachdev and Bhatt<cit.> to study
spin systems and applied to spin ladders<cit.>,
bilayers<cit.>, intrisically dimerized spin 
systems<cit.> and the `Kondo necklace'<cit.>.
It was also applied to the PAM<cit.>, as well as 
antiferromagnetic (AF) ordering in the planar 
KLM<cit.>, 
a discussion of its different AF phases<cit.> and the band structure
in the AF phase<cit.>.



It is by nature a strong-coupling theory which should work best in the
(unphysical) limit of J/t≫ 1. However, as will be shown below by comparison
to numerical results, there is some reason to hope that the theory retains its
validity down to t/J≈ 1 which may be sufficient to discuss
magnetic ordering phenomena.



§ HAMILTONIAN


We consider eigenstates of the single cell exchange term
h_n. Introducing the matrix 4-vector
γ=(τ_0,τ) with τ_0=1
the state 4-vector β_n=(s_n, t_n) 
is<cit.>

    |β_n⟩ =1/√(2) ∑_σ,σ' 
    c_n,σ^† (γiτ_y)_σ,σ' f_n,σ'^†
    |0⟩.

These are the
singlet (|s_n⟩) with energy
-3J/4 and the three components of the triplet
(|t_n⟩) with energy J/4.

The single-cell states with an odd number of electrons 
(which have energy 0) are 

    |a,n,σ⟩   =    f_n,σ^†|0⟩ ,
    
    |b,n,σ⟩   =    c_n,↑^† c_n,↓^†
    f_n,σ^†|0⟩.

We now rewrite the KLM as a Hamiltonian for
bosons and fermions which correspond to these single-cell eigenstates.
More precisely, if a given cell n is in one of the states
(<ref>) with 2 electrons we consider it as occupied by a boson, 
created by the respective operator 4-vector (s_n^†,t_n^†) 
whereas if it is in one of the states (<ref>) with a single (three) 
electrons we consider it as occupied by a fermion created by 
a_n,σ^† (b_n,σ^†). The latter correspond
to the `bachelor spins' in the U/t=∞ Hubbard model
to which the KLM reduces<cit.> for J/t→∞  and 
N_c  N. 

In order for this representation to make sense each cell must be occupied
by precisely one of these particles resulting in the constraint
(to be obeyed for each n)

    s_n^† s_n^ + ∑_σ( a_n,σ^† a_n,σ^ +
    b_n,σ^† b_n,σ^)+ t_n^†·t_n^ = 1.

On the other hand, each of the basis states (<ref>) and
(<ref>) obeys the constraint (<ref>) exactly,
so that this constraint is `built in' into the theory.
In terms of the bond particles the exchange term is

    H_J   =   ∑_n ( 3J/4∑_σ 
    (b_n,σ^† b_n,σ^ + a_n,σ^† a_n,σ^)
    + J  t_n^†·t_n^ )
                                          
     - 3NJ/4.

On the other hand, we might also write

    H_J   =   ∑_n (  -3J/4 s_n^† s_n^
    + J/4 t_n^†·t_n^ ).

As long as the constraint (<ref>) holds these two
forms are equivalent - we wil continue to use (<ref>).
From now on we take a fermion operator with omitted
spin index to denote a two-component column vector, e.g.

    c_n   =   ([ c_n,↑; c_n,↓ ]).

In this notation the representation of the electron annihilation
operator in terms of the bond particles is

    c_n=1/√(2) :( (s_n+  t_n·τ) 
    iτ_y a_n^† - (s_n^† -  t_n^†·τ) 
     b_n ):

where :…: denotes normal ordering.
Up to the numerical prefactors the form of this equation
follows from the requirement that both sides be covariant
spinors and the fact that t and t^† 
are vector operators. 


The representation of the kinetic energy H_t is obtained by 
substituting (<ref>) into (<ref>). One obtains
H_t=H_1+H_2+H_3+H_4 with






    H_1   =   ∑_m,n t_m,n/2 ∑_σ ( b_m,σ^† b_n,σ^ -
    a_m,σ^† a_n,σ^ ) s_n^†  s_m^ 
    -∑_m,n t_m,n/2 [ ( 
    b_m,↑^† a_n,↓^† - b_m,↓^† a_n,↑^† ) s_m^  s_n^ +H.c. ],
    
    H_2   =   ∑_m,n t_m,n/2 ∑_σ ( b_m,σ^† b_n,σ^ -
    a_m,σ^† a_n,σ^ ) t_n^†·t_m^ 
    +∑_m,n t_m,n/2 [ ( 
    b_m,↑^† a_n,↓^†
    - b_m,↓^† a_n,↑^† ) t_m^·t_n^ + H.c ],
    
    H_3   =    - ∑_m,n t_m,n/2 [ ( π^†_m,n·( s_m^ t_n^ - t_m^ s_n^) + H.c. )
    +( b_m,n - a_m,n )·
    (t_n^† s_m^ + s_n^† t_m^ ) ],
    
    H_4   =   ∑_m,n t_m,n/2 [ [ iπ^†_m,n·
    ( t_m^×t_n^ )
    + H. c.  ]
    - i ( b_m,n - a_m,n )·
    ( t_n^†×t_m^ ) ],


with the following vectors formed from the fermions:

    b_m,n   =   ∑_σ,σ' b_m,σ^† τ_σ,σ' b_n,σ'^,
    
    a_m,n   =   ∑_σ,σ' a_m,σ^† τ_σ,σ' a_n,σ'^,
    
    π^†_m,n   =   ∑_σ,σ' b_m,σ^† (τ iτ_y)_σ,σ' a_n,σ'^†.

Strictly speaking the individual terms in this Hamiltonian
have to be `site-wise normal ordered' e.g.
b_m,σ^† b_n,σ^ t_n,α^† t_m,α^→ b_m,σ^† t_m,α^ t_n,α^† b_n,σ^
but since this normal ordering always involves commutation of
a fermion and a boson neither nonvanishing commutators nor
Fermi signs will arise.
The number of electrons - including the localized
f electrons - in the system is

    N_e   =   2 ∑_n ( s_n^† s_n^ + t_n^†·t_n^ )
    
                          
    +∑_n,σ( 3 b_n,σ^† b_n,σ^ + a_n,σ^† 
    a_n,σ^)
       =   2N + ∑_n ∑_σ( b_n,σ^† b_n,σ^ - 
    a_n,σ^† a_n,σ^)
       =   ∑_n ∑_σ( b_n,σ^† b_n,σ^ + 
    a_n,σ^ a_n,σ^† ),

where (<ref>) was used to obtain the second line.
The Hamiltonian (<ref>)+(<ref>) together with the 
constraint (<ref>) provides an exact representation of the
KLM. On the other hand, it
is compliated and impossible to solve even approximately e.g. by diagrammatic
methods, due to the constraint (<ref>),
which is equivalent to an infinitely strong Hubbard-like
repulsion between the bond particles. The whole formulation
in terms of bond particles will be useful only if we can identify
the fermions and triplet bosons as approximate quasiparticles and
spin excitations of the system and find a way 
to extract a sufficiently simple yet accurate theory for these.
A considerable simplification becomes possible by making use of
the fact - to be verified below - that over large regions of parameter space 
the densities of fermions and triplets are relatively small so that the vast
majority of cells is in the singlet state. Thus, if one can get
rid of the singlets by either considering them as
condensed or by re-interpreting the singlet as the 
`true vacuum state' of a cell, one retains a theory for a system of
fermions and bosons which in principle are still subject to the infinitely 
strong repulsion implied by the constraint (<ref>) but which have
a low density so that relaxing the constraint may be a reasonable approximation.
Put another way, by using the bond particles one can trade
the constraint (<ref>) which refers to a dense system of electrons - the
density of f electrons is 1/cell - for a constraint like
(<ref>) without singlets
which refers to a system of particles with a relatively low density.
In the following, we explore two possible approximation 
schemes to  `get rid of the singlets' and compare the results to numerical 
calculations. It should also be noted that while we will not
do so in the following, it is in principle
possible to deal with strong repulsion in a low density
system by well-known field theoretical methods<cit.>. For
the case of bond bosons in spin systems this has in fact been carried
out explicitely by Kotov et al.<cit.> and 
Shevchenko et al.<cit.>.
We will compare the results from bond particle theory to numerical results
for the paramagnetic state in a 
one-dimensional chain with only nearest neighbor hopping -t and
n_e=2 - i.e. the 1D Kondo insulator<cit.>,
throughout t is the unit of energy.



§ MEAN FIELD THEORY


As a first approximation we study the Hamiltonian in mean-field
approximation. This approximation was applied previously to spin 
systems<cit.> and to antiferromagnetic ordering in the 
planar KLM<cit.>.
Since we are interested in the paramagnetic phase
we initially drop the terms H_3 and H_4.
H_3 describes  pair creation and propagation processes
whereby a single triplet-boson is absorbed or emitted. In mean-field
theory this term would contribute only in a state where the triplets
are condensed<cit.> i.e. a magnetically ordered 
state<cit.>.
Similarly, H_4 describes pair creation and propagation processes
whereby two triplets coupled to a vector are emitted/absorbed. The resulting 
vector-like order parameters would be important to describe a state with
incommensurate or spiral magnetic order but vanish in a rotationally invariant 
state.

In the remaining terms H_J+H_1+H_2 the singlets are assumed to be condensed 
whence the corresponding operators can be replaced by a real number,
s_n^†,s_n^→ s. The condensation amplitude
s now is a freely variable internal parameter of the system, to be 
determined by minimization of the Helmholtz free energy.
The constraint  (<ref>) is replaced by the global
constraint

    ∑_n ( s^2 + t_n^†·t_n^ +
    ∑_σ (b_n,σ^† b_n,σ^ + a_n,σ^† a_n,σ^)
     )   =   N.

We perform a Hartree-Fock factorization of the
quartic terms in H_2 and add the constraints
(<ref>) and (<ref>)
using the Lagrange multipliers λ and μ, respectively.
We call the resulting Hamiltonian H_MF and have
H_MF = H_F+ H_B + Nc. The fermionic Hamiltonian is

    H_F   =     e_0 ∑_n,σ ( b_n,σ^† b_n,σ^ +
    a_n,σ^† a_n,σ^ )
            +∑_m,n t̃_m,n ∑_σ ( b_m,σ^† 
    b_n,σ^ -
    a_m,σ^† a_n,σ^ )
              -∑_m,n [Δ_m,n ( 
    b_m,↑^† a_n,↓^†
    - b_m,↓^† a_n,↑^† ) + H.c.]
              -μ ∑_n,σ ( b_n,σ^† b_n,σ^ +
    a_n,σ^ a_n,σ^† )

where e_0=3J/4-λ and

    t̃_m,n   =   t_m,n/2 s^2 +ζ_m,n
    Δ_m,n   =   t_m,n/2 s^2 + η_m,n.

with
 ζ_m,n=t_m,n/2 ⟨t_n^†·t_m^ ⟩,
η_m,n=-t_m,n/2 ⟨t_m^·t_n^ ⟩.
We consider a translationally invariant and isotropic state and accordingly
assume that expectation values such as ζ_m,n depend only on
|R_m - R_n| so that

    H_F    =    e_0 ∑_k,σ ( b_k,σ^† 
    b_k,σ^ -a_-k,σ̅^ a_-k,σ̅^† )
          + ∑_k,σ (t̃_k -μ)
    (  b_k,σ^† b_k,σ^ +
    
    a_-k,σ̅^ a_-k,σ̅^† ) 
                  -∑_k,σ sign(σ)(Δ_k 
    b_k,σ^† 
    a_-k,σ̅^† + H.c.)
                           - 2∑_k (t̃_k -e_0),

with the Fourier transform of the hopping integral

    t̃_k   =   ∑_r e^i k·r t̃_r
       =   ∑_α z_αt̃_αγ_α( k)

and an analogous definition of Δ_k.
Here α denotes shells of symmetry-equivalent
neighbors of a given site, z_α the number of neighbors
belonging to a shell, and γ_α the respective tight-binding
harmonic.
This can be solved by the unitary transformation










    α_k^†   =      
    u_k b_k^† + v_k 
    iτ_y a_-k^,
    β_k^†   =   
    -v_k b_k^† + u_k 
    iτ_y a_-k^,

so that

    H_F   =   ∑_k,σ ( E_α,k α_k,σ^†α_k,σ^
    + E_β,k β_k,σ^†β_k,σ^ )
                                           
    - 2∑_k (t̃_k -e_0).

Here E_ν,k=t̃_k± W_k - μ
and α corresponds to the lower of the two energies.
Thereby

    W_k   =   √( e_0^2 + Δ_k^2),
    
    u_k   =   -√(W_k-e_0/2W_k),
    
    v_k   =   -Δ_k/√(2W_k(W_k-e_0)).

By virtue of the unitarity of (<ref>) it follows that the
electron number (<ref>) becomes

    N_e   =   2 ∑_ k,σ ( α_k,σ^†α_k,σ^
    + β_k,σ^†β_k,σ^ ).

The volume of the quasiparticle Fermi surface therefore corresponds
to both, conduction electrons and f electrons.
Despite the fact that all basis states have precisely
one f electron per cell, so that these are strictly localized, 
the f electrons do
contribute to the Fermi surface volume as if they  were 
itinerant<cit.>.

The bosonic Hamiltonian is (with J̃=J-λ)

    H_B   =   J̃ ∑_n  t_n^†·t_n^
    + ∑_m,n ζ̃_m,n t_n^†·t_m^
                   + ∑_m,n( η̃_m,n t_m^·t_n^ + H.c.),
    ζ̃_m,n   =   t_m,n/2 ∑_σ ⟨ b_m,σ^† b_n,σ^ - a_m,σ^† a_n,σ^ ⟩,
    η̃_m,n   =   t_m,n/2 ⟨ b_m,↑^† a_n,↓^†
    - b_m,↓^† a_n,↑^† ⟩.

Fourier transformation gives

    H_B   =   ∑_k (J̃+ζ̃_k)
     t_k^†·t_k^
     + ∑_k ( η̃_k t_k^·t_-k^
    + H.c.),

with ζ̃_k and η̃_k defined
as in (<ref>).
This can be solved by the ansatz



τ_k^†= ũ_k t_k^† + ṽ_k t_-k^
and H_B becomes

    H_B   =   ∑_k ω_k τ_k^†·τ_k^
    
    + 3/2 ∑_k (ω_k - (J̃ + ζ̃_k)).

Thereby

    ω_k   =   √((J̃+ζ̃_k)^2 - 4 η̃_k^2)
    ũ_k   =   2η̃_k/√(2ω_k(J̃+ζ̃_k-ω_k)),
    ṽ_k   =   √(J̃+ζ̃_k-ω_k/2ω_k).

The additive constant is Nc with

    c   =   -3J/4+λ (1-s^2)+μ n_e -Θ,
    Θ   =   ∑_α 2z_α/t_α ( ζ_αζ̃_α - 2
     η_αη̃_α ),

and the Helmholtz free energy becomes


    F    =    -2/β ∑_k ∑_ν∈{α,β} log(1+e^-β E_ν,k)
    + 3/β ∑_k log(1-e^-βω_k) 
     - 2∑_k (t̃_k - e_0) 
    
    + 3/2 ∑_k (ω_k - (J̃+ζ̃_k))
    + Nc.

Had we used the alternative form (<ref>) for H_J we would have obtained
the same expression with λ-3J/4→λ.
Minimizing F with respect to λ gives
1 - s^2 - n_F - n_B =0 where

    n_F    =   1/N∑_k,σ ⟨ b_k,σ^† 
    b_k,σ^ + a_k,σ^† a_k,σ^⟩
    =2/N∑_k[ 1 - e_0/W_k ( f(E_α,k) - f(E_β,k) )],
    
    n_B    =   1/N∑_k ⟨t^†·t^⟩
    = 3/N∑_k [J̃ + ζ̃_k/2ω_k (βω_k/2)
    -1/2],


are the densities of fermions and bosons, respectively.
Minimization with respect to the ζ and η parameters
in H_F and H_B gives the self-consistency equations

    ζ_α   =   3t_α/2N∑_k γ_α,k[J̃+ζ̃_k/2ω_k(βω_k/2) - 1/2],
    η_α   =   3t_α/2N∑_k γ_α,k η̃_k/ω_k(βω_k/2),
    ζ̃_α   =   t_α/N ∑_k γ_α,k [
    f(E_α,k) + f(E_β,k) -1],
    η̃_α   =   t_α/N ∑_k γ_α,k Δ_k/2W_k [ f(E_α,k) - f(E_β,k) ].

The above equations can also
be derived `directly', by evaluating the respective
thermal averages with the bosonic and fermionic mean-field Hamiltonian.
Differentiation with respect to s gives the additional condition

    λ   =   1/N ∑_k ϵ_k [ f(E_α,k) + 
    f(E_β,k) -1 ]
    
                  - 1/N ∑_k ϵ_k η_k/W_k [ f(E_α,k) - f(E_β,k) ],

where ϵ_ k is the noninteracting dispersion.
The resulting set of coupled equations can be solved numerically thereby using 
Broyden's algorithm<cit.> for better convergence.



For the Kondo-insulator with nearest-neighbor hopping particle hole-symmetry
results in ζ_ k=ζ̃_ k=0. Figure <ref> shows the
remaining parameters, η, s and λ as functions
of J/t. The parameter s reaches zero for J_min/t≈ 0.1173 
and there is no solution for smaller values of J/t.
The reason is that even for s=0 the parameters Δ_ k and
η̃ are finite and in fact
increase for small J/t so that the resulting density of particles,
n_B+n_F, exceeds 1 at J_min and (<ref>)
can no longer be fulfilled. This may be a consequence of the
fact that the bond particle formulation of the Kondo lattice ultimately
is a strong-coupling theory which is justified best for 
J/t→∞. It should also be noted that once 
the particle density n_B+n_F approaches unity, the bond particle theory is 
highly unreliable anyway. 







The mean-field expectation values η_ k and 
η̃_ k are small for J/t>1,
the parameter λ is relatively large and negative. 



Figure <ref> shows the bands E_ν, k
for the fermions and the dispersion ω_ k of the bosons.
The smallness of η̃ results in a small bandwidth for the bosons,
whereas the relatively large and negative λ results in a large
bandgap for the fermions which stays approximately constant
for J/t<1, as well as a considerable upward shift of the triplet dispersion. 
The band structure is consisten with the hybridization picture
with extended `heavy' band portions and is roughly consistent with numerical 
results results for the 1D PAM<cit.> and KLM<cit.> 
although the size of the gap comes out too large.
QMC has also shown a well-defined weakly dispersive and gapped mode in the 
dynamical spin correlation function of the PAM<cit.>, roughly 
consistent with the mean-field boson dispersion. 
The boson dispersion is symmetric with 
respect to k=π/2 whereas DMRG calculations find the maximum of the 
dispersion of the lowest triplet state at k=0, the minimum 
at k=π<cit.>. We define the quasiparticle gap, 
Δ_QP=E_0^(N+1)+E_0^(N-1) - 2E_0^(N) which we approximate by the band 
gap i.e. in 1D Δ_QP=E_β(k=0)-E_α(k=π).

So far we have ignored the terms H_3 and H_4 because a mean-field treatment
of these terms would result in some type of magnetic order.
To study the contribution of H_3 and H_4
as well as the unfactorized remainder of H_2 to the ground state energy 
at least approximately, we treat these terms in 2^nd order perturbation 
theory in analogy to Møller-Plesset perturbation theory<cit.>.
Since the Kondo insulator has a finite gap in its excitation spectrum
this is probably a good approximation.
More precisely, we take the mean-field Hamiltonian
H_MF=H_F+H_B+Nc as the unperturbed Hamiltonian -
its ground state is the product
of the ground states of H_F and H_B. The perturbation is

    H̃_1    =     H_2 + H_3 + H_4 - (H_2,MF - NΘ)

where H_2,MF is the mean-field factorized
form of H_2 i.e. the terms in H_F+H_B which are
∝η, ζ, η̃, ζ̃.
H̃_1 thus is a sum of terms of the form

    ∑_m,nt_m,n/2      ( O_m,n^(F) O_m,n^(B) 
    - O_m,n^(F) ⟨ O_m,n^(B)⟩. 
          .         - ⟨ O_m,n^(F)⟩ O_m,n^(B)
    +⟨ O_m,n^(F)⟩ ⟨ O_m,n^(B)⟩ )

where O_m,n^(F) (O_m,n^(B))
contain only fermion (boson) operators and ⟨ ..⟩
denotes expectation values in the mean-field ground state
(which are zero for terms arising from H_3 and H_4). Then 
H_MF+H̃_1 is the complete Hamiltonian (with added constraints)
and the first order correction ⟨H̃_1⟩=0.  
All matrix elements of H̃_1 between the mean-field ground state
|GS⟩ and states which contain either only a fermionic excitation - 
such as β_ k,σ^†α_ k,σ'^|GS⟩ - 
or only a bosonic excitation - such as
τ_ q,x^†τ_- q,y^†|GS⟩ - are zero.
It follows that in 2^nd order perturbation theory we may as well
take the perturbation to be H̃_1=H_2+H_3+H_4
but consider only intermediate states which contain both, 
a fermionic and a bosonic excitation.

We begin with H_3, which can be rewritten as

    H_3    =    -∑_n  (t_m,n/2  t^†_n· A_n  + H.c.)
     A_n    =   ∑_m [ ( b_m,n- a_m,n) s_m +
    s_m^† (π_m,n - π_n,m) ].

A considerable simplification comes about by noting that
since the mean-field expectation value η̃ is small,
resulting in an almost flat triplet dispersion,
ω_ k≈J̃, 
we may neglect all terms in the triplet Hamiltonian 
other than the energy term 
J̃∑_n  t_n^†· t_n^.
The ground state then is the vacuum for triplets and only
the terms ∝ t^†_n contribute to the energy correction.
In the Kondo insulator the operators A_n, being quadratic in the 
fermions, can only excite a quasiparticle from the lower to the upper 
quasiparticle band, say from momentum p to momentum q.
A typical state which couples to the ground state in this way
would be form β_ q,σ^†α_ p,σ^ t_n,z|GS⟩.



The unperturbed energy of this state is 
E_β, q + J̃ - E_α, p + E_0,
because the triplet which is created along with the particle-hole
pair contributes the energy J̃ if we neglect the dispersion of the 
triplets. The matrix element for the transition can be evaluated
by using (<ref>) and is 

    s/2N e^i( p- q)· R_n m_ p, q

with

    m_ p, q   =   ϵ_ q u_ p (u_ q+v_ q) -
    ϵ_ p u_ q (u_ p-v_ p).

The correction to the energy/site due to creation of a single
triplet then is

    δ E_0^(1)   =    -3s^2/2N^2 ∑_ p, q m_ p, q^2/E_β, q + J̃ - E_α, p.



We proceed to the correction due to H_2 and H_4. The parts
which give a nonvanishing result when acting onto the vacuum for triplets are

    H_2'    =    ∑_m,n t_m,n/2( a_n,↓b_m,↑ -
    a_n,↑b_m,↓ )  t^†_m· t^†_n
    
    H_4'    =    -i ∑_m<n t_m,n/2( π_m,n - π_n,m )
    · ( t^†_m× t^†_n )

The states which can be reached have the form
β_ q,σ^†α_ p,σ'^ t_n,x^† t_m,x'^†|GS⟩.
Proceeding as above and
specializing to a 1D chain with only nearest neighbor hopping -t
we find for the energy shift due to the creation of two
triplets

    δ E_0^(2)   =   -9t^2/N^2 ∑_ p, q (1 + cos( p+ q)) u_ p^2 u_ q^2/E_β, q + 2J̃ - E_α, p.

The different contributions to the ground state energy are shown in
Figure <ref>.
The perturbation correction is quite small which indicates that the use of
perturbation theory is adequate. 
The energy shift due to creation of a single triplet, 
δ E_0^(2) is small but still of order 0.1t, 
whereas δ E_0^(3) is negligible, of order 10^-2t.




Lastly, we discuss an improved calculation of the triplet dispersion
ω_ k. While the mean-field calculation predicts these to
be almost dispersionless the term H_3 gives a more substantial dispersion.
We make the variational ansatz for a z-triplet-like excitation
with momentum q

    |Ψ_ q⟩   =   [ a_ q t_ q^† + 1/√(2N)∑_ k,σ,σ'
    b_ q, k β_ k+ q,σ^†τ_σ,σ'α_ k,σ'^ ]|0⟩,

with variational parameters a_ q and b_ q, k.
We obtain the triplet frequency

    ω̃_ q   =   J̃ +s^2/2N ∑_ k V^2( k,q)/ω̃_ q-(E_β, k+q -
    E_α, k)
    
    V( k,q)   =   
    (ϵ_ k+q-ϵ_ k) v_ k+q v_ k
                  + ϵ_ k+q u_ k+q v_ k
    + ϵ_ k u_ kv_ k+q,

where ϵ_ k is the noninteracting dispersion of the
conduction electrons.
Numerical evaluation shows that this always takes its minimum
at k=π and the maximum at k=0 - consistent with
DMRG<cit.>. The energy ω̃_q=π thus
is the energy of the lowest S=1 excitation, called the spin gap, 
Δ_s. The bandwidth of the spin excitations is 
W_s=ω̃_q=0-ω̃_q=π.
Figure <ref> compares the dependence of the ground state energy per site,
the quasiparticle and spin gap and the bandwidth of the spin excitations
on t/J for the 1D Kondo insulator to results obtained by DMRG<cit.>. 
As expected, the perturbation correction improves the 
ground state energy/site, which is reasonably close to the numerical values
for t/J<1. For all other quantities the calculated energies
deviate from the numerical results already for relatively large J/t. 
One might 
wonder if this is the consequence of the additional approximation
to neglect η̃, but even for J/t=1 where the boson bandwidth
is quite small (see Figure <ref>) the deviation for the
spin gap is already substantial.
To conclude this section we discuss the differences to the previous mean-field treatment 
in Ref.<cit.>. In this work, both the singlet and
the triplet operators in (<ref>) were replaced by c-numbers,
s_n^†,s_n^→ s and t_n^†,t_n^→t e^iQ·R_n, thus reducing (<ref>) to a quadratic form
from the outset,
and then minimizing F with respect to s and t.
The dynamics of the spin excitations thereby was not studied.



§ RENORMALIZED ENERGY OF FORMATION


We consider a different approximation scheme to account for
the constraint whereby we consider sites occupied by
a singlet as `empty'. This is equivalent to working in a fictitious
Hilbert space for the fermionic particles a^†_n,σ and
b^†_n,σ as well as the triplets t^†_n,α,
whereby the states in this fictitious Hilbert space correspond to those of 
the physical Kondo lattice according to the rule

    ∏_i∈ S_aa_i,σ_i^† ∏_j∈ S_bb_j,σ_j^† ∏_l∈ S_tt_l,x_l^†|0⟩→                           
    ⊗_i∈ S_a|a,i,σ_i⟩ ⊗_j∈ S_b|b,j,σ_j⟩ ⊗_l∈ S_t|t_x_l,l⟩ ⊗_n∈ S_s |s_n⟩.

S_a, S_b and S_t denote the set of sites occupied by a hole-like fermion,
an electron-like fermion or a triplet, respectively, and
S_s=(S_a∪ S_b∪ S_t)^C the set of remaining sites.
In other words, all sites not occupied by a fermion or triplet
are filled up with `inert' singlets. The Hamiltonian - and all other operators 
in the bond particle representation - then can be obtained from
(<ref>) by replacing all singlet operators by unity. 
Only the form (<ref>) of the exchange term can be used.
For this representation to
make sense we again have to impose the constraint that no two particles
of any type occupy the same site, because the resulting state could not
be translated meaningfully to a state of the physical Kondo lattice
via (<ref>). Assuming that the density of bond particles is small,
however, we relax again this constraint. 

On the other hand, if the constraint were rigorously enforced, presence of 
any one particle -  be it a^†_n,σ, b^†_n,σ or 
t^†_n,α - at a given site n would prevent all remaining 
terms in the Hamiltonian which involve creation or annihilation of any other 
particle at site n from acting, resulting in a loss of kinetic
energy. The constraint 
thus increases the cost in energy for adding a fermion or boson.
Accordingly, in (<ref>) the energy for adding a fermion therefore should be
e_0=3J/4+κ rather than 3J/4 and the energy for adding
a boson J̃=J+κ rather than J, where κ is some as 
yet unspecified loss of kinetic energy.
Actually κ may be expected to be different for fermions and bosons.
We will
discuss possible estimates for κ later on. It should
also be noted that such an increase of the 
energies of formation of the particles would reduce their densities and thus 
make relaxing the constraint of no double occupancy an even better 
approximation.
The mean-field theory outlined in the previous section and the
approximation scheme discussed in the present section mimick the constraint
in different ways: mean-field theory amounts to a 
Gutzwiller-like downward renormalization
of the hopping integrals whereas the present scheme amounts to a
higher energy of formation of the particles.
Collecting all terms which become quadratic when we drop the singlets
we obtain the noninteracting Hamiltonian


    H_0   =   ∑_n ( e_0 ∑_σ (b_n,σ^† b_n,σ^ 
    + a_n,σ^† a_n,σ^) + J̃ t_n^†·t_n^ ) +∑_m,n t_m,n/2 ∑_σ ( b_m,σ^† b_n,σ^ -
    a_m,σ^† a_n,σ^ )
              
    -∑_m,n t_m,n/2 [ ( 
    b_m,↑^† a_n,↓^†
    - b_m,↓^† a_n,↑^† ) + H.c. ] - 3NJ/4.

The interaction part of the Hamiltonian is the sum of

    H_2   =   ∑_m,n t_m,n/2 ∑_σ ( b_m,σ^† b_n,σ^ -
    a_m,σ^† a_n,σ^ ) t_n^†·t_m^ 
    +∑_m,n t_m,n/2 [ ( 
    b_m,↑^† a_n,↓^†
    - b_m,↓^† a_n,↑^† ) t_m^·t_n^ + H.c ]
    
    H_3   =    - ∑_m,n t_m,n/2 [ ( π^†_m,n·( t_n^ - t_m^) + H.c. )
    +( b_m,n - a_m,n )·
    (t_n^† + t_m^ ) ],
    
    H_4   =   ∑_m,n t_m,n/2 [ [ iπ^†_m,n·
    ( t_m^×t_n^ )
    + H. c.  ]
    - i ( b_m,n - a_m,n ) 
    ( t_n^†×t_m^ ) ].


The part H_0 was used in Refs.<cit.>.
Due to particle-hole symmetry the extra Lagrange multiplier introduced 
in these Refs. to enforce consistency of
the c-like spectral weight with N_e is not necessary here.
We now proceed as in the case of mean-field theory, that means
first diagonalize H_0 to obtain the band structure and treat the interaction 
terms in the same approximation as there, i.e. in perturbation theory for the 
ground state energy and using the variational ansatz (<ref>) for the spin 
excitations.
The fermionic part H_F again can be diagonalized by the unitary 
transformation (<ref>) with the result

    H_F   =   ∑_k,σ ( E_α,k α_k,σ^†α_k,σ^
    + E_β,k β_k,σ^†β_k,σ^ )
                                          
    - 2∑_k (ϵ_k/2 - e_0), 
    
    E_ν,k   =   ϵ_k/2± W_k - μ,
    
    W_k   =   √( e_0^2 + (ϵ_k/2)^2),
    
    u_k   =   -√(W_k-e_0/2W_k),
    
    v_k   =   -ϵ_k/2/√(2W_k(W_k-e_0)).

The noninteracting ground state for the bosons is the bosonic
vacuum i.e. the bosons do not contribute to the
ground state energy E_0 in this approximation.
The Helmholtz Free energy is

    F    =    -2/β ∑_k ∑_ν∈{α,β} log(1+e^-β E_ν,k) 
                  + N(2e_0 -3J/4) + μ N.

We can obtain the expectation value of the kinetic energy
by multiplying all hopping integrals by a parameter χ: 
t_m,n→χ t_m,n and forming 
∂ F/∂χ|_χ=1, with the result:

    k   =   1/N ∑_ k [ ϵ_ k( 1-ϵ_ k/2W_ k ) f(E_α, k)
    .
    
          .                 + ϵ_ k( 1+ϵ_ k/2W_ k ) f(E_β, k)].

Treating the interaction part H_2+H_3+H_4 in 2^nd order perturbation theory
a slight modification occurs. Since the ground state is the vacuum for
bosons and the vector operators a, b and π
have zero expectation value in the fermionic ground state (which is spin 
singlet) only intermediate states which contain both, a bosonic and a 
fermionic excitation, can be reached from the ground state by acting with
H_3 or H_4. The energy shifts due to such doubly excited intermediate
states take the form (<ref>), (<ref>), but with s→ 1 
in (<ref>). In addition, however, the term H_2 also has a
nonvanishing matrix element with states of the form
t_m^†· t_n^† |GS⟩, because the
fermionic factor of the corresponding term is a singlet, which does have
a nonvanishing ground state expectation value.
For the 1D chain with nearest neighbor hopping
this gives an additional contribution to
δ E_0^(2) of -3 I^2/2J̃, whereby

    I    =   -t^2/N ∑_q cos^2(q)/W_q.

Numerical evaluation shows that this contribution is quite substantial
and obviously this replaces the energy gain due to
the mean-field factorized terms in the previous section.
Finally, the equation for the dispersion of the spin excitation takes the
form (<ref>), again with s→ 1.




Lastly, we consider the value of κ,
the correction to the energies of formation of the fermions and bosons.






We switch to a phenomenological approach and
approximate κ≈ x k, i.e. a dimensionless
parameter x times the kinetic energies of the Fermions/site, given
in (<ref>). With x fixed, 
κ has to be determined self-consistently for each J/t. 
We neglect the loss of kinetic energy of bosons which is 
reasonable for J/t>1 where the boson density is low. In fact, as will be
shown in a moment, we can obtain good agreement with numerics over the whole 
range J/t>1 by choosing x≈ 0.4 independent of JJ/t.
Varying 0.2<x<0.6 thereby does not deteriorate the 
agreement significantly, so that also an x which varies with J/t - which
is actually what one might expect - would give similar results.

To begin with, Figure <ref> shows the various contributions to the ground 
state  energy/site obtained with this choice of x and demonstrates that at 
least for J/t>1 the perturbation correction is small as it should be. 
Figure <ref> compares characteristic energies of the system as functions 
of t/J to numerical results. For larger t/J>1 the agreement is
poor, in particular for t/J ≈ 1 the spin excitation energy 
ω̃_q from the variational ansatz becomes negative around 
q=π indicating the failure of the calculation. Accordingly, results for 
the spin gap Δ_s and spin excitation bandwidth W_s are shown only up 
to this value of t/J.
W_s comes out quite good for t/J< 1 - the Figure also shows the
very good agreement between DMRG and  SE for W_s at t/J=0.25 in (d).
DMRG and SE also agree very well for Δ_s so that we do not show 
the SE results in (c). Figure <ref> also shows results obtained by
perturbation expansion in t/J<cit.>. 
As one might have expected the DMRG results and bond particle theory approach 
these for t/J→ 0. The ground state energy is reproduced
remarkably well be the perturbation expansion, but all other
characteristic energies deviate substantially from the
perturbation expansion for t/J→ 1. This shows that despite
being a strong coupling theory by nature, bond
particle theory does go beyond simple perturbation theory.
Figure <ref> shows the same characteristic
energies but now plotted versus J/t in the range J/t>1. It is obvious
that in this range the agreement between bond particle theory
and numerics is quite good.
Figure <ref> shows the densities of fermions and bosons,
n_F and n_B. Thereby n_B is
obtained from n_B=∂ E_0/∂J̃ where
E_0 is the ground state energy per site
including the second order perturbation 
correction. The main contribution thereby comes from 
(<ref>). The data points for t/J=1 are DMRG 
results<cit.>. The densities are small for t/J<1 and bond 
particle theory somewhat underestimates the densities of the particles
at t/J=1. For larger t/J the densities increase rapidly and
the sum n_B+n_F exceeds 1 at t/J≈ 2, indicating the
breakdown of bond particle theory.
Figure <ref> shows the band structure of the fermions and the
dispersion of the bosons obtained from the variational ansatz.
For J/t=2 the dispersions can be compared to results obtained
by SE taken from Ref. <cit.>.
While the quasiparticle gap is approximately correct,
the `heavy' part of the band structure has too much dispersion
as compared to SE, and the bandwidth is slightly overestimated.
The dispersion of the spin excitations is reasonably correct
for J/t=2, both the spingap, the overall form of the dispersion
and the bandwidth compare quite well with the SE result.






The combined DMRG and SE data in Figure <ref> suggest
that there is a crossover between two regimes at around t/J≈ 1: 
the Δ_s/J vs. J/t curve drops rapidly for J/t>1 but then bends 
sharply at J/t≈ 1 and Δ_s/J is small but finite for J/t<1. 
Similarly, the band width of the spin excitations, W_s/J increases
with decreasing J/t in the range J/t>1 but then must drop sharply at 
J/t≈ 1. This may indicate a crossover from a strong coupling
regime for J/t>1 where the system apparently can be described well
by the bond particle theory, to a weak coupling regime for J/t<1
where maybe mean-field theories work better.
It should also be noted that the present theory must fail in the limit 
J/t→ 0 not only because the bond particle density increases
sharply but also because for because J/t→ 0 the quasiparticle gap 
vanishes, whereas the energy e_0 - which determines the 
magnitude of the gap - cannot approach zero for any x>0.
All in all Figure <ref> indicates that for t/J<1 the bond-particle 
Hamiltonian (<ref>) with suitably renormalized
energies of formation gives a reasonably correct description of the
low-energy elementary excitations of the 1D Kondo insulator.



§ SUMMARY AND DISCUSSION


In summary we have derived an exact representation of the Kondo lattice
model in terms of bond particles: fermions corresponding to
unit cells with an odd number of electrons and bosons corresponding to
cells with two electrons coupled to a singlet or triplet.
Thereby the constraint to have
precisely one f electron/cell, which considerably complicates the solution of
the KLM, is fulfilled automatically and replaced by
the constraint to have precisely one bond particle per site.
If the singlet bosons are considered as condensed or the singlet is
defined as the vaccum state of a cell, this constraint becomes an infinitely 
strong Hubbard-like repulsion, but for a relatively dilute system of particles, 
so that it may be justified to relax it
(for a system of low density even an infinitely strong repulsion can be treated 
diagrammatically<cit.>).
The requirement of low particle density is indeed fulfilled for J/t>1.

We have discussed two schemes to approximately
incorporate effects of the remaining Hubbard repulsion between bond particles 
into their Hamiltonian. First, mean-field theory, where the singlets are taken 
as condensed, amounts to a Gutzwiller-like downward renormalization of 
all hopping integrals. Second, a scheme where the singlet is considered as
the vacuum state of a site and the constraint is mimicked by adding the 
loss of kinetic energy, which incurs due to the 
blocking of a site by a bond particle, to the energy ascribed
to the respective particle. Approximating this loss of kinetic energy as 
the kinetic energy of fermions per site times a phenomenological constant
of order unity allowed to reproduce numerical results ontained by 
density matrix renormalization group
and series expansion calculations for the 1D Kondo insulator in the range 
J/t>1 with good accuray.
Thereby relatively simple techniques were used
- 2^nd order perturbation theory for the ground state energy and
the simplest possible variational wave function for the triplet dispersion -
to produce these results. 
The good agrement with numerics in the range J/t>1 for a variety of
quantitiesis then is a strong indication that in this parameter range 
the triplets and fermions indeed correspond to the approximate elementary 
excitations and this is the main result of the present paper.
Despite being a somewhat lengthy
expression the bond particle Hamiltonian with renormalized particle
energies appears to be useful for quantitative calculations.
Of course the phenomenological approach used here is somewhat unsatisfactory
and a more rigorous calculation following Refs. <cit.>
would be desirable.

The question then arises, whether J/t>1 is a sufficient range of validity
to discuss magnetic ordering and quantum critical points in the KLM.
For the 2D square lattice with nearest neighbor hopping it is
known that at n_e=2 and T=0 antiferromagnetic ordering occurs for
J/t ≤ J_c/t=1.45<cit.>, i.e. for relatively large J/t
(although with the information at hand we cannot say much about the range 
of vailidity of the bond particle description in higher dimensions).
As pointed out by Sachdev and Bhatt<cit.> bond particle theory 
gives a rather natural description of magnetic ordering, namely
the condensation of triplets into a momentum corresponding
to the magnetic ordering wave vector. Applying this description
of antiferromagnetic ordering in the KLM has already produced encouraging
results: using the mean-field version of bond particle theory
with condensed triplets, Jurecka and Brenig found J_c/t=1.5, quite
close to the exact value. This is even more remarkable in that even 
numerical methods appear to have difficulties to accurately reproduce
J_c/t: VMC gives J_c,1/t=1.7<cit.>, DCA gives 
J_c,1/t=2.1<cit.>, and DMFT gives 
J_c,1/t=2.2<cit.>.
Moreover, in Ref. <cit.> it was shown that using unrenormalized 
energies of formation (i.e. κ=0) bond particle theory for the planar 
KLM did give the too large J_c/t=2.20 but reproduced the phase diagram of 
the model in the 
(J/t,n_e)-plane obtained by VMC<cit.> and
DMFT<cit.> quite well if J was measured in units of J_c - i.e. if the 
phase diagram was plotted in the (J/J_c,n_e) plane - so that the error in 
J_c/t cancelled to some degree. This is 
encouraging in that the phase diagram of the KLM in 2D is quite 
intricate, comprising the paramagnetic and two antiferromagnetic phases 
with different Fermi surface topology, with various 1^st and 2^nd order
transitions between them. Also, the band structure in the AF phase and its
change with J/J_c as obtained by DCA<cit.> could be
reproduced in this way<cit.>. It should also be noted that in
the above bond particle calculations antiferromagnetic order appears
without an additional Heisenberg exchange between f spins, that means it
comes about solely by the interaction mediated by the conduction electrons.



Stewart
G. R. Stewart, Rev. Mod. Phys. 56, 755 (1984).
kondoinsulators
P. A. Lee, T. M. Rice, J. W. Serene, L. J. Sham, and J. W. Wilkins, 
Comm. in Condensed Matter Phys. 12, 99 (1986).
StewartII
G. R. Stewart,Rev. Mod. Phys. 73, 797 (2001).
loenireview
H. v. Löhneysen, A. Rosch, M. Vojta, and P. Wölfle,
Rev. Mod. Phys. 79, 1015 (2007).
Steglichreview 
Q. Si and F. Steglich, Science 329, 1161 (2010).
SchriefferWolff
J. R. Schrieffer and P. A. Wolff, Phys. Rev. 149, 491 (1966).
Yosida
K. Yosida, Phys. Rev. 147, 223 (1966).
Varma
C. M. Varma and Y. Yafet,
Phys. Rev. B 13, 2950 (1976).
Gu
O. Gunnarsson and K. Schönhammer, Phys. Rev. B 28, 4315 (1983).
YoshimoriSakurai
A. Yoshimori and A. Sakurai,
Progr. Theor. Phys. Supp. 46, 162 (1970).
ReadNewns
N. Read and D. M. Newns, J. Phys. C 16, 3273 (1983).
Keiter
H. Keiter and J. C. Kimball, Int. J. Magn. 1, 233 (1971).
Kuramoto
Y. Kuramoto, Z. Phys. B 53, 37 (1983).
Coleman
P. Coleman, Phys. Rev. B 29, 3035 (1984).
Bickers
N. E. Bickers, Rev. Mod. Phys. 59, 845 (1987).
RG
K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
Bethe
N. Andrei, K. Furaya, and J. H. Lowenstein, 
Rev. Mod. Phys. 55, 331 (1983).

Oshikawa
M. Oshikawa, Phys. Rev. Lett. 84, 3370 (2000).

LacroixCyrot 
C. Lacroix and M. Cyrot, Phys. Rev. B 20, 1969 (1979).
Lacroix
C. Lacroix, Solid State Commun.  54, 991 (1985).
AuerbachLevin 
A. Auerbach and K. Levin, 
Phys. Rev. Lett.  57, 877 (1986).
Burdinetal 
S. Burdin, A. Georges, and D. R. Grempel, 
Phys. Rev. Lett.  85, 1048 (2000).
ZhangYu 
G.-M. Zhang and L. Yu, Phys. Rev. B 62, 76 (2000).
Lavagna 
M. Lavagna and C. Pepin, Phys. Rev. B 62, 6450 (2000).
Senthil 
T. Senthil, M. Vojta, and S. Sachdev,
Phys. Rev. B 69, 035111 (2004).
Global 
M. Vojta, Phys. Rev. B 78,  125109 (2008).
ZhangSuLu 
G.-M. Zhang, Y.-H. Su, and L. Yu,
Phys. Rev. B 83, 033102 (2011).
Nilsson 
J. Nilsson, Phys. Rev. B 83, 235103 (2011).
Montiel
X. Montiel, S. Burdin, C. Pepin, and A. Ferraz,
Phys. Rev. B 90, 045125 (2014).

MillisLee
A. J. Millis and P. A. Lee, Phys. Rev. B  35, 3394 (1987).
ReadNewns1
D. M. Newns and N. Read, Adv. Phys. 36, 799 (1987).

RiceUeda 
T. M. Rice and K. Ueda, Phys. Rev. Lett. 55, 995 (1985).
Fazekas 
P. Fazekas and E. Müller-Hartmann, Z. Phys. B 85, 285 (1991).

yuwhite
C. C. Yu and S. R. White, Phys. Rev. Lett. 71, 3866 (1993).
The DMRG data for the ground state energy are not published in this reference 
but seem to be given only in Ref. <cit.>.
MC1
S. Moukouri and L. G. Caron, Phys. Rev. B. 52,15723(R) (1995).
MC2
S. Moukouri and L. G. Caron, Phys. Rev. B. 54,12212 (1996).
Mutou
T. Mutou, N. Shibata, and K. Ueda, Phys. Rev. Lett. 81, 4939 (1998). 
Smerat
S. Smerat, U. Schollwöck, I. P. McCulloch, and H. Schoeller
Phys. Rev. B 79, 235107 (2009).

J1
M. Jarrell, H. Akhlaghpour, and Th. Pruschke
Phys. Rev. Lett. 70, 1670 (1993).
J2
M. Jarrell, Phys. Rev. B 51, 7429 (1995).
J3
A. N. Tahvildar-Zadeh, M. Jarrell, and J. K. Freericks,
Phys. Rev. Lett.  80, 5168 (1998).
peterspruschke 
R. Peters and T. Pruschke
Phys. Rev. B 76, 245101 (2007).

Carsten
C. Gröber and R. Eder, Phys. Rev. B. 57,12659(R) (1998).
Assaad
F. F. Assaad, Phys. Rev. Lett. 83, 796 (1999).
series
Z.-P. Shi, R. R. P. Singh, M. P. Gelfand, and Z. Wang
Phys. Rev. B 51, 15630(R) (1995).
seriesexp
W. Zheng and J. Oitmaa, Phys. Rev. B 67, 214406 (2003).
Trebst
S. Trebst, H. Monien, A. Grzesik, and M. Sigrist, 
Phys. Rev. B 73, 165101 (2006).

WatanabeOgata
H. Watanabe and M. Ogata, Phys. Rev. Lett. 99, 136401 (2007).
Asadzadeh
M. Z. Asadzadeh, F. Becca, and M. Fabrizio, 
Phys. Rev. B 87, 205144 (2013).
Kubo
K. Kubo, J. Phys. Soc. Jpn. 84,  094702 (2015).

Tsunetsugu
H. Tsunetsugu, Y. Hatsugai, K. Ueda, and M. Sigrist,
Phys. Rev. B 46, 3175 (1992).
White
 J. A. White, Phys. Rev. B 46, 13905 (1992).
Tsutsui
K. Tsutsui, Y. Ohta, R. Eder, S. Maekawa, E. Dagotto, and J. Riera,
Phys. Rev. Lett. 76, 279 (1996).

MartinAssaad
L. C. Martin and F. F. Assaad,
Phys. Rev. Lett. 101, 066404 (2008).
MartinBerxAssaad 
L. C. Martin, M. Bercx, and F. F. Assaad,
Phys. Rev. B 82, 245105 (2010).

Lenz
B. Lenz, R. Gezzi, and S. R. Manmana,
Phys. Rev. B 96, 155119 (2017).
Doniach
S. Doniach, Physica B 91, 231 (1977).
RKKY
M. A. Ruderman and C. Kittel
Phys. Rev. 96, 99 (1954); 
T. Kasuya, Progress of Theoretical Physics, 16, 45 (1956);
K. Yosida, Phys. Rev. 106, 893 (1957).
K1 
Q. Si, S. Rabello, K. Ingersent, and J. Smith, 
Nature 413, 804 (2001).
K2 
Q. Si, S. Rabello, K. Ingersent, and J. L. Smith, 
Phys. Rev. B 68, 115103 (2003).
K3 
L. De Leo, M. Civelli, and G. Kotliar, 
Phys. Rev. Lett. 101, 256404 (2008).
K4 
L. Zhu and Q. Si, Phys. Rev. B 66, 024426 (2002).
K5 
J. X. Zhu, D. R. Grempel, and Q. Si, 
Phys. Rev. Lett. 91, 156404 (2003).
K6 
P. Sun and G. Kotliar, 
Phys. Rev. Lett. 91, 037209 (2003).
K7 
S. J. Yamamoto and Q. Si, 
Phys. Rev. Lett. 99, 016401 (2007).
K8 
M. T. Glossop and K. Ingersent, 
Phys. Rev. Lett.  99, 227203 (2007).
K9 
J.-X. Zhu, S. Kirchner, R. Bulla, and Q. Si, 
Phys. Rev. Lett. 99, 227204 (2007).
K10 
Q. Si, J. H. Pixley, E. Nica, S. J. Yamamoto,
P. Goswami, R. Yu, and S. Kirchner,
J. Phys. Soc. Jpn. 83, 061005 (2014).
K11 
E. Abrahams, J. Schmalian, and P. Wölfle,
 Phys. Rev. B 90, 045105 (2014).
K12
R. Peters and N. Kawakami,
Phys. Rev. B 92, 075103 (2015).
K13
R. Peters and N. Kawakami,
Phys. Rev. B 96, 115158 (2017).

SachdevBhatt
S. Sachdev and R. N. Bhatt, Phys. Rev. B 41, 9323 (1990).
Gopalan
S. Gopalan, T. M. Rice, and M. Sigrist,
 Phys. Rev. B 49, 8901 (1994).
vojta1
M. Vojta and K. W. Becker, Phys. Rev. B 60, 15201 (1999).
vojta2
S. Ray and M. Vojta, Phys. Rev. B 98, 115102 (2018).
sushkov
O. P. Sushkov, Phys. Rev. B 60, 3289 (1999).
Park
K. Park and S. Sachdev, Phys. Rev. B 64, 184510 (2001).
Siahatgar
M. Siahatgar, B. Schmidt, G. Zwicknagl, and P. Thalmeier,
New J. Phys. 14, 103005 (2014).
Oana
R. Eder, O. Stoica, and G. A. Sawatzky, Phys. Rev. B. 55, R6109 (1997);
R. Eder, O. Rogojanu, and G. A. Sawatzky, Phys. Rev. B. 58, 7599 (1998).
JureckaBrenig
C. Jurecka and W. Brenig, Phys. Rev. B 64, 092406 (2001).
KotovHirschfeld
V. N. Kotov and P. Hirschfeld,
Physica B: Condensed Matter, 312-313, 174 (2002).
afbf
R. Eder, K. Grube, and P. Wróbel, Phys. Rev. B 93, 165111 (2016).
afbf1
R. Eder and P. Wróbel, Phys. Rev. B 98, 245125 (2018).
Nozieres
P. Nozieres, Eur. Phys. J. B 6, 447 (1998).

FW
See, e.g., A. L. Fetter and J. D. Walecka,
Quantum Theory of   Many   Particle   Systems
(McGraw-Hill,   New   York, 1971).
Kotov
V. N. Kotov, O. Sushkov, ZhengWeihong, and J. Oitmaa, 
Phys. Rev. Lett. 80, 5790 (1998).
Shevchenko
P. V. Shevchenko, A. W. Sandvik, and O. P. Sushkov
Phys. Rev. B 61, 3475 (2000).
uedareview
H. Tsunetsugu, M. Sigrist, and K. Ueda
Rev. Mod. Phys. 69, 809 (1997).
Broyden
D. D. Johnson, Phys. Rev. B. 38, 12807, (1988).
MP
C. Møller and M. S. Plesset, Phys. Rev. 46  618 (1934).

