# Mutual information in heavy fermions systems

###### Abstract

A key notion in heavy fermion systems is the entanglement between conduction electrons and localized spin degrees of freedom. To study these systems from this point of view, we compute the mutual information in a ferromagnetic and anti-ferromagnetic Kondo lattice model in the presence of geometrical frustration. Here the interplay between the Kondo effect, the RKKY interaction and geometrical frustration leads to partial Kondo screened, conventional Kondo insulating and antiferromagnetic phases. In each of these states the mutual information follows an area law, the coefficient of which shows sharp crossovers (on our finite lattices) across phase transitions. Deep in the respective phases, the area law coefficient can be understood in terms of simple direct product wave functions thereby yielding an accurate measure of the entanglement in each phase. The above mentioned results stem from approximation-free auxiliary field quantum Monte Carlo simulations.

Introduction.— The Kondo effect describes the screening of a spin-1/2 magnetic impurity embedded in a metallic environment Hewson. At high temperatures the spin degree of freedom is decoupled from the conduction electrons and below the Kondo scale a many body entangled state of the spin and conduction electrons emerges. To quantify entanglement between a bipartition A and B of a system, one traces out the degrees of freedom B to obtain a reduced density matrix, , the Renyi entropy of which, , corresponds to the entanglement entropy. Taking one subsystem to be the impurity spin, and the other the conduction electrons, the Kondo effect can be elegantly characterized by the transfer of thermal entropy at high temperatures to entanglement entropy in the ground state Sorensen07. The energy scale at which this transfer from thermal to entanglement entropy occurs is the Kondo temperature.

In the presence of a lattice of spins Kondo coupled to conduction electrons, corresponding to heavy fermion systems Coleman07_rev, the above picture breaks down. In fact spins can now interact through the indirect Ruderman-Kittel-Kasuya-Yosida (RKKY) exchange interaction RKKY, and thereby compete with Kondo screening. Comparing energy scales, it becomes apparent that Kondo screening dominates when the exchange interaction between the local spins and the conduction electrons, , is positive and large, and that the RKKY interaction dominates at small values of . The intricate interplay between these two effects on non-frustrated lattices leads to quantum phase transition (QPT) between disordered and ordered magnetic phases Doniach77; Assaad99a. The Kondo effect can be switched off by considering , thereby promoting magnetically ordered phases Assaad-KLQMC. In addition, geometrical frustration is found to be of experimental relevance in many heavy fermion materials such as CePdAl, PrIrO, YbAgGe, YbAlC, YbPtPb akito16; nakatsuji06; kim_frustkondo08; sengupta_frustkondo10; kato_frustkondo08, where quantum phases do not easily fit into aforementioned cases. Geometrical frustration can lead to so called partial Kondo screened (PKS) phases where frustration is alleviated by selective spacial screening local spins SAG-18; Aulbach2015. The essence of all aforementioned states can be captured by direct product variational wave functions from which one can directly assess the degree of entanglement between the spins and conduction electrons. Entanglement entropies, although not presently experimentally accessible in heavy fermion systems, lend themselves to an experimental measure in systems of cold atoms IMPTLRM-15, which, in turns, allow to realize Kondo lattice models GHGXJYZDLR-10. Alternatively, entanglement properties have been recently proposed to be experimentally studied by engineering the so-called entanglement Hamiltonian in cold atom systems DVZ-18. In this context, Ref. PTA-18 introduces a numerically exact method to determine the entanglement Hamiltonian in interacting models of fermions.

In this letter we investigate a Kondo lattice model Hamiltonian amenable to negative sign free quantum Monte Carlo (QMC) simulations that provides specific realizations of the states discussed above. Using recently developed methods to compute the Renyi entropies Grover13 with the auxiliary field QMC, we compute the mutual information and show that, deep in the respective phases, the numerical value of the area law coefficient can be well understood in terms of the product state wave function description of the phase supplemented by fluctuations, if necessary. Furthermore, we observe a singular behavior of the area law coefficient across phase transitions.

Model.— We consider the generalized Kondo model on the honeycomb lattice introduced in SAG-18 with Hamiltonian

(1) |

