Moment screening in the correlated Kondo lattice model

The magnetic correlations, local moments and the susceptibility in the correlated 2D Kondo lattice model at half filling are investigated. We calculate their systematic dependence on the control parameters J_K/t and U/t. An unbiased and reliable exact diagonalization (ED) approach for ground state properties as well as the finite temperature Lanczos method (FTLM) for specific heat and the uniform susceptibility are employed for small tiles on the square lattice. They lead to two major results: Firstly we show that the screened local moment exhibits non-monotonic behavior as a function of U for weak Kondo coupling J_K. Secondly the temperature dependence of the susceptibility obtained from FTLM allows to extract the dependence of the characteristic Kondo temperature scale T* on the correlation strength U. A monotonic increase of T* for small U is found resolving the ambiguity from earlier investigations. In the large U limit the model is equivalent to the 2D Kondo necklace model with two types of localized spins. In this limit the numerical results can be compared to those of the analytical bond operator method in mean field treatment and excellent agreement for the total paramagnetic moment is found, supporting the reliability of both methods.


Introduction
The understanding of strong correlations in mixed valent and heavy f-electron compounds is mostly based on two generic models described by the Anderson lattice Hamiltonian or, in the special case of almost integer valence, by the Kondo lattice Hamiltonian [1,2]. The latter model which will be the subject of this work represents an extreme limit where correlations in the f-orbitals are taken as infinitely large due to the Coulomb integral U f ≫ t with t denoting the hopping energy of conduction (c) electrons. On the other hand the Coulomb interaction U of conduction electrons (and also the one between c-and f-electrons) is completely neglected in this model. The U = 0 Kondo and Anderson models have the advantage of a simple and meaningful mean field solution [2,3] with the constraint of only singly occupied f orbitals implemented by an auxiliary boson in mean field approximation. In the lattice this leads to hybridized quasiparticle bands with an exponentially reduced hybridization gap. Close to the half filled case with one c-and f-electron per site they have an effective mass m * ≫ m b much larger than the conduction electron band mass m b [1]. The mean field Kondo lattice model may be merged with local density band structure calculations in the renormalized band theory [4] leading to a powerful method to calculate realistic Fermi surfaces for Ce and Yb heavy fermion compounds, including crystalline electric field effects.
It has been known that a strong Coulomb repulsion among the conduction electrons significantly alters the electronic properties of a metal. First and foremost, the kinetic energy as measured by the (effective) band width will be reduced and eventually vanish at a metal-to-insulator transition. Second, the trend towards localization is accompanied by the appearance of magnetic correlations.
The central focus of the present paper is the influence of the above-mentioned correlation effects on the screening of local moments. Of particular interest is the question how the Kondo energy scale is affected by conduction electron correlations. We adopt a two-dimensional Hubbard model at half-filling for the conduction electrons where correlation effects have been found to affect strongly the electronic properties. As the interacting conduction electron Hamiltonian, i. e., the two-dimensional Hubbard model cannot be solved exactly for an infinitely extended system we extract the relevant information from suitably chosen finite clusters. This approach seems justified considering the local character of the quantities to be investigated. The relevance of the cluster approach is assessed by comparing its results to the predictions of a constrained mean-field theory for an infinite system in the limit of extremely strong Coulomb repulsion.
Attempts to include the effect of correlations between conduction electrons which may become important when the latter originate from d-orbitals have sofar been mostly limited to the Kondo impurity models using various analytical techniques like perturbation theory [5,6] and 1/N f expansion for small U, a Schrieffer-Wolff approach [7] and RVB ansatz [5] for large U, a scaling approach [8] and also numerical techniques like NRG [9]. In the impurity problem it was concluded that the Kondo energy scale T * increases with U [5,6].
Concerning concentrated systems with magnetic moments at every lattice site the situation is rather controversial. In the case of non-interacting conduction electrons (U=0) the electronic properties will be determined by the competition between the energy gain due to (local) Kondo singlet formation and to developing (long-range) magnetic correlations. It is to be expected that the subtle balance between these two tendencies will be affected by conduction electron repulsion.
The correlated Anderson lattice problem was treated with a Gutzwiller variational method [10] and the 1D correlated Kondo lattice with DMRG approach [11]. It was further investigated with QMC simulations and a fermionic bond operator method [12] used before in the uncorrelated case [13]. The strongly correlated (U ≫ t) 'Kondo necklace' limit of the Kondo lattice was treated using exact diagonalization [14]. It was suggested [10] that the Kondo scale actually decreases for increasing U, in contrast to the impurity results. The competition of singlet formation and induced inter-site RKKY coupling and its signature in thermodynamics was studied in Ref. [15] using ED and FTLM methods for finite clusters of the 2D uncorrelated Kondo lattice and in Ref. [16] for 1D chains in the strongly correlated limit.
The present work which uses the unbiased numerical ED and FTLM methods as well as the analytical bond operator technique has two clear objectives: i) Firstly to investigate the characteristic T=0 on-site and inter-site correlations and local paramagnetic moment µ 2 loc systematically in the whole (J K , U ) parameter range of the correlated Kondo lattice model, where J K is the local antiferromagnetic exchange coupling. In particular the total local paramagnetic moment is an excellent measure for the formation or breakup of the Kondo singlet state as function of control parameters [15]. We will show that it exhibits an unexpected non-monotonic behaviour that has not been reported before. ii) Secondly we calculate the finite temperature behaviour of susceptibility and specific heat with FTLM on a small tile of the square lattice. From the maximum position we may extract the characteristic 'Kondo' temperature scale T * and in particular we investigate its systematic dependence on the correlation strength U. We will resolve the ambiguity mentioned before and show that T * increases monotonically with U. The ED and FTLM calculations will be performed on small Kondo clusters with open boundary conditions which have a more realistic one-particle density of states than the one with periodic boundary conditions [15].
We stress from the outset that in the context of the FTLM calculations on a finite cluster the meaning of T * is that of an average energy of low energy singlet-triplet excitations obtained from the maximum in the temperature dependent susceptibility and specific heat. This is indeed the way how an estimate of T * in a real concentrated Kondo compound may be obtained. It is not the same as the genuine single ion Kondo temperature T K in the continuum limit which is exponentially small compared to the hopping energy. This difference between lattice and single ion Kondo scales appears already within approximate analytical theories [2].
Nevertheless a comparison of ED numerical results for finite clusters with analytical results obtained by the bond operator method for the extended lattice and in the large U limit is useful and legitimate. Since the former are exact it allows one to judge how far the large U model approximation is justified and give an estimate of the reliability and accuracy of the bond operator approximation. This theory has previously been used for quantum critical properties of the localized model in the large U limit [17,18,19]. A generalized fermionic version has also been applied to the finite U case [13]. In this comparison we focus on the paramagnetic local moment size and on-site correlations as function of U and on the finite temperature susceptibility.
We note that the investigation of finite Kondo atom clusters on nonmagnetic surfaces is also of potential experimental interest. Sofar single magnetic Co adatoms embedded in metallic Cu n clusters [20] and two site Co-Co Kondo clusters connected by Cu n chains [21] have been investigated. However in these cases of a metallic substrate it is the simple conduction electron density of the substrate surface that plays the dominant role for the Kondo properties. The present investigation would be relevant for correlated Kondo clusters weakly adsorbed on an inert insulating substrate without chemical bonding of substrate and Kondo atoms. This has, to our knowledge, not been experimentally realized until now. Nevertheless it is clear that the investigation of finite size Kondo clusters is of great interest in itself and not only in view of the bulk Kondo lattice materials.

