Matrix Product Study of Spin Fractionalization in the 1D Kondo Insulator

The Kondo lattice is one of the classic examples of strongly correlated electronic systems. We conduct a controlled study of the Kondo lattice in one dimension, highlighting the role of excitations created by the composite fermion operator. Using time-dependent matrix-product-state methods we compute various correlation functions and contrast them with both large-N mean-field theory and the strong-coupling expansion. We show that the composite fermion operator creates long-lived, charge-e and spin-1/2 excitations, which cover the low-lying single-particle excitation spectrum of the system. Furthermore, spin excitations can be thought to be composed of such fractionalized quasi-particles with a residual interaction which tend to disappear at weak Kondo coupling.


I. INTRODUCTION
Kondo insulators are an important class of quantum material, which historically, foreshadowed the discovery of heavy fermion metals and superconductors [1].These materials contain localized d or f-electrons, forming a lattice of local moments, immersed in the sea of conduction electrons [2][3][4][5].Remarkably, even though the high temperature physics is that of a metallic half-filled band, at low temperatures, these materials transition from local moment metals, to paramagnetic insulators.In the 1970s, theorists came to appreciate that the origin of this behavior derives from the formation of local singlets through the action of an antiferromagnetic exchange interaction between electrons and magnetic moments [2,3,6,7], a model known as the Kondo lattice Hamiltonian.
The Kondo lattice model contains a tight-binding model of mobile electrons coupled antiferromagnetically to a lattice of local moments via a Kondo coupling constant J.The deceptive simplicity of this model hides many challenges.Perturbative expansion in J, reveals that the Kondo coupling is marginally relevant, scaling to strong-coupling at an energy scale of the order of Kondo temperature T K ∼ W e −1/Jρ .Moreover, the localized moments, with a two-dimensional Hilbert space, do not allow a traditional Wick expansion of the Hamiltonian, impeding the application of a conventional field-theoretic methods.
The strong-coupling limit of this Hamiltonian, in which J is much larger than the band-width, J/t c 1 provides a useful caricature of the Kondo insulator as an insulating lattice of local singlets.In the 1980s [8][9][10][11], new insight into the Kondo lattice was obtained from the large-N expansion.Here, extending the spin symmetry from the SU(2) group, with two-fold spin degeneracy, to a family of models with N fold spin-degeneracy allows for an expansion around the large-N limit in powers of 1/N .The physical picture which emerges from the large-N expansion accounts for the insulating behavior in terms of a fractionalization of the local moments into spin-1/2 excitations, S j → f † jα ( σ/2) αβ f jβ which hybridize with conduction electrons [7,[9][10][11][12] to form a narrow gap insulator.However, the use of the large-N limit provides no guarantee that the main conclusions apply to the most physically interesting case of N = 2.
In this paper we use matrix product state methods to examine the physics of the one dimensional Kondo insulator.Our work is motivated by a desire to explore and contrast the predictions of the strong coupling and large-N descriptions with a computational experiment, taking into account the following considerations: • Traditionally, Kondo insulators are regarded as an adiabatic evolution of a band-insulating groundstate of a half-filled Anderson lattice model.We seek to understand the insulating behavior, which is akin to a "large Fermi surface", from a purely Kondo lattice perspective, without any assumptions as to the electronic origin of the local moments.
• What are the important differences between the excitations of a half-filled Kondo insulator and a conventional band insulator?
• Many aspects of the Kondo lattice suggested by the large-N expansion, most notably the formation of composite fermions and the associated fractionalization of the spins, have not been extensively examined in computational work.In this respect, arXiv:2302.09701v1[cond-mat.str-el]20 Feb 2023 our work complements the recent study of Ref. [13], highlighting the mutual independence of the conduction electron and composite fermions through the matrix structure of the electronic Green's function.We extend this picture even further by examining the dynamical spin susceptibility, providing evidence for fractionalization of the spin into a continuum of quasi-particle excitations.

A. Past studies
Our work builds on an extensive body of earlier studies of the 1D Kondo lattice that we now briefly review.The ground-state phase diagram of this model was first established by Tsunetsugu, Sigrist, and Ueda [14], who established the stability of the insulating phase for all ratios of J/t c , while also demonstrating that upon doping, the 1D Kondo insulator becomes a ferromagnet.More recently, the 1D KL has been studied using Monte Carlo [15][16][17], density matrix renormalization group (DMRG) [18][19][20][21][22], bosonization [23,24], strong-coupling expansion [25] and exact diagonalization [26].Additionally, renormalization and Monte-Carlo methods have also been used to examine the p-wave version of the 1D Kondo lattice, which exhibits topological end-states [27][28][29].
The Kondo insulator can be driven metallic by doping, which leads to a closing of charge and spin gaps, forming a Luttinger liquid with parameters that evolve with doping and J/t c [21,26].Both the insulating phase at half-filling and the doped metallic regime are nontrivial, as the k F , extracted from spin and charge correlation functions corresponds to a large Fermi surface, which counts both the electrons and spins v FS /π = n e +1.The weak-coupling J/t c 1 regime at finite doping continues to be paramagnetic, but the strongly-coupled J/t c 1 regime becomes a metallic ferromagnetic state for infinitesimal doping.In this regime, the spin-velocity goes to zero, characteristic of a ferromagnetic state [30], as inferred from spin susceptibility.
The excitation spectrum of a one dimensional Kondo lattice at half-filling was first studied by Trebst et al. [31] who employed a strong coupling expansion in J/t c to examine the one and two-particle spectrum.Their studies found that beyond t c /J > 0.4, the minimum in the quasi-hole spectrum shifts from k = π to k < π.Furthermore, they extracted the quasi-particle weights showing that Z → 0 right at t c /J = 0.4 when the dispersion is flat around k ∼ π.Smerat et al. [32] used DMRG to compute the quasi-particle energy and lifetime to verify these results and extend them to partial filling.They pointed out that the exchange of spin between conduction electrons and localized moments leads to formation of "spinpolarons", here referred to as "composite fermions".

B. Motivation and summary of results
The appearance of an insulator in a half-filled band goes beyond conventional band-theory and requires a new conceptual framework.A large body of work, dating back to the 1960s recognized that there are two ways to add an electron into a system containing localized moments [33][34][35][36], either by direct "tunneling" an electron into the system, formally by acting on the state with the conduction electron creation operator c † σ , or by "cotunneling" via the simultaneous addition of an electron and a flip of the local moment at the same site Both processes change the charge by e and the spin by one half.The object created by F † has also alternately referred to as a "composite fermion" or a "spin-polaron" [32].Here we will employ the former terminology, introducing as the composite fermion creation operator: F † β transforms as a charge e and spin S = 1/2 fermion, and with the above normalization the expectation value of its commutator with the conduction electron operator vanishes {c α , F † β } = 0, while that of commutator with itself is unity {F α , F † β } = δ αβ , in the strong Kondo coupling limit.
Co-tunneling lies at the heart of the Kondo problem, and insight into its physics can be obtained by observing that in the interaction, the object that couples to electron in the Kondo interaction is a composite fermion, for In certain limits, such as the large J limit and large-N limit, F behaves as a physically independent operator, suggesting that the Kondo effect involves a hybridization of the conduction electrons with an emergent, fermionic field.The large-N limit accounts for the emergence of the independent composite fermions as a consequence of a fractionalization of the local moments, and in this limit, both the composite fermion and the local moments are described in terms of a single f -electron field, Though the Kondo Lattice is a descendent of the Anderson lattice, it exists in its own right.In particular, rather than the four-dimensional Hilbert space of an electron at each site, the spins have a two-dimensional Hilbert space.If there are "f" electrons they are by definition, Z = 0 quasiparticles, as there is absolutely no localized electron Hilbert space.Field theory and DMRG of single impurity provide a clue: the presence of many body poles in the conduction self-energy can be interpreted in a dual picture as the hybridization of the conduction electrons with fractionalized spins.
One of the key objectives of this work is to shed light on the quantum mechanical interplay between the composite fermion, the conduction electron and the possible fractionalization of local moments in a spin-1/2 1D Kondo lattice (1DKL).This is achieved by carrying out a new set of calculations of the dynamical properties of the Kondo lattice while also comparing the results with those of large N mean-field theory and strong coupling expansions about the large J limit.In each of these methods, we evaluate the joint matrix Green's function describing the time evolution of the conduction and composite fermion fields following a tunneling or cotunneling event.
Matrix-product states are ideally suited to one dimensional quantum problems, permitting an economic description of the many-body ground-state with sufficient precision to explore the correlation functions in the frequency and time domain.Here, we take advantage of this method to compute Green's function matrix between conduction electrons and composite fermions and to compare the spin correlation functions of the local moments and composite fermions.For simplicity, we limit ourselves to zero temperature T = 0.At the one-particle level we find that by analyzing the Green's function matrix between c and F , we are able to show that these operators define a hybridized two-band model, in agreement with the large-N limit.The evolution of our computed one-particle spectrum with t c /J is consistent with earlier strong-coupling expansions.Rermarkably, the shift in the minimum of the quasiparticle dispersion seen in the strong-coupling expansion at t c /J = 0.4 [31] can be qualitatively accounted for in terms of the evolution of the hybridization between conduction electrons and composite fermions.
Moreover, by calculating the dynamical spin susceptibility using MPS methods, and comparing the results with mean-field theory, we are able to identify a continuum in the spin excitation spectrum that is consistent with the fractionalization of the local moments into pairs of S = 1/2 excitations.Our strong coupling expansion coincides with the matrix product state calculation in the large J limit and we also see signs of the formation of S = 1 paramagnon bound-states below the continuum: a sign of quiescent magnetic fluctuations.

II. MODEL AND METHODS
The model we consider is deceptively simple.It is given by the one-dimensional Kondo lattice Hamiltonian where c † i,σ creates an electron of spin σ at site i and t controls the electron tunneling between sites.The operator S i is an immobile S = 1/2 spin located at site i and FIG. 1. Diagrammatic representation of calculating the Green function by MPS methods.The MPS |0 (blue) is the ground state found using DMRG.The small circles in (a) represent the single-site operators |F † x 2 (magenta circle) and cx 1 (green circle).They can be placed at any sites x1 and x2 (though requiring separate computations), giving access to the Green function in real space.The time evolution operator (orange rectangle) is split into two halves, each half approximated by a sequence of unitary gates (dark green rectangles) using a Trotter approximation.The Green functions is found by computing the overlap of the two independent time-evolved wavefunctions.In the bottom right region, we demonstrate the procedures during every step in the time evolvolution (red shaded region).The gate tensors are contracted with MPS tensors, followed by a singular value decomposition (SVD) to reorganize the tensors back into MPS form but with a increased bond dimension.
(c † i σ c i )• S i is a Heisenberg coupling between the spin moment of an electron at site i and the spin S i .In the limit of large J/t c the half-filled ground-state is composed of a product of Kondo singlets at every site, a state that is self-evidently an insulator.The challenge then, is to understand how this state evolves at finite J/t c .

A. Matrix Product State Methods
The primary tool we will use to study the properties of the Kondo lattice model will be matrix product state (MPS) tensor networks.An MPS is a highly compressed representation of a large quantum state as a contraction of many smaller tensors and is the seminal example of a tensor network.In contrast to other numerical or analytical approaches, MPS methods work well for both weakly and strongly correlated electronic systems and do not suffer from a sign problem as in the case of quantum Monte Carlo methods.The main limitation of MPS is that they are only efficient for studying one-dimensional or quasione-dimensional systems.
The two key MPS techniques we use in this work are the density matrix renormalization group (DMRG) al- gorithm for computing ground states in MPS form [37], and the time-evolving block decimation (TEBD) or Trotter gate method for evolving an MPS wavefunction forward in time [38][39][40].Our implementation is based on the ITensor software [41].
Our computational approach is illustrated at a high level in Fig. 1, using the example of computing After computing an MPS representation of the ground state |0 using DMRG, we act with F † x2 on one copy of |0 and with c † x1 on another copy of |0 .The first copy is evolved forward in time by acting e −iHt/2 using a Trotter decomposition of the time evolution operator, and the second is evolved similarly but acting with e iHt/2 .Finally, the Green function is computed from the overlap of the resulting MPS.We give additional technical details of our computational approach in Appendix H.

B. Strong coupling expansion
An insight into the nature of ground state and its elementary excitations can be obtained in the strong Kondo coupling limit.At t c /J K = 0 the ground state is a product state of local Kondo singlets The ground state is Here, ⇑ and ⇓ refer to the spins (magnetic moments) and ↑ and ↓ refer to conduction electrons.The spin-1/2 excitations corresponding to addition/removal of electrons and spin-1 excitations of changing local singlets into triplets.At finite t c /J K electrons hop to nearby sites, creating holon-doublon virtual pairs [Fig.2(b)].Consequently, the vacuum contains short-lived holon-doublon pairs, which lead to short-range correlations.

III. COMPOSITE FERMIONS: SINGLE-PARTICLE PROPERTIES A. More details on the composite fermion operator
To characterize the single-particle excitations of the Kondo lattice, observe that acting on the ground state by the operators c † ↑ and c † ↓ S + each create charge-1, spin-1/2 excitations.However, instead of c † ↓ S + , we will find that the composite fermion operator is the more natural operator to consider.One motivation is that F † σ transforms under the S = 1/2 representation of SU (2).A more intuitive motivation is that the spinelectron interaction term in the Kondo lattice Hamiltonian Eq. ( 5) can be written as ( σ is the operator which couples to electronic excitations.The factor of 2/3 in Eq. ( 7) has been chosen so that the commutator (see Appendix A for the proof) is unity in the strong coupling limit J/t c 1.The second line spoils the canonical anti-commutation of F operators, however, in the strong coupling limit J/t c 1 the expectation value of the second term is zero in the ground state, indicating that {F α , F † β } = δ αβ has canonical anti-commutation on average.Within the triplet sector, the expectation value of the anticommutator becomes δ αβ /9 and within the holon/doublon manifold, the expectation value of the anticommutator depends on the state of the magnetic moment.The overlap between the original c and the composite F electrons is The right-hand side has zero average (but finite fluctuations) in the strong-coupling ground state, suggesting that c and F create independent excitations in average.However, F σ and c σ overlap due to quantum fluctuations, motivating us to compute the full Green's function matrix involving both operators to study their associated excitations in a controlled way.This approximately particle-like behavior of the composite-fermion F σ has strong resemblance to the twoband model of heavy-fermions obtained in the large-N mean-field theory.In such a model the spin is represented using Abrikosov fermions S = 1 2 f † σf and the constraint f † f = 1 is applied on average using a Lagrange multiplier.Within mean-field theory, the Kondo interaction leads to a dynamic hybridization between f -electrons and c-electrons [c.f.Eq. ( 3)] The similarity of the two results suggests F σ ∼ f σ , implying the fact that the spin is fractionalized into spinons.
In fact we can define In the rest of this section we will confirm the picture outlined above by computing the full Green function using two approaches.We first use time-dependent matrix product state techniques on finite systems, then carry out a strong-coupling analysis to shed further light on the results.

B. Matrix Spectral Function
To examine the independence of the c and F fields, it is useful to combine them into a spinor allowing us to define a retarded matrix Green's function where θ(t) is the step function.G defines a matrix of amplitudes for the c and F fields.The G cF component determines the amplitude for a composite F to convert to a conduction electron.We are primarily interested in the properties of a translationally invariant Kondo lattice, with momentumspace Green's function In our numerical calculations, we estimate this Green's function using the expression for a translationally invariant system simply applied to finite size Green's function G(x i , x j ; t).We then perform a discrete Fourier transform on G to obtain where ∆t = T /N t is the spacing of the N t time-slices over the total duration T of the time evolution, t j = j∆t and the frequencies are sampled at the values J=0.9Although we independently compute the four components of G(k, ω), the kinematics of the Kondo lattice imply that they are not independent, which provides us a means to test and interprete our calculations.From the Heisenberg equations of motion of the conduction electron operators in the translationally invariant limit, where c (k) = −2t cos(k) is the dispersion of the conduction electrons and F kα = L −1/2 x e −ikx F xα is the Fourier transform of the composite fermion.It follows that When we transform these equations into the frequency domain, replacing i∂ t → z = ω + iη, we see that G cc and where we have suppressed the (k, z) label on the propagators, and g c = [z − c (k)] −1 is the bare conduction electron propagator.Although these equations closely resemble the Green's functions of a hybridized Anderson model, with hybridization 3J/2, we note that G F F represents a composite fermion.
From these results, it follows that without any approximation, the inverse matrix Green's function can be written in the form where g F (k, z) is the one-particle irreducible composite Green's function, determined by This quantity corresponds to the unhybridized composite fermion propagator.By reinverting (20) we can express the original Green's functions in terms of g F (k, z) as follows These are exact results, which even hold for a ferromagnetic, J < 0, Kondo lattice.By calculating G and inverting it, we can thus check the accuracy of our calculation, and we can extract the irreducible F propagator g F (k, z).
From this discussion, we see that the G R matrix offers information about both the individual excitations and their hybridization.If the Kondo effect takes place, i.e if J > 0 is antiferromagnetic, then we expect the formation of an enlarged Fermi surface, driven by the formation of sharp poles in the composite fermion propagator g F .For example, in the special case where the Green's function g F develops a sharp quasiparticle pole, then we expect allowing us to identify V = Z(3J/2) as an emergent hybridization.

C. Spectral Functions: Numerical Results
The spectral function is associated with the Green's function by The set of (q, ω) values for which the spectral function has a maximum is the analogue of a band structure for an interacting system.We show the spectral functions computing using MPS for the cases of J/t c = 2 and J/t c = 0.9 in Fig. 3 and Fig. 4 respectively.Figs.
3(e) and 4(e) show the quantity where the denominator is (2, 2) entry of 2-by-2 matrix (G R ) −1 .This quantity can be interpreted as the Green's function of the unhybridized F electrons.
The most striking feature of spectral functions in Figs. 3 and 4 is the sharp and narrow bands indicative of longlived and dispersing quasiparticles.For the larger Kondo coupling of J/t c = 2, the spectra consist of two cosine dispersion curves − cos(k) ± ∆E/2 shifted to positive and negative frequencies.For the smaller Kondo coupling of J/ = 0.9, the dispersion can be thought to arise as the hybridization of a dispersing band (mostly c content) and a localized band (mostly F content).It is apparent that in both cases, the dispersion can be approximately reproduced using a two-band fermionic model.Assuming that this is so, the quantity Im[g F (k, ω − iη)] shown in panels 3(d) and 4(d) can be interpreted as the bare dispersion the putative F fermion would need to have in order to reproduce the observed spectral functions.In both cases, a non-zero dispersion is discernible which is more significant in the J/t c = 0.9 case.Since in absence of Kondo interaction, composite fermions are localized F i F † j ∝ δ ij , this bare dispersion is naturally associated with dynamically generated magnetic coupling between the spins due to RKKY interaction which gives rise to dispersing spinons.

D. Comparison with strong coupling and mean-field
It is natural to expect some of the numerical results to match those obtained in the strong Kondo coupling limit J K /t c 1. When t c = 0 the decoupled sites each have the spectrum where the quantum numbers S and n are the total spin and charge at that particular site.Creating or annihilating a particle from the ground state has the energy cost of E 1 = 3J/2.To understand how the ground state and single-particle excited states evolve for a finite t, we have carried out a perturbative analysis for the full 2 × 2 Green's function in the Appendix D and found that to lowest orders in t c /J, where V = E 1 .The eigenenergies are which confirms the picture of two hybridized bands.Here z is the complex frequency and k = −2t c cos k is the bare dispersion of the conduction electrons.Note that to this order, the dispersion of the bare F band is not captured in agreement with previous results [31].
The quasi-particle spectrum (28) has the same form as in the large-N mean-field theory, with the difference that the value of V is determined from self-consistent mean-field equation [see Appendix F].We have plotted A F F (k, ω) spectra in Fig. (5) along with predictions from strong-coupling expansion and mean-field theory.Overall, a good agreement is found albeit deviations start to appear at lower Kondo coupling of J = 0.9.
One artifact of the mean-field theory is that the hybridization V is systematically underestimated which can be traced back to the relation between F and f in Eq. ( 3) and (10).For example, at the strong coupling limit, the mean-field theory predicts V = J/2 [see Appendix F].In order to get an agreement, we had to re-scale V → 3V /2 when comparing mean-field results to numerical results on the interacting system.One can alternatively motivate this rescaling by viewing the mean-field theory as an effective model, where the V parameter in the mean-field is "renormalized" from the bare hybridization.

E. Evolution of Single Particle States
A vivid demonstration of the particle nature of F excitations can be seen by a calculating the motion of a composite fermion wavepacket.Here, it proves useful to take account of the spatially dependent normalization of the composite fermions, defining a normalized composite fermion as follows FIG. 5. A comparison of AF F (k, ω) with the dispersion from mean-field theory (solid line) and strong-coupling expansion (dotted line).Left and right panels are J/tc = 2 and J/tc = 0.9, respectively.For J/tc = 0.9, there is a noticeable deviation of the perturbative or mean-field results from the numerical results near k = 0 for the upper band and near k = π for the lower band.
where the normalization, is calculated from the measured expectation value of the anticommutator Here the second expression follows from particle-hole symmetry.This normalization guarantees that the expectation value of the anti-commutator is normalized In the ground-state, Z(x) is a constant of motion, which with our definition of F xσ (7), is unity in the strong-coupling limit.However, at intermediate coupling, Z(x) becomes spatially dependent near the edge of the chain.
Consider a wave-packet where is a normalized wave-packet centered at y with momentum k 0 .The time-evolution of this one-particle-state will give rise to a state of the form where the . . .denotes the many-particle states that lie outside the Hilbert space of one conduction and one composite fermion.Taking the overlap with the states f † xnσ |GS and c † xnσ |GS , the coefficients of the wavepacket can be directly related to the Green's functions as follows and similarly where E g is the ground-state energy and Using the Green's functions computed from the MPS time-evolution, we can thus evaluate the time-evolution of the wave-packet.Figure 6 shows the evolution of the probability density |φ f (x, t)| 2 + |φ c (x, t)| 2 of an initial Gausian wavepacket for two values of J/t = 2 and J/t = 0.9.In the former case the composite fermion wavepacket moves ballistically until it is scattered by the boundary of the system.In the J/t = 0.9 case, however, the wave-packet undergoes significant dispersion and decay with distance, appearing to "bounce" long before reaching the wall.One possible origin of this effect, is the break-down of the Kondo effect in the vicinity of the wall, due to a longer Kondo screening length ξ = v F /T K .

F. Interpretation of single-particle results
One of the most remarkable aspects of this comparison, is the qualitative agreement between the spectral functions derived from the matrix product, strong coupling and large N expansion.In all three methods, we see that the description of the spectral function requires a two-band description.Our matrix product simulation shows that the composite fermion propagator g F contains sharp poles at k = ±π/2, ω = 0, which reflect a formation of composite fermion bound-states, as if the F fields behave as sharp bound-states.
The single-particle excitation spectrum exhibits a coherent two-band fermionic model which continues to low J/t c .This suggests that the composite F -excitations, behave as bound-states of conduction electrons and spin flips of the local moments, forming an emergent Fock space that is effectively orthogonal to that of the conduction electrons, so that c and F fermions are effectively independent fields.In effect, the microscopic Hilbertspace of the spin degrees of freedom has morphed into the Fock-space of the F-electrons.
In the large N limit, the composite fermion F is synonymous with a fractionalization of the local moments into half-integer spin fermions, moving under the influence of an emergent U (1) gauge field that imposes the constraints.From the single-particle excitation spectrum alone, aside from hybridization with conduction electrons, these emergent F fermions appear to be free excitations: the comparison with mean-field theory suggests that the original spin is fractionalized to S F = F † σ 2 F .As shown in A, in the strong coupling regime S F ∼ 1 3 S. How accurate and useful is this picture?If F electrons are indeed free beyond one-particle level, their higherorder Green's functions (including two-point functions of S) would factorize into spin-1/2 fermions.We now investigate this possibility.

IV. COMPOSITE FERMIONS: TWO-PARTICLE PROPERTIES AND SPIN SUSCEPTIBILITY
Next, we turn to the two-particle spectrum and focus on the spin susceptibility which can be probed experimentally.This function satisfies the sum-rule dωχ S (q, ω) = 2π S z = 0 for any q.Fig. 7 shows χ S (q, ω) for two values of J/t, computed from the Fourier transform of χ(x 1 , x 2 , t) and using the same Fourier transform procedure as in Eq. (H9).A broad incoherent region and at least one sharp dispersing mode (at low positive frequency) is visible.The latter is more pronounced at higher J/t = 4 compared to J/t = 1.8.A spin-flip creates a localized triplet.Since only the total magnetization is conserved, the triplet can move in the lattice forming a coherent magnon band.However, in this interacting system, the magnon can decay into many-body states and the reduced weight of the coherent band is compensated by the incoherent portion of the spectrum.
In the previous section, based on the behavior of F particles we conjectured that the spin S is proportional to S F = 1 2 F † σF .The relationship S F = 1 3 S is in fact correct in the strong Kondo coupling limit (appendix A).To test its validity beyond this limit, we compare χ S (q, ω) with 9χ F (q, ω) defined in terms of composite fermions: and involves four-point functions like We see that the two are exactly equal, proving the relation S ∼ 3 S F at least within the two-particle sector.J = 0.9 χ S (q, ω) J = 0.9 χ F (q, ω) J = 2 χ F (q, ω) The spin susceptibility χS(q, ω) and (c,d) the composite fermion susceptibility 9χF (q, ω) for two values of J/tc = 1.8 and J/tc = 4 case.The two are nearly identical with minor differences at small momenta and high frequency.
However, while this relation seem to hold, fractionalization as seen in 1D Heisenberg AFM requires the fourpoint function χ F to be expressible in terms of the convolution of two single-particle propagators.To examine this possibility, we compare the spin susceptibility χ S (q, ω) with the mean-field spin susceptibility χ M F (q, ω) computed from the convolution of two f-electron propagators Fig. 8.The mean-field dynamical susceptibility contains a continuum of excitations bordered by two sharp lines that result from the indirect gap between the f-valence and f-conduction bands (lower sharp line) and the cvalence and c-conduction bands (upper sharp line) of the fractionalized Kondo insulator.A particularly marked aspect of the mean-field description in terms of fractionalized f-electrons is the continuum at q ∼ 0 which stretches from the hybridization gap (2V ) out to the half bandwidth of the conduction band.At finite q this continuum evolves into a characteristic inverted triangle-shaped continuum.At strong coupling, J/t = 2 χ S (q, ω) contains a sharp magnon peak, and the triangle-shaped continuum is absent.This is clearly different from χ M F (q, ω).However at weaker coupling J/t = 0.9, the MPS susceptibility is qualitatively similar to the mean-field theory, displaying the triangle-shaped continuum around q ∼ 0 and a broadened low energy feature that we can associate with the indirect band-gap excitations of the f-electrons.It thus appears that at strong-coupling, the f-electrons are confined into magnons, whereas at weak-coupling the spins have fractionalized into heavy fermions.

A. Strong coupling and mean-field perspective
To gain further insight into the dynamical spin susceptibility, we discuss the two-particle sector from both strong coupling and mean-field perspectives.It is useful J = 0.9 TN J = 0.9 MF J = 0.9 RPA FIG. 8. Spin susceptibility χS(q, ω).The first rwo shows tensor network results for (a) J/t = 2 and (b) J/t = 0.9.This is compared with large-N mean-field theory (MF) results in (c) J/t = 2 and (d) J/t = 0.9, and random-phase approximation (RPA) results in (e) J/t = 2 and (f) J/t = 0.9.The parameters in (e) and (f) are U /t = −2 and U /t = −0.75sin(q/2).
A generally q-dependent interaction between quasi-particles within RPA captures both magnon branch and the details of the correlation function at q ∼ 0.
to generalize the Hamiltonian of Eq. ( 5) by including a Coulomb repulsion U > 0, i.e.
which favors one electron per site.We assume U is small enough so that the ground state is smoothly connected to the original problem with U = 0.The starting point is that at the strong coupling, all sites are singlets, and therefore the relation holds.This means that S can be replaced by − 1 2 c † σc in χ S defined in Eq. ( 37), creating the following strongcoupling picture: A S + spin-flip can be considered as a creation of a local doublon-holon spin-triplet T + pair at the same site.Such a state has energies around E 2 = 2E 1 as shown in Fig. , 9.Under time-evolution, the doublon and holon can move around and recombine at site n where the T + triplet is annihilated.Such a T + triplet is described by Including the U interaction, each holon or doublon costs an energy E 1 +U/2 and E 2 → 2E 1 +U .By acting on this with the Hamiltonian H |F + = E |F + and projecting the result to within the two-particle excitations, we find that the wavefunction ψ(n 1 , n 2 ) obeys the first-quantized Schrödinger equation This is a two-particle problem, where the particles interact via the U = −J K − U term.Note that a repulsive/attractive interaction among electrons is an attractive/repulsive interaction among doublon and holon.
In the usual regime (U ≥ 0) the interaction U < 0 is attractive.While a continuum of excited states exists, the ground state is a stable magnon boundstate between doublon and holon, with a correlation length that diverges as U → 0. The continuum is essentially a fractionalized magnon into doublon and holon pairs as can be seen in the U = 0 case.It is natural to expect that due to interactions not considered, the doublon-holon pair decay into the ground state.For an attractive U ≤ −J K /2, the interaction between doublon and holon U > 0 is repulsive, rendering the bound-state highly excited and unstable.
The eigenstates |F + can be used to compute the spinsusceptibility χ S .The result is shown in Appendix E. The result at the strong coupling J/t = 2 contains a magnon band in good agreement with MPS results.This indicates that while a spin-flip has fractionalized into a doublon-holon pair, there are residual attractive forces in a Kondo insulator that bind the two.
On the other hand, at the weak coupling J/t = 0.9 limit, the strong-coupling analysis is incapable of reproducing MPS results around q ∼ 0. This suggests that U renormalizes to zero in the small momentum limit.As seen in Fig. 8(b) in the weak-coupling limit, the strong magnon-like resonance in the MPS dynamical spin susceptibility, broadens and merges with the triangular feature around q ∼ 0, in sharp contrast to the strong coupling results and more closely resembling the mean-field theory [Fig.8(d)].
To capture the doublon-holon interaction and the magnion band, the mean-field theory be improved by including a momentum-dependent residual interaction U (q) between f-electron quasi-particles within a randomphase-approximation (RPA) framework.The resulting susceptibility can be written as The RPA results is shown in Fig 8(e,f) for a constant U in the strong coupling J/t = 2 case and U (q) ∼ q in the weak-coupling J/t = 0.9 case, both in good agreement with the MPS results.
Overall, these results indicate that the low-lying spin-1 charge-neutral excitation of the ground state can be regarded as fractionalized into spin-1/2 charge-e single particle excitations that have some residual attraction in one-dimension, forming a magnon branch in the dynamical spin-susceptibility.The disappearance of this magnon branch at q → 0 in the weak-coupling regime and its comparison with RPA suggests that at long distances the residual interaction disappears, leading to deconfined quasi-particles.

V. CONCLUSION
By contrasting strong coupling, mean-field theory and matrix product calculations of the dynamics of the one dimensional Kondo insulator, we gain an important new perspective into the nature of the excitations in this model.There are a number of key insights that arise from our results.
Firstly, we have been able to show that the composite fermion, formed between the conduction electrons and localized moments behaves as an independent fermionic excitation, giving rise to a two-band spectrum of charge e, spin-1/2 excitations, with hybridization between the electrons and the independent, composite fermions.Our results are remarkably consistent with the mean-field treatment of the Kondo insulator.
By contrast, our examination of the dynamical spin susceptibility paints a more nuanced picture of the multiparticle excitations.At strong-coupling, we can explicitly see that the triplet holon and doublon combination created by a single spin-flip form a bound magnon, giving rise to a single magnon state in the measured dynamical susceptibilty.Thus at strong coupling, the spin excitation spectrum shows no sign of fractionalization.On the other hand, it can be easily checked that spin-singlet charge-2e excitations are always deconfined.Essentially, two dobulons (or two holons) can never occupy the same site, very much as same-spin electrons avoid each other due to Pauli exclusion, and thus do not interact.
However, at weaker coupling, the dynamical susceptibility calculated using MPS methods, displays a dramatic continuum of triplet excitations with an inverted triangle feature at low momentum, characteristic of the direct band-gap excitations across a hybridized band of conduction and f-electrons, and high momentum feature that resembles the indirect band-gap excitations of heavy f-electrons.These results provide clear evidence in support of a fractionalization picture of the 1D Kondo in-sulator at weak coupling.Based on these results, it is tempting to suggest that there are two limiting phases of the 1D Kondo insulator: a strong coupling phase in which the f-electrons are confined into magnons, and a deconfined weak-coupling phase where the local moments have fractionalized into gapped heavy fermions The emergence of a continuum in the spin-excitation spectrum at weak coupling may indicate that that the confining doublonholon interaction at strong coupling, either vanishes, or changes sign at weak coupling, avoiding the formation of magnons.

A. Further Directions
It would be very interesting to extend these results to two dimensions.The strong-coupling analysis of the composite fermion Green's function and the doublonholon bound-states can be extended to higher dimensions, where it may be possible to calculate a critical J at which confining doublon-holon bound-state develops.Further insight might be gained into the two-dimensional dimensional Kondo insulator using matrix-product states on Kondo-lattice strips, or alternatively, by using fully two dimensional tensor-network approaches or sign-free Monte-Carlo methods [13].
B. Discussion: Are heavy fermions in the Kondo lattice fractionalized excitations?
The 1D Kondo lattice is the simplest demonstration of Oshikawa's theorem [42]: namely the expansion of a Fermi surface through spin-entanglement with a conduction electron sea.Traditionally, the expansion of the Fermi surface in the Kondo lattice is understood by regarding the Kondo lattice as the adiabatic continuation of a non-interacting Anderson model from small, to large interaction strength [43].Yet viewed in their own right, the "f-electron" excitations of the Kondo lattice are emergent.
Our calculations make it eminently clear that in the half-filled 1D Kondo lattice, the f-electrons created by the fields form an emergent Fock space of low energy, charge e excitations that expand the Fermi sea from a metal, to an insulator.Less clear, is the way we should regard these fields from a field-theoretic perspective.From the large-N expansion it is tempting to regard heavyfermions as a fractionalization of the localized moments, S j → f † jα σ 2 β f jβ .Our calculations provide support for this picture in the weak-coupling limit of the 1D Kondo lattice, where we see a intrinsic dispersion of the underlying F electrons, reminiscent of a spin liquid, and a continuum of S = 1 excitations in the dynamical spin susceptibility.
Yet the use of the term "fractionalization" in the context of the Kondo lattice is paradoxical, because the excitations so-formed are self-evidently charged.Fieldtheoretically, the spinons transform into heavy fermions, acquiring electric charge while shedding their gauge charge via an Anderson-Higgs effect that pins the internal spinon and external electromagnetic gauge fields together.
Why then, can we not regard the f-electrons of the Kondo lattice as "Higgsed"-fractionalized excitations?This is because the classical view of confinement [44] argues that confined and Higgs phases are adiabatic limits of single common phase: i.e. the excitations of a Higgs phase are confined.
Yet on the other hand, we can clearly see the one and two-particle f-electron excitations, born from the localized moments, not only in the large N field theory, but importantly, in the matrix productstate calculations of the 1D Kondo lattice.Moreover, a recent extension of Oshikawa's theorem extends to all SU (N ) Kondo lattices [45], suggests that the large N picture involving a fractionalization of spins into heavy fermions is a valid description of the large Fermi surface in the Kondo lattice.How do we reconcile these two viewpoints?Further work, bringing computational and analytic techniques together, extending our work to higher dimensions will help to clarify these unresolved questions.
Y. K. acknowledges discussions with E. Huecker.This research was supported by the U. S. National Science Foundation division of Materials Research, grant DMR-1830707 (P. C. and also, Y. K during the initial stages of the research).

Appendix A: F commutation relations
The composite fermion operators have the expression which means The same factor of 2/3 appears in strong coupling expansion of multi-channel lattices [46].The anti-commutation relations are We use the identities where This excited state has energy E λ = E 0 + 3J, so the second-order correction to the strong-coupling ground state energy is leading to the energy E g /N = −3J/2 − t 2 /3J.The correction to the wavefunction is Note that {c nσ , F † nσ } = σS z n , where σ = ± for σ =↑, ↓.Assuming k is a good quantum number we can find in terms of single doublon and holon states defined as Using these and the spectral representation we find the Green's function where V ≡ 3J/2.This can be written as So to lowest order in t we get and we can interpret V as a hybridization between the conduction and composite f-electron.The real-space Green's function can be computed from G(z, k) via where When taking Fourier transform Appendix F: Mean-field theory Representing the spin in Eq. ( 1) with fermionis S αβ = f † α f β along with a constraint f † α f α = 1 and decoupling the resulting four-fermion interaction using a Hubbard-Stratonovitch transformation, we arrive at Due to π-periodicity of the tan 2α k , we are free to choose either the period 2α k ∈ (0, π) or 2α k ∈ (−π/2, π/2).We choose the former interval, because the angle evolves more continuously in the Brillion zone.Therefore, The relation between Kondo couling and the dynamic hybridization is given by Assuming f = 0, c = −2t c cos k, in the continuum limit, where K(k) is the complete elliptic integral of the first kind.The strong-coupling (large V ) limit of this integral is V → J/2.We can use this mean-field theory to compute the retarded Green's function as well as (anti-)time-ordered Green's functions These together with G R = θ(t)(G T − G T ).
Appendix G: Two-particle excitations -Random phase approximation In the non-interacting limit, the only contribution to Eq. (E16) is the disconnected part coming from Wick's contraction χ S (q, τ ) = −T τ c n↓ (τ )c † 0↓ T τ c † n↑ (τ )c 0↑ (G1) For non-interacting systems, we get the usual result Therefore, we could assume that this is just the non-interacting G 0 q (τ ) but multiplied by the factor e −(3J K +2U )τ .Furthermore, the hopping of holons and doublons is exactly the same.So, we propose χ dis.S (q, ω + iη) where This is shown in the figure.However, the magnon band is missing, even if we include RPA: Both term are there, but for ω > 0 only the first term contributes to the imaginary part χ f f (q, ω > 0) = π k cos 2 α k sin 2 α k+q × (I6) There were some numerical errors previously.Next, we plot this function assuming f = 0. We can also do RPA on this.We define χ RP A (q, ω) = 1 χ −1 (q, ω) − W (I7) Appendix J: Parallelization of MPS Calculations The Green's function at a given time is a matrix defined in (x 1 , x 2 ) domain.The calculation of each entry (x 1 , x 2 ) involves independent time-evolution calculations and overlaps of different wave functions |u (x, t) s and |u (x, t) s.These wavefunctions originate with creation and annihilation operators acting on different sites x).So we can parallelize these calculations and significantly reduce the time to solution.For each time slice or value of t, the computation contains two parts, the time evolution and measurement.
The time evolution of |u (x, t) s and |u (x, t) s are independent of each other and consume approximately same amount of time, which can run in different threads with minor data exchange.In total, there are O(N ) number of wave functions, which can be parallelized with no overhead cost, and scale well with increasing number of threads.
The measurement of the Green's functions matrix involves calculating the overlap of |u (x, t) s at different sites x 1 and x 2 .Both x 1 and x 2 run from 1 to N .And the computation of these overlaps are independent, which can be computed with O(N 2 ) threads.
The time evolution step takes the dominant amount of time, because each time evolution requires application of series of gates and repeating singular value decomposition to keep bond dimension increasing, which contributes to a large prefactor before O(N ).Though the measurement scales as O(N 2 ), the overlap operation is much faster.

FIG. 2 .
FIG. 2. Strong coupling diagram.(a) Ground state is comprised of local singlets between spins and conduction electrons.(b) The hopping term in the Hamiltonian creates doublon-holon pairs whose corresponding spins are together in a singlet states, i.e. charge-0, spin-singlet admixture.

2 FIG. 3 .
FIG. 3. Numerical results for the spectral function for J/tc = 2. (a) The conduction electron component of spectral function (b) The composite f fermion component of spectral function.(c) The line plots of both c(blue) and f(red) fermion spectral functions.(d) and (e) The real and imaginary part of gF (k, ω) defined in Eq. (25)

FIG. 4 .
FIG. 4. Numerical results for the spectral function for J/tc = 0.9.(a) The conduction electron component of spectral function.(b) The composite f fermion component of spectral function.(c) The line plots of both c(blue) and f(red) fermion spectra functions.(d) and (e) The real and imaginary part of gF (k, ω) defined in Eq.(25)

FIG. 9 .
FIG. 9.The spin-flip is equivalent to creation of a triplet doublon-holon pair at the same site.The pair can move together as a magnon or decay into fractionalized doublon and holon.The former has lower energy if U = JK + 2U > 0.

FIG. 10 .
FIG. 10.Graphical illustration of what is done here.(a) |ψg is the strong coupling ground state.(b) the result of acting with Ht on |ψg .(c) Acting with cjσ or Fjσ creates a holon.(d) The holon moves due to Ht moving around singlets, but also high-energy triplets.(f) The final result after projection to low-energy singlet sector.
F1) where c = −2t c cos k and f = 0 and the Lagrange multiplier λ imposes the constraint on average.At p-h symmetry, considered here, λ = 0.The Hamiltonian (F1) can be diaganalized using a SO(2) rotationc kσ f kσ = cos α k − sin α k sin α k cos α k
there will be virtual doublon-holon pairs |C n+1,n whose corresponding spins are in a singlet state |S n+1,n .