where the first sum extends over the nearest-neighbor sites and describes the hopping of conduction electrons, , giving rise to the well-known semimetallic band dispersion CNGPNG-09, the second sum accounts for the Kondo screening interaction between c-electrons and spin 1/2 local moments , while in the third term the next-nearest-neighbor antiferromagnetic interaction between localized spins encodes frustration effects. This model can be solved without encountering the negative sign problem SAG-18; here and in the following we fix . Due to the antiferromagnetic coupling this half-filled Kondo lattice model on the honeycomb lattice with geometric frustration exhibits novel PKS phases alongside the conventional Kondo insulator (KI) and antiferromagnetically ordered phases SAG-18.

For a given subsystem of conduction electrons and of local spins , the mutual information between and is

(2) |

where is the -th Renyi entropy for a subsystem . Here we take as a compact subset of conduction electron sites, and as the corresponding local spin sites coupled to the subset . Assuming the ubiquitous area law for the entanglement entropy ECP-10, results are proportional to the size of the boundary shared between and :

(3) |

The mutual information is also defined in terms of the von Neumann entanglement entropy, corresponding to the limit in Eq. (2). In this case satisfies OhyaPetzBook

(4) |

where the numerator represents the connected correlation of two arbitrary operators and acting on the subsystem and , respectively, and is the norm of an operator . According to Eq. (4), bounds all mutual correlations of operators in and , thus providing an operator-independent entanglement measure. Due to the above bound, captures both high and low energy scales.

Here, we shall consider the mutual information for Renyi index . The coefficient introduced in Eq. (3) is the main quantity investigated in this work. In the presence of Kondo screening essentially counts the number of Kondo singlets formed between and , such that in the limit Eq. (3) holds exactly with .

Method.— We have investigated the Hamiltonian of Eq. (1) by means of auxiliary field QMC Blankenbecler81; White89; AF_notes simulations, using the method of Ref. SAG-18 which, in essence, consists in a fermion representation of localized spins obtained via Lagrange multipliers. We refer to Ref. SAG-18 for more details on the formulation. A similar technique can be used to simulate the canonical ensemble WAPT-17. Simulations have been performed using the ALF package ALF. To compute the Renyi entropies we have used a method introduced in Ref. Grover13, and also used in Grover13; ALPT-13; DP-15; DP-16, which allows to formulate the reduced density matrix within auxiliary field QMC. Beside Renyi entropies, the technique can be exploited to unbiasly determine the entanglement Hamiltonian PTA-18. Ref. PTA-18b provides a short review of computational approaches to entanglement in interacting fermionic systems.

Results.— In Fig. 1 we reproduce the rich phase diagram of the model for , . At finite the model has a reduced U(1) spin symmetry corresponding to spin rotations around the -axis, as well the point group and translation symmetries of the Honeycomb lattice. The Kondo insulating (KI) phase breaks no symmetries, the antiferromagnetic (-AFM) phase breaks the U(1) spin symmetry and the -PKS breaks nematically the point group and has reduced translation symmetry. The -PKS differs from the -PKS in that it additionally breaks the U(1) spin symmetry SAG-18. In Fig. 1 we also show three lines on the phase diagram where we analyze the mutual information. Moreover, we study the entanglement for , , where the model favors the formation of an effective spin Heisenberg antiferromagnet in the strong coupling limit. In all the QMC data presented here we have simulated unit cells corresponding to 144 orbitals.