Model definition and single particle spectrum
The correlated Kondo lattice model ,sometimes called Kondo-Hubbard (KLU) model contains three terms. As in the non-interacting model there is the hopping term for conduction electrons c iσ (∼ t) and the Kondo exchange term (∼ J K ) to localized spins S i at every site. Their localization may be considered as result of their large repulsion U f . As mentioned before the Coulomb repulsion U for conduction electrons is usually neglected in the lattice model, although well studied for the Kondo impurity case. It is, however, very interesting to include its effect, both because it is physically present, in particular in intermetallic 3d-4f compounds and also because it allows to study theoretically the continuous crossover from the metallic Kondo screening case to the insulating case of interacting local spin dimers. The model is given by where n iσ = c † iσ c iσ . In 2D the first. n.n. hopping term leads to the single particle spectrum given by the dispersion ǫ k = −t(cos k x + cos k y ). In the continuum limit (N → ∞) it leads to a well known density of states (DOS) function shown in Fig. 1 by the full line which has a logarithmic singularity a ǫ = 0 corresponding to half filling. As was done in Ref. [15] for the uncorrelated KL model (U = 0) we will use the exact diagonalisation with Lanczos method for finite square lattice tiles. For the systematic classification of finite size tiles on the square lattice we refer to [22]. Due to the large number of eight states per site {S z = ± 1 2 } ⊗ {n = 0; n = 1(τ z = ± 1 2 ); n = 2} in the KLU model only tiles with eight sites can be used (see inset of Fig. 1). This precludes also finite size scaling such as may be performed in the 1D case [23,11]. Furthermore as in Ref. [15] we restrict to the half filled case with the number of conduction electrons given by n c /N = 1 (N=8). Due to small cluster size (N) n c can have only discrete values and therefore the next possible value n c /N = 3/4 in the singlet spin sector is already far from the interesting half filled case.
As discussed before [15] one may choose both open and periodic boundary conditions for the tiles (OBC and PBC respectively). Their single particle spectrum differs greatly. While for PBC most of the states are put directly at the Fermi energy, giving little resemblance to the continuum DOS the spectrum is more realistic for OBC with more even distribution of energies (Fig. 1). It was shown in the case U = 0 [15] that this leads to a qualitatively different moment screening or singlet formation for small J K . In this work we use exclusively OBC. Naturally, since there is a finite gap between occupied and unoccupied conduction electron states for half filling (Fig. 1), the true continuum behavior of the infinite Kondo lattice cannot be literally described, the numerical ED results presented in the following rather describe finite 'Kondo molecules'. As mentioned in the introduction they may in fact be artificially generated by adsorption of 3d clusters on nonmetallic surfaces. On the other hand we will focus mostly on the local moment and on-site or n.n. correlations. At least in the strong correlation limit (U/t ≫ 1) it turns out that the local quantities obtained numerically for the eight-site cluster agree rather well with the infinite lattice results (Sec. 5).

Numerical determination of the local moment and thermodynamic properties
To determine the ground state and thermodynamic properties of our system, we diagonalize the Hamiltonian iteratively on a small tile, applying the Lanczos algorithm. We also use the low-lying states and their energies to construct thermodynamic expectation values, in particular to calculate the partition function and its derivatives [24]. We stress that this unbiased method is exact and reliable for finite clusters and can therefore be used as a reference to compare with other numerical or analytical methods. We note that the large number of (eight) states per site prevents the use of larger clusters and the application of a finite size scaling procedure such as described in detail in Ref. [22] for a two-dimensional S = 1/2 model with just two states per site. Therefore using the Lanczos method for the present model is limited to tiles with eight sites (a maximum of 10 sites seems attainable). Even then the largest subspace with total spin S tot z = 0 has dimension 739 162 (38 165 260 for N = 10 sites). Nevertheless, as the comparison with analytical results will show, the local moment and on-site correlations for the small tile give a realistic approximation to the extended Kondo lattice behavior.
The total paramagnetic local moment is given by It contains the spin fluctuations of both the itinerant and localized spins and their antiferromagnetic on-site correlations due to the Kondo coupling. This is a central quantity for investigating the influence of correlations on the Kondo effect because the Coulomb repulsion U tends to localize the τ spins and thus influences strongly their singlet formation with S spins. For the uncorrelated case (U = 0) it was shown previously that the on-site Kondo coupling also induces effective RKKY inter-site coupling between the localized spins. This mechanism is still present for non-zero U but it has to compete with the effective superexchange. Therefore, in addition to the total moment we will also investigate the on-site singlet correlation S KS = τ i · S i and nearest neighbor (i, j) correlation function S ij = S i · S j . Due to the smallness of the cluster it is only reasonable to calculate it up to next nearest neighbors. The same would be true for the correlation function S i · τ j characterizing the 'screening cloud' of a given localized spin S i which is not evaluated here because its characteristic length scale is much larger than the cluster size. In the case of a single Kondo impurity in the continuum limit it extends over a range of ξ = v F /T K where v F is the Fermi velocity and T K the single ion Kondo temperature. Note however that even in this simple case the screening cloud and its length scale ξ have not yet been observed in reality [25].
With the FTLM approach [24] it is possible to calculate the temperature dependence of various important thermodynamic quantities like specific heat C(T ) = N A k BĈ (T ) and uniform susceptibility χ(T ) = N A µ 0 (gµ B ) 2χ (T ) per mole as cumulants of simple operators, e. g., where N A is the Avogadro number, µ 0 is the magnetic permeability, µ B and k B are the Bohr magneton and the Boltzmann constant, g is the gyromagnetic ratio, and N is the size of our tiles (or number of sites). For simplicity, we assume the same gyromagnetic ratio for both the τ and the S spins. The thermal expectation values are to be calculated with the statistical operator In this work we restrict to zero magnetic field, and as we are working with a finite system size N , we can safely set S tot z T = 0 in the first equation. Numerically, for each set of Hamiltonian parameters {t, U, J K }, the calculation of thermal averages involves two averaging procedures: Firstly, a random starting wave function is chosen, and the Lanczos algorithm yields up to O(100) extremal eigenvectors and eigenvalues. Secondly between 100 and 500 of these Lanczos runs with different starting wavefunctions are performed. All eigenvalues and eigenfunctions obtained in this way are subsequently used to calculate the traces over the statistical operator as indicated in Eqs. (3), leading to the desired thermal expectation values. For details we refer to Ref. [24].  Although we do not consider the magnetization M (per site) of the correlated Kondo model explicitly, we may get some information about finite-field properties by calculating the third-order susceptibility defined through the expansion for an applied magnetic field B. Here is given by a higher order cumulant according tô This quantity is a measure of the nonlinearity of magnetization at low field strength B. It has been discussed previously for a localized spin model [26] and plays a significant role in the discussion of some heavy fermion compounds [27].