We first consider a scan for , where the model reduces to a standard Kondo lattice model on the honeycomb lattice with the conventional RKKY driven AFM phase and KI phase. A fit of for up to choices of , to the right-hand side of Eq. (3) allows to extract the coefficient shown in Fig. 3, for two inverse temperatures , , and as a function of the minimum subsystem size taken into account in the fits. More technical details on the chosen subsystems , are reported in the Supplemental Material SM. We observe consistent results for the fitted values of . In the KI phase conduction and localized electrons are paired into a spin singlet, such that for the ground state approaches a product of single-site singlet wave functions shown in Fig. 2(c), giving a entanglement entropy per pair. Indeed, in the KI phase reaches an asymptotic value for . Interestingly, a change of concavity in the plot of occurs around the QPT between the AFM and KI phases, at SAG-18. In the AFM phase at the conduction electron local moment aligns antiparallel to the spin and is reduced in magnitude due to charge fluctuations. The essential features of the AFM phase can be captured by a product wave function, in which the total local moment of conduction electrons and spins form a Néel order as depicted in Fig. 2(a). Entanglement between spins and conduction electrons originates from sub-leading Kondo screening as argued in Ref. Assaad-KLQMC and also from long wave length spin wave fluctuations of the Néel order parameter. In Fig. 3 we also show for a ferromagnetic coupling . Different than the case, here the Kondo coupling favors the formation of a spin triplet in the ground state of the model, without screening. Starting from the limit , where on each lattice site the spin singlet state is projected away, a finite large value of gives rise to an antiferromagnetic exchange term between states on the honeycomb lattice. Thus, in this situation the system reduces to a Heisenberg model. Its antiferromagnetic ground state is well captured by the semiclassical large- expansion, and consists in a Néel state, illustrated in Fig. 2(b). Since such a ground state is primarily built on states, the entanglement between conduction electrons and localized spins is expected to be small. This observation is clearly reflected in the coefficient shown in Fig. 3, whose values for are substantially smaller than for . For instance, we find for and for . The computed value of grows on reducing , but appears to saturate to a value significantly smaller than the limiting value found for .

In order to investigate the -PKS phase we compute the mutual information as a function of for fixed . As illustrated in Fig. 1, this second scan crosses the conventional KI phase and the -PKS phase. Due to an enlarged unit cell expected in the -PKS phase SAG-18, we have in this case considered only three possible subsystems , , with equal size , , . In view of the limited amount of available data, and in order to reliably study the coefficient we have considered three possible linear fits: a fit including all data, a fit disregarding the smallest size , and a fit disregarding the largest size . In Fig. 4(a) we show the corresponding results, for two inverse temperatures , . Despite fluctuations larger than the error bars, indicating sizable corrections to Eq. (3), we observe a clear trend in , which decreases from to . Moreover, the curve shows again a change in the curvature at a value of approximately consistent with the onset of a QPT between the KI phase and the -PKS phase, located at SAG-18 (see also Fig. 1). For , saturates to a plateau . Such a value, which approaches the limit at fixed is in fact -dependent, as shown by the computation of along the path (3) on Fig. 1. In Fig. 4(b) we show the resulting , whose curves are qualitatively similar to the case of Fig. 4(a), but saturate to a large value for large .

The emergence of an area law with a coefficient confirms the mechanism of partial screening in the -PKS phase, emerging from the competition of the antiferromagnetic interaction and the Kondo coupling. To understand the structure of the ground state we consider the limit , , i.e., the atomic limit in which the honeycomb lattice decomposes into two independent triangular sublattices. In this limit the Hamiltonian (1) commutes with the -component of the total spin operator on the site , , with . Accordingly, one expects the ground state to have , such that its wave function is constructed from the states and , whereas states obtained with the remaining base vectors and are gapped. In the Hilbert space spanned by the Hamiltonian (1) is, up to a constant,

(5) |

where the operators , are defined by , and satisfy the commutation rule of the SU(2) algebra. Therefore, in the atomic limit and close to the ground state the Hamiltonian (1) reduces to that of an antiferromagnetic transverse-field Ising model, on a triangular lattice. Beyond the atomic limit, the effective low-temperature Hamiltonian acquires additional interactions which we have already anticipated in Eq. (5). However since for a finite, but large, , and still the additional states are gapped, we expect the ground state to be well representable in the basis. It is know that the above Ising model has a three sublattice structure and that the sign of the 6-fold clock term in the Landau Ginzburg functional determines if the ground state will be hierarchical, , or uniform Blankschtein84; Moessner01b,

(6) |

Here, , , label the three sites of the sublattice structure. Numerical calculations in Ref. SAG-18 point to the uniform ground state, so we use this product wave function to account for our mutual information results. Here, controls the partial polarization of the sites and and the direct product wave function is illustrated in Fig. 2(d). As we discuss in Ref. SM, for the state of Eq. (6), , which takes values between and , thus including the plateaus found in Fig. 4.