Discussion of numerical results
First we discuss results for the total local paramagnetic moment presented in Fig. 2(a). Two counteracting trends determine its size: On one hand the increase in U localizes the τ -spins and leads to an increase of τ 2 i from 3/8 for U = 0 to 3/4 for U/t ≫ 1. On the other hand for finite J K the Kondo term establishes the AF on-site singlet correlation τ i · S i < 0. Therefore in the case J K =0 the moment increases monotonically from 9/8 = 1.125 to 2 · (3/4) = 1.5 while for any finite J K in the large U limit when τ i becomes localized the moment will decrease due to singlet formation. For moderate J K there is an initial increase of µ 2 loc with U due to the first correlation effect and eventually a decrease due to the effect of J K when τ spins become localized at larger U . In between a maximum in µ 2 loc develops as function of U . For larger J K the singlet formation effect is so strong that it overwhelms the increase in τ 2 i and therefore no maximum appears beyond J K /t=0.1. This effect is completely dominated by local correlations and should not be strongly influenced by the tile size. The corresponding local moment dependence on J K for various U is shown in Fig. 2(b). For the uncorrelated U = 0 case the singlet formation leads to a continuous decrease of µ 2 loc with J K . Increasing U facilitates this formation due to the localization of τ spins and the decrease becomes progressively steeper as function of J K . In the limit U → ∞ an arbibtrary small J K will lead to the singlet ground state.
The singlet formation may also be monitored directly by the on-site correlation τ i · S i which is shown in Fig. 3(a). Starting from zero at J K = 0 it becomes increasingly antiferromagnetic for growing J K until the singlet value −3/4 is reached. Again the latter is approached more rapidly with increasing correlation strength U . The nearest neighbor induced AF spin correlations are presented in Fig. 3(b). For U=0 these correlations are of the induced RKKY type [15]. They decrease with increasing J K because their evolution is impeded by the increasing singlet formation. For a fixed but small J K the antiferromagnetic inter-site correlation first increases for small U due to the reduction of single and double occupancies and for larger U it falls off again due to the reduction of the superexchange 4t 2 /U with increasing U. For larger J K this effect is more pronounced as function of U.
A recurrent topic in previous analytical theories of the correlated Kondo impurity model is the U-dependence of the Kondo energy scale. From a practical viewpoint it is frequently taken as the maximum position T * of magnetic specific heat or susceptibility which corresponds to an average singlet-triplet excitation energy. In the present context of finite size tiles the U or J K dependence of the characteristic temperature scale T * can be conveniently obtained from the FTLM results according to Eq. (3) for both quantities. Although it is not identical to the single ion Kondo temperature T K of the impurity model one may expect a comparable qualitative dependence on U or J K may be comparable.
The contour plot of C V (T ) in the T-U plane is shown in Fig. 4(a),(b) for J K /t = 0, 2 respectively. The former corresponds to the case of the pure Hubbard model. It illustrates that the high temperature peak at U = 0 due to uncorrelated charge fluctuations splits into a lower temperature spin fluctuation peak associated with J = 4t 2 /U and a broad continuum at very high temperatures. The charge fluctuation peak at U=0 correspond to the excitation energy of ≈ 0.75t between the highest occupied and lowest unoccupied states around ǫ = 0 in Fig. 1. The behavior as function of U may be viewed as a precursor of the thermodynamic Mott-Hubbard transition for a finite tile size. When J K is turned on the specific heat peak is dominated by the lowest singlet triplet excitation energies. It increases linearly for small U (Fig. 4(b)) and reaches a plateau around U ≈ 6t corresponding to the highest excitation energies in Fig. 1.
To get rid of the influence of charge fluctuations we also calculated the spin susceptibility with FTLM. The results are presented in Figs. 5(a),6 as function of J K , U respectively. In the former a contour plot shows the evolution of susceptibility maximum or Kondo temperature scale T * (black dots) as function of J K for constant U/t = 4. For comparison the T * maxima for the uncorrelated U = 0 Kondo lattice tile are also given in Fig. 5(b) by the crosses. In both cases T * increases linearly with J K in the strong coupling limit J K /t ≫ 1. For small J K the values for U = 0 are considerably below the results for finite U which indicates the increase of the Kondo temperature scale T * with U. This is also seen directly in Fig. 6 where T * (black dots) can be seen to increase linearly with U for small U and reach a plateau for U > 6t similar as for the specific heat peak before. In that figure we also included the peak position of the third order susceptibilityχ (3) (T ) (white dots) which also increases with U. This quantity is experimentally accessible [27]. It peaks at a systematically lower temperature than the first order susceptibility. The reason is that it characterizes the strongest nonlinear increase of the local moment with field which happens precisely in the temperature region where the screened local moment and susceptibility at zero field drops to zero (see Fig. 6). For comparison the size of the spin gap from analytical calculation in the large U limit (Eq. (13)) shown by white crosses is seen to be sandwiched between the above FTLM values for the first order (black dots) and third order (white dots) susceptibility peak positions.

Bond operator treatment of the strongly correlated 'Kondo necklace' limit
A deeper insight into the phases and excitations of low dimensional quantum magnets requires the application of both numerical and analytical techniques. As demonstrated e.g. by the spin wave analysis of the frustrated 2D spin systems ( Ref. [22] and references cited therein) analytical results for the extended system, even if approximate or only available in limiting cases , are very helpful to understand the systematics of numerical ED results for finite clusters. It is therefore perfectly legitimate to proceed in a similar way for Kondo lattice models. Because we focus on the paramagnetic phase we will, however, use the bond operator approach as the appropriate analytical technique for comparison in the large U limit. In the following we will derive such analytical results in the limit of large conduction electron correlations (U ≫ 2zt, z=4 is the coordination number) at half filling (n c = 1). In this case the charge fluctuations may be eliminated from the hopping and Hubbard terms leading to a pure exchange term of the now also localized conduction electron spins τ i . Thus the appropriate large U limit of the model in Eq. (1) is a localized 2D spin Hamiltonian [28] given by Here the kinetic hopping term in Eq. (1) is replaced by an effective inter-site spin exchange J = 4t 2 U in the strong correlation limit of conduction electrons. The above Hamiltonian is of the generalized 'Kondo-necklace' type (in 2D) originally studied for 1D chains with only the xy inter-site terms included [29]. It has, however, later been extended to higher dimension and including all components in the intra-and intersite exchange with possible uniaxial anisotropies for both [17,18,19]. The model may also be viewed as an asymmetric bilayer Heisenberg model [30] with τ i and S i spins residing on different layers and only the former coupled by inter-site exchange J. These are generic models to describe quantum phase transitions between a singlet ('Kondo') phase favored by the second term and an antiferromagnetically ordered phase favored by the first term. The transition occurs when the control parameter J/J K is larger than a the value (J/J K ) c defining the quantum critical point (QCP). In 2D we have in the present isotropic model (J/J K ) c = 0.88 [18,19]. Such transitions are frequently found in f-electron compounds where the control parameter may be varied by pressure or doping (i.e., chemical pressure). The above model allows to study the characteristic quantum critical behavior around the QCP disregarding the charge fluctuations. It has been investigated using numerical methods like Monte Carlo (MC) simulations [31,32], exact diagonalization methods [14] and dynamical mean field theory (DMFT) and also analytical methods like bond operator approach in mean field [17,18,19] or hard-core boson treatment [33]. In this approach the Kondo necklace Hamiltonian in Eq. (6) is mapped to a model of interacting singlet (s) and triplet (t α , α = x, y, z) bosons by the bond-operator transformation [34]. These bosons describe the singlet and triplet states |s = s † |0 and |t α = t † α |0 (α = x, y, z) of the pair of spins (τ i , S i ) coupled by J K at every site i. In terms of bosonic operators the spins are given by where α, β, γ = x, y, z and ǫ αβγ is the totally antisymmetric tensor. The conventional spin and boson commutation rules are fulfilled. The restriction to physical states (only one boson per site) is expressed by the constraint s † i s i + α t † i,α t i,α = 1 which may be implemented either on the mean field level or by hard core boson technique. The former is chosen here.