Conclusions and Discussion.— The bound on the von Neumann entanglement entropy based mutual information presented in Eq. 4 implies that for our choice of and the mutual information picks up IR as well as UV physics. Deep in a phase, or in a sink in the Renormalization Group parlance, where the correlation length is finite one can model the ground state with a direct product wave function that provides a modeling of all energy scales. With this wave function the coefficient of the area law of the mutual information can be computed. In the Kondo lattice model considered here, we have an explicit realization of three phases, AFM, PKS and KI. The simple direct product wave functions for these three phases presented in Fig. 2 have an area law coefficient set by for the Néel representation of the AFM phase, of for the strong coupling KI and of bounded by for the PKS phase. The KI value of is exact in the strong limit and is very well reproduced when is comparable to half the bandwidth, . Deep in the PKS phase shows a dependent plateau upon enhancing . The value of this plateau depends on the degree of partial Kondo screening as described by in Eq. 6. Clearly, the Néel wave function shows no entanglement, but of course does not capture fluctuations leading to Goldstone modes and to entanglement. For these phases, both for positive and negative values of the mutual information remains small but does not vanish. It is however noticeable that is larger in the AFM phase at than at . We understand this difference in terms of sub-leading Kondo screening present for the antiferromagnetic model but absent for the ferromagnetic one. Computations of the single particle spectral function Assaad-KLQMC confirm this point of view.

On our finite lattices we have observed clear crossovers in the area law coefficient of the mutual information across phase transitions. How this behavior reflects the criticality of the transition as well as possible corrections to the area law at critical points is left to future investigations.

From the technical point of view, the present calculation of the mutual information does not lead to a noticeable increase in computational effort. It can used as a standard observable independent measure to pick up QPTs and thereby map out phase diagrams in various heavy fermion systems and beyond. This is very similar to unsupervised machine learning algorithms aiming at automatically mapping out phase diagrams Broecker17. As mentioned above it also provides further information on phases. In this context we note that this quantity has recently been used Hofmann18 to validate the understanding of a Kondo breakdown transition.

Acknowledgments.— We would like to thank T. Grover and J. S. Hofmann for insightful discussions. F.P.T. thanks the German Research Foundation (DFG) through Grant No. AS120/13-1 of the FOR 1807. T.S. thanks the DFG for financial support from grant number AS120/14-1. F.F.A. thanks the DFG through SFB 1170 ToCoTronics. We gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for allocation of CPU time on the SuperMUC computer at the Leibniz Supercomputing Center.

bu1.aux\@input@bu1.bbl\@inputbu1.aux\@inputbu.aux

## I Supplemental Material

## Ii Computation of the mutual information

In order to compute the coefficient introduced in Eq. (3) we have sampled choosing different subsets of lattice sites. As explained in the main text, is a compact subset of lattice sites of the conduction electrons and contains the local spin sites coupled to . In Fig. 5 we illustrate the choices of , used here. In the presence of nonvanishing antiferromagnetic coupling , due to the enlarged unit cell SAG-18 we use only the subsystems No. , , and , shown in Fig. 5. As we discuss in the main text, for a given value of and , the coefficient has been obtained by a linear fit of as a function of the size . In Figs. 7, 7 we show two examples of such a procedure, in the KI phase and in the -PKS phase.

## Iii Renyi entropies in the -PKS phase

For a product wave function Ansatz like that in Eq. (6), and for any choice of and , consisting in a set of conduction electron sites and in the set of the corresponding coupled local spin sites, the reduced density matrix of is a pure state. Hence, and , such that, in order to determine is it sufficient to compute the entanglement entropy of the localized spins . On each triangular unit cell, the Ansatz of Eq. (6) results in a factorized reduced density matrix for

(7) |

with as given in Eq. (6) and

(8) |

Since, apart from a site index, the single-site reduced density matrices and are identical, in Eq. (8) and in the following, with a slight abuse of notation, we have dropped the site indexes , in . Using the definition and , the reduced density matrix for the sites and is readily computed as

(9) |

Using Eq. (8) and Eq. (9) the second Renyi entropy for the single triangular unit cell is

(10) |

Finally, the coefficient defined in Eq. (3) corresponding to the Renyi entropy of Eq. (10) is

(11) |

In Fig. 8 we plot as a function of , for . The minimum value of is , and corresponds to , : for such values of the state is a product of a conduction electrons and a localized spins states. Viceversa, a maximum value of is found for , corresponding to the maximally entangled spin singlet state.

bu2.aux\@input@bu2.bbl\@inputbu2.aux\@inputbu.aux