Ground state energy and triplon excitations
In mean field approximation the bond operator transformation leads to a bilinear bosonic Hamiltonian where µ is the chemical potential introduced to ensure the constraint ands = s is the mean field singlet amplitude. Furthermore we define t kα = 1 √ N i exp(ikR n )t iα and γ k = (cos k x + cos k y ). Here we restrict to the nonmagnetic phase where the triplet amplitudet ≡ 0. In this cases is always close to one. This Hamiltonian may be diagonalised by the Bogoliubov transformation where ω k is the threefold degenerate (α = x, y, z) triplon dispersion and E 0 (µ,s) the ground state energy. Minimization of the latter leads to selfconsistency equations for the chemical potential µ and singlet amplitudes given by 3 2N

Diagonalization leads to the bosonic triplon Hamiltonian
The smooth dependence of µ,s on the control parameter J K /J is shown in the inset of Fig. 7. It also extends continuously across the QCP into the antiferromagnetic region [19]. The bond operator method also gives an explicit expression for the singlet-triplet gap which is defined as ∆ t = ω Q with Q = (π, π) [18]. Using Eqs. (10,11) and setting J = 4t 2 /U we obtain This expression for the gap may be compared to the Kondo temperature scale T * from the susceptibility maximum obtained in FTLM finite cluster calculations.

The paramagnetic effective local moment
In the numerical calculation the central quantity in the correlated Kondo model is the local paramagnetic moment µ 2 loc which is a direct measure of the singlet formation as function of J K /t and U/t. This also extends to the large U limit where the local moment may be calculated analytically within bond operator approach. To gain a better understanding of the numerical results a comparison with the analytical method for large U is helpful. In spirit this is similar to the comparison of ED results for the J 1 − J 2 Heisenberg model with analytical spin wave calculations [22]. The elementary excitations here are, however, gapped triplon modes in the paramagnetic ('Kondo singlet ') regime rather than spin waves in the antiferromagnetic broken symmetry state.
First we transform the expression for µ 2 loc given in Eq.
(2) to bond operator basis. Using the defining relations in Eq. (7) with θ α = 1 2 (t † × t) α we derive the operator identity where we used the spin space isotropy or degeneracy of triplon modes for the expectation value. Using the explicit expressions for θ α , defining u k = cosh(φ k ), v k = sinh(φ k ) and performing a Fourier transformation we obtain a concise form of the moment: The positiveness of the moment squared is ensured by Schwartz' inequality. Using Eq. (10) we can write Inserting these explicit expressions into Eq. (15) we obtain the final result for the paramagnetic moment as where d k , f k are defined in Eq. (10). From this closed expression the moment may be obtained by the 2D momentum integration. It is worthwhile to consider the expression of µ 2 loc in the strong Kondo coupling limit J/J K ≪ 1 or (4t 2 /U J K ≪ 1). We define ǫ k = 2f k d k and use ω k = d k 1 − ǫ 2 k . Expanding Eq. (2) in terms of ǫ k we get Which may be further evaluated to give Where in the last expression we replaced J = 4t 2 U . This equation demonstrates that the local moment is decreasing with increasing J K (note that µ ands also depend on J K /J ) and with increasing correlation U. This is what the full numerical calculation of Eq. (17) discussed below indeed confirms.

Spin correlations and high temperature susceptibility
To obtain a more detailed insight into the Kondo singlet formation and the influence of Coulomb correlations on it we have previously also calculated the evolution of on-site and next-neighbor spin correlations in the ED approach. It is desirable to calculate them with bond operator approximation in large U limit for comparison with numerical results. First we consider the on-site spin correlation function S KS (i) = τ i · S i between localized and itinerant spins and also the complementary partial local moments τ 2 i and S 2 i . Transforming to bond operator representation and using the isotropy we obtain The equality τ 2 i = S 2 i is only valid in the localized limit but does not hold in the original KLU model defined in Eq. (1) for small U. The site index i is suppressed i.f. because these local quantities are uniform. Expressing the triplet operators t x , t † x in terms of triplon eigenmodes by using Eq. (9) this leads to from Eq. (16) we finally obtain If we take the sum of these expression according to Eq. (2) one obtains again the total local moment. The total moment in the large-U (localized spin) limit calculated in Sec. 5.2 is a zero temperature quantity determined by quantum fluctuations in the ground state. On the other hand the effective moment of a localized spin system is obtained from the high temperature behavior of the susceptibility. The latter may be calculated from a high temperature expansion. In this section we want to investigate how these effective moments are related and how well the high temperature expansion agrees with the FTLM results in the large-U limit. First we give the high temperature expansion for the total local moment at finite temperature T : where is the thermal expectation value of A formed with the Kondo necklace Hamiltonian. Expanding the statistical operator for large k B T we obtain The prefactor is the paramagnetic moment of uncorrelated spins, where the factor 2 is due to the presence of two spins τ and S per site with τ = S = 1/2. The parenthesis gives the first order high temperature correction of the local moment which depends only on the local Kondo exchange J K . Likewise we can calculate the high temperature uniform susceptibility, which was defined previously aŝ where the same g-factor has been assumed for both τ -and S-spins. Evaluating the thermal average in high temperature approximation as before and using τ = S one obtains for the high temperature susceptibilitŷ 2S(S + 1) 1 − S(S + 1) 3 Formally, we can regard the last expression as the first two terms of a high-temperature expansion of a Curie-Weiss type susceptibilitŷ with an effective noninteracting Curie susceptibilityχ 0 (T ) and a Weiss temperature Θ given, respectively bŷ i.e. k B Θ = J for S = 1 2 , z = 4. We write the high-temperature behavior of the susceptibilityχ(T ) in this suggestive form to stress that the temperaturedependent effective moment screening is due to the Kondo coupling J K , and the Weiss temperature Θ due to the effective magnetic intersite exchange J. One may conjecture that this expression may be valid beyond its formal expansion regime. This may be checked by comparing Eq. (26) with the unbiased FTLM results forχ(T ) andμ 2 loc (T ).

Comparison with numerical results
It is an interesting question whether the previous analytical mean field results for the large U limit can be qualitatively compared to the numerical ED results for finite tiles in Sec. 4. The most obvious quantity to check this is µ 2 loc . The comparison is shown in Fig. 7. The inset of this figure gives the dependence of µ ands on J/J K as obtained from the self-consistent solution of Eq. (12) almost up to the quantum critical value (J/J K ) c = 0.88 where the singlet-triplet spin gap vanishes and AF order would set in. We stay in the paramagnetic parameter range throughout this work. The AF region of the 2D Kondo necklace region has been explored in Refs. [18,19].
The main Fig. 7 shows the dependence of the local moment µ 2 loc 1 2 and its square on J/J K from bond-operator mean field theory (full black line) and from numerical results for the eight site cluster model in the large U limit with J = 4t 2 /U . Here large U means large as compared to the tight binding bandwidth W = 2zt = 8t. The numerical results for various U/t in this limit are given as dashed lines. The agreement of analytical and numerical results is surprisingly good. This proves two points: Firstly the mean-field bond operator technique gives reliable on-site properties in the ground state. This is in agreement with the observation [19] that values for (J/J K ) c that determine the quantum critical point between Kondo singlet and AF phase are well reproduced by that theory. Secondly the agreement points to the fact that the numerical finite size effects on a local quantity like µ 2 loc are apparently quite moderate. Noticeable deviations in the numerical and analytical results appear at larger J/J K when the quantum critical point to AF order is approached or when U/t becomes too small. The comparison may also be made for the on-site spin correlation functions S KS (i) and the partial moments τ 2 and S 2 which are equal in the large U limit. The latter is shown in Fig. 8a with numerical results from ED giving the proper local moment value 3/4 and the results from bond operator theory lying below. Their difference increases with increasing J, i.e. decreasing U . On the other hand the onsite AF Kondo correlation presented in Fig. 8b shows the analytical result lying above the numerical ED value by a similar amount. The latter depends on the effective coordination number of the considered site in the finite cluster. The total moment µ 2 loc is the sum of these individual contributions and because of the opposite sign of the differences the add up to zero approximately. This explains why numerical ED and analytical bond operator results for µ 2 loc in Fig. 7 agree so well despite the small cluster size. Including the previous results [18,19] on quantum critical properties of the Kondo necklace model we can conclude that the energetics and local correlations are well described by the mean field bond operator method. However, it should be noted that the description of inter-site correlations and their U dependence is well beyond this approximation.
Finally we come to the high temperature expansion for local moment and susceptibility in the large U limit and its comparison to the FTLM results. According to Eq. (26) the ratio T χ(T )/μ 2 loc (T ) should be proportional to the Curie Weiss factor in this equation. This comparison is shown in Fig. 9. It is seen that the result of the high temperature expansion lies considerably below the curve obtained from FTLM. Part of this discrepancy may be due to the fact that in FTLM the effective (average) coordination number in the eight site tile is smaller than z = 4. Therefore, in order to compare with high temperature expansion we may consider an effective Curie Weiss temperature Θ eff as fit parameter in Eq. (26). Then the temperature dependence of χ(T ) from FTLM is well reproduced by the form of the high temperature expansion results in Fig. 9.

Summary and conclusion
In this work we investigated the local moment screening and spin correlations and in particular thermodynamic properties of the correlated Kondo lattice or Kondo-Hubbard model. We have focused primarily on the systematic dependence of local moment, spin correlations and susceptibility on the control parameters J K and U and secondly on the U-dependence of the characteristic temperature scale T * . We used unbiased and exact numerical techniques like ED and FTLM for small tiles for the whole range of U/t and J K /t as well as approximate analytical bond-operator method for the large U-limit of the extended Kondo necklace model which contains only localized spins. In this limit both methods give excellent agreement on the dependence of the Kondo-screened total paramagnetic local moment (Fig. 7) on U or J = 4t 2 /U . Our ED and FTLM investigation may also be of particular relevance for finite nano-clusters of Kondo atoms adsorbed on surfaces.
A first central result obtained for smaller U and J K is a clear non-monotonic behaviour of µ 2 loc on correlation strength is observed in the ED results. This nonmonotonic U-dependence is also observed in the on-site and inter-site spin correlations where the latter are of mixed RKKY and superexchange character. It is the result of competing effects of conduction electron localization by U and local singlet formation due to J K .
The second central conclusion concerns the dependence of the Kondo temperature scale T * on correlation strength U for which controversial results have been reported previously. Our FTLM susceptibility and the analytical results presented in Fig. 6 show that it increases monotonically with U and reaches a plateau in the large U limit. In this limit the Kondo scale corresponds to the spin-gap for triplon excitations at the AF wave vector Q = (π, π).
Furthermore a high temperature expansion for the susceptibility of the Kondo necklace model leads to an effective Curie-Weiss type expression where the local susceptibility is modified only by the Kondo term and the Curie-Weiss temperature is only due to the superexchange term. The resulting temperature dependence is similar to that of FTLM results in the large U limit and if the effective Curie-Weiss temperature is used as a fit parameter a quantitative agreement over large temperature region is obtained.