{'fraction_non_alphanumeric': 0.053130850842521324, 'fraction_numerical': 0.01945773524720893, 'mean_word_length': 4.436218335343788, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum many-body states and Green functions of nonequilibrium electron-magnon systems: Localized spin operators vs. their mapping to Holstein-Primakoff bosons

The operators of localized spins within a magnetic material commute at different sites of its lattice and anticommute on the same site, so they are neither fermionic nor bosonic operators. Thus, to construct diagrammatic many-body perturbation theory, the spin operators are usually mapped to the bosonic ones with Holstein-Primakoff (HP) transformation being the most widely used in magnonics and spintronics literature. However, to make calculations tractable, the square root of operators in the HP transformation is expanded into a Taylor series truncated to some low order. This poses a question on the range of validity of truncated HP transformation when describing nonequilibrium dynamics of localized spins interacting with each other or with conduction electron spins. Here we apply exact diagonalization techniques to Hamiltonian of fermions (i.e., electrons) interacting with HP bosons vs. Hamiltonian of fermions interacting with the original localized spin operators in order to compare their many-body states and one-particle equilibrium or nonequilibrium Green functions. The Hamiltonian of fermions interacting with HP bosons gives incorrect ground state and electronic spectral function, unless large number of terms are retained in truncated HP transformation. Furthermore, tracking nonequilibrium dynamics of localized spins over longer time intervals requires progressively larger number of terms in truncated HP transformation. Finally, we show that recently proposed [M. Vogl et al., Phys. Rev. Research 2, 043243 (2020); J. K\"{o}nig et al., SciPost Phys. 10, 007 (2021)] resummed HP transformation resolves the trouble with truncated HP transformation, while allowing us to derive an exact (manifestly Hermitian) Hamiltonian consisting of finite and fixed number of boson-boson and electron-boson interacting terms.


I. INTRODUCTION
The concept of spin waves was introduced by Bloch [1] as a disturbance in the local magnetic ordering of ferromagnetic materials. In the spin wave, the expectation value of localized spin operators precess around the easy axis with the phase of precession of adjacent expectation values varying harmonically in space over the wavelength λ. The quanta of energy of spin waves behave as quasiparticles termed magnons each of which carries energy ω and spin . As regards terminology, we note that in spintronics and magnonics [2] literature it is common to use "spin wave" for excitations described by the classical Landau-Lifshitz-Gilbert (LLG) equation [3] within numerical micromagnetics [4] or atomistic spin dynamics [5], while "magnon" is used for quantized version of the same excitation. In other subfields of condensed matter physics, terms "spin waves" and "magnons" are sometimes used to distinguish between long-and short-wavelength excitations, respectively, or both names are used interchangeably [6].
The second-quantization description of magnons was introduced by Holstein and Primakoff (HP) [7] by mapping the localized spin operatorŜ i on site i of the lattice to bosonic operatorŝ Hereâ † i (â i ) creates (annihilates) HP boson on site i and satisfies the bosonic commutation relations where 1 is the unit operator in the infinite dimensional Hilbert space of bosons. The HP boson number operator, n i =â † iâ i whose eigenvalues and eigenstates are defined arXiv:2104.07657v2 [cond-mat.str-el] 1 Jul 2021 byn i |n = n|n , measures how much the localized spin deviates away from the ground state [where the ferromagnetic ground state with the z-axis as the easy axis is assumed in Eq. (1)]. Thus, the creation of one HP boson is equivalent to removing of one unit of spin angular momentum from the ground state [see Fig. 1(c) for illustration and Sec. II I for technical details]. The textbook literature [8,9] is typically focused on band structure of noninteracting magnons (which can also be topologically nontrivial [10,11]), so it discusses only the lowest-order truncation of the original HP transformation in Eq. (1) while retaining the terms in the Hamiltonian that are up to the quadratic order in the bosonic operators. This effectively assumes low-density limit n i /2S 1 achieved at, e.g., sufficiently low temperatures [12] and/or large S 1 in which HP bosons can be treated as noninteracting. Taking into account higher order terms in the Hamiltonian generated by Eq. (3), as well as in the Taylor expansion of the square root in Eq. (1), produces higherthan-quadratic terms in the bosonic operators which describe boson-boson interactions [6,9,[11][12][13][14][15] leading to renormalization of magnon energy, magnon decay (one magnon decays into two) [6], coalescence (two magnons coalesce into one), four-magnon interactions, decays into four magnons and other higher order processes [16].
Since bosonic operatorsâ † i ,â i act on an infinitedimensional Hilbert space, but the physical Hilbert space corresponding to a single localized spin on site i is spanned by only 2S +1 states, the extra unphysical states are decoupled from the physical ones by the square root in Eq. (1). Such exact HP transformation in Eq. (1) splits the infinite dimensional Hilbert space spanned by boson number states {|n } n∈N into two sectors-physical states {|n } n=0,...,2S ; and all the unphysical ones {|n } n>2S [see also Eq. (36)]. Those sectors cannot be connected bŷ S + i andŜ − i operators. However, when the square root in Eq. (1) is expanded in power series and then truncated (see Secs. II E and II F) to any finite order N T , the physical and unphysical subspaces become coupled. In addition, canonical commutation relations for the spin operators are then satisfied only approximately, resulting in artificial breaking of rotational symmetries that may be present in the original Hamiltonian [17,18].
Retaining higher order terms in the Taylor series expansion of Eq. (1) is necessary to study, e.g., equilibrium properties at increasing temperature [16,19,20] or nonequilibrium dynamics [13][14][15]21]. For example, Dyson [19] calculated how magnetization of the Heisenberg model of a three-dimensional ferromagnet decays with temperature, M (T )/M (0) = 1 − c 0 T 3/2 − c 1 T 5/2 − c 2 T 7/2 −c 3 T 4 +O(T 9/2 ), where T 3/2 is the so-called Bloch law for noninteracting magnons with parabolic energymomentum dispersion; second and third term also stem 2 considered in this study) which interact with each other via the nearest-neighbor Heisenberg exchange interaction of strength JH . At t = 0, nonequilibrium dynamics can be initiated by flipping the localized spin on site i = 1. The 1D quantum many-body system in (b) is the same as (a) but with smaller number of N = 3 sites whose spin-S localized spins additionally interact with conduction electron spins (blue arrows) via the sd exchange interaction of strength J sd . The conduction electrons hop between the sites with the hopping parameter γ where "half-filled" (i.e., Ne = 3) tight-binding chain is used. Panel (c) illustrates the reduction of eigenvalues S z i of the z-componentŜ z i of localized spin-5 2 operator on site i = 1 by Nmag units, which is equivalent to creating Nmag Holstein-Primakoff bosons on site i = 1 once the localized spins are mapped [7] to bosonic operators. from noninteracting magnons but with nonparabolic dispersion on a discrete lattice; and magnon-magnon interactions start manifesting at order T 4 . Such calculations require many-body perturbation theory (MBPT) [22,23] which is most easily formulated in terms of bosonic or fermionic operators. For such operators, the Wick theorem [24] for their averages over the noninteracting system makes it possible to expand properties of the interacting system into the Feynman diagrammatic series of perturbation order O(g n ) where g is the strength of fermionboson interaction. On the other hand, spin operators which commute on different sites and anticommute on the same site are neither fermionic nor bosonic operators, and there is no Wick theorem for spin operators.
Since higher order terms in the power-series expansion of square root in the HP transformation in Eq. (1), conjectured by Kubo [25] to be only an asymptotic series, lead to cumbersome MBPT [26,27], a plethora of other mappings of original localized spin operators to bosonic or fermionic operators has been proposed. For example, one has a choice to map to Dyson-Maleev bosons [19,28], Schwinger bosons [21], fermions [29,30], Majorana fermions [31], sypersymmetric operators [32] and exotic particles called semions [33]. We note that mapping of localized spin operators to bosonic or fermionic operators can be evaded altogether within the path integral formulation by using spin coherent states, but that approach leads to topological terms associated with the Berry phase so that even in path integral formalism mapping to bosonic operators is preferred [21]. Although the Dyson-Maleev transformation evades usage of the square root of operators in Eq. (1), it generates Hamiltonian that is no longer manifestly Hermitian. The Schwinger transformation require the introduction of auxiliary fields. The Jordan-Wigner transformation [29] or mapping to Majorana fermions [31] are exact, but they work only for S = 1/2 operators.
These drawbacks have prompted very recent reexaminations [17,18] of HP transformation to find possible nonperturbative replacements of the Taylor series of the square root in Eq. (1) which can be written as a polynomial in bosonic operators, while ensuring no coupling between physical and unphysical subspaces as well as manifestly Hermitian bosonic Hamiltonian. Although such polynomial expressions do not reproduce exactly the canonical commutation relations for the spin operators, the extra terms generated turn out to be unimportant because they do not couple physical and unphysical subspaces of the bosonic Hilbert space, i.e., they act solely on the unphysical subspace [17,18].
The MBPT calculations of equilibrium magnon properties based on Dyson-Maleev vs. truncated HP transformation have been carefully compared in the literature over many decades [26,27]. On the other hand, much less is know about the range of validity [34,35] of truncated HP transformation when describing nonequilibrium dynamics of localized spins, including situations where additional interactions with conduction electrons are present. The electron-localized-spin interactions are frequently encountered in quantum transport phenomena in spintronics. The nonequilibrium MBPT [22,23] for such problems is virtually always conducted using truncated HP transformation, as exemplified by theoretical and computational modeling of inelastic electron tunneling spectroscopy in magnetic tunnel junctions [36]; spin-transfer [37][38][39] and spin-orbit torques [40]; ultrafast demagnetization [41]; and conversion of magnonic spin currents into electronic spin current (or vice versa) at magnetic-insulator/normal-metal interfaces [42][43][44]. Similarly, truncated HP transformation is typically chosen for problems in quantum magnonics, such as for nonequilibrium dynamics of localized spins within magnetic insulators [13][14][15]45]; their interaction with external electromagnetic fields [12,46]; and analysis of coherence of magnon quantum states [14,47]. One can expect that truncated HP transformation will eventually break down at sufficiently long times (as confirmed Figs. 4 and 5) when higher-order terms in the expansion of the square root in HP transformation become important. Such breakdown then precludes [37] accurate tracking of nonequilibrium dynamics of localized spins, which can be driven far from their initial direction (along the easy axis) and eventually reversed by, e.g., spin-transfer torque [48,49]. Surprisingly, rigorous analysis of such breakdown is lacking.
Instead, current-driven magnetization reversal via spin-transfer torque [48] is standardly modeled by the LLG equation [50], which is combined in a multiscale fashion with some type of steady-state [51] or timedependent quantum transport calculations [52][53][54][55][56][57][58] considering single-particle quantum Hamiltonians for electrons. Thus, such hybrid quantum-classical theories are justified only in the classical limit → 0 and for large localized spins S → ∞ (while S × → 1) [3,57,59], as well as in the absence of entanglement [3,49,60] between quantum states of localized spins. For example, in the emerging concept of quantum spin torque [49,[60][61][62][63], describing transfer of angular momentum between spins of flowing electrons and localized spins in situations [64] where the latter must be described by quantummechanical operators, the whole system of electrons and localized spins can only be modeled by a quantum manybody Hamiltonian [as exemplified by Eqs. (4) and (72)].
In this study, we apply exact diagonalization techniques [65] to quantum many-body Hamiltonians defined on a one-dimensional (1D) chain of N sites hosting fermionic (for electrons) and localized spin operators, or fermionic and bosonic (obtained by mapping the original localized spin operators) operators. By comparing their many-body quantum states and Green functions (GFs), both in equilibrium and in nonequilibrium, makes it possible to precisely delineate the range of validity of truncated HP transformation. We consider 1D quantum Heisenberg ferromagnetic spin-S XXX chain hosting localized spins which interact via the nearestneighbor exchange interaction of strength J H , as illustrated in Fig. 1(a), where both spin S = 1 (as the "ultraquantum" limit) and S = 5/2 (as in, e.g., Fe 3+ valence state with five 3d electrons coupled by Hund's rule into the high spin state forming a localized S = 5/2 moment) are employed. Naively, the eigenvalue ofŜ 2 i being S 2 (1 + 1/S) suggests that quantum effects become progressively less important for S > 1, but they exist for all S < ∞ vanishing as 1/2S in the classical limit [66]. The nonzero electron hopping γ between the sites, where N = 3 sites are chosen when electrons are present as illustrated in Fig. 1(b), means that such 1D chain models a ferromagnetic metal (FM). Its N e conduction electrons [we consider half filled lattice, so N e = 3 for systems in Fig. 1(b)] interact with localized spins via sd exchange interaction [67] usually considered in spintronics. From the viewpoint of the physics of strongly correlated electrons, the model illustrated in Fig. 1(b) can also be interpreted as the Kondo-Heisenberg chain [68].
For technical reasons (i.e., exponential increase of the size of matrix representation of Hamiltonian), we consider 1D chains of N ≤ 7 sites while concentrating on generic features which are not bound to one dimension or small number of electrons and localized spins considered. In fact, artificial atomic chains have also been realized experimentally using ferromagnetically [69] or antiferromagnetically [70] coupled few Fe atoms on a substrate, where magnons along the chain were excited and detected via atom-resolved inelastic tunneling spectroscopy in a scanning tunneling microscope [69]. It is also worth recalling that small clusters (composed of, e.g., 2-8 lattice sites) in 1, 2, and 3 spatial dimensions-hosting electrons interacting with each other via the on-site or nearest-neighbor Coulomb interaction [71][72][73] (as described by "pure" and extended Hubbard models [71], respectively); or electrons interacting with bosons [74][75][76]-have played an important role in testing approximation schemes for quantum many-body problem against numerically exact benchmarks in different subfields of condensed matter and atomic-molecularoptical physics. Furthermore, the advent of numerically exact algorithms and supercomputers has led to recent re-examination of many physically motivated simplifications and approximations developed earlier in quantum many-body theory for condensed matter systems (such as Migdal-Eliashberg theory for electron-phonon systems [77]; partial summation of classes in Feynman diagrams in MBPT [78]; and existence of Luttinger-Ward functional of dressed one-particle Green function [79]) in order to draw boundaries of parameters for which their complete breakdown ensues. Our study proceeds in the same spirit, where we explicitly delineate "breakdown" times-in Fig. 4 for pure localized spins and in Fig. 5 for localized spins interacting with conduction electrons-at which widely used in spintronics and magnonics truncated versions of the HP transformation in Eq. (1) inevitably break down by generating quantum time evolution which starts to substantially deviate from the exact one obtained by using the original localized spin operators.
The paper is organized as follows. In Sec. II we introduce different versions of quantum many-body Hamiltonian describing systems in Fig. 1 and their matrix representations, as well as a procedure to obtain the exact oneparticle double-time-dependent retarded and lesser GFs. In particular, subsection II E introduces an infinite power series expansion of the HP transformation from Eq. (1) and its truncation, while subsection II F provides a brief summary of recently proposed [17,18] resummation of truncated HP transformation. The time evolution of quantum many-body states of a spin chain with no electrons (N e = 0) is employed in Sec. III A to examine the range of validity of truncated HP transformation when tracking time evolution of localized spins in the presence of magnon-magnon interaction and different number of initially excited magnons N mag . Then in Sec. III B we introduce electrons into 1D chain to examine the range of validity of truncated HP transformation when track-ing time evolution of localized spins in the presence of both magnon-magnon and electron-magnon interactions. In the same Sec. III B, we additionally employ resummation [17,18] of truncated HP transformation to derive quantum many-body Hamiltonian [Eq. (72)] for electronmagnon systems in terms of fermionic and bosonic operators whose usage reproduces numerically exact result from calculations based on the original localized spin operators. In Secs. III C and III D we compare ground state and electronic spectral function (or "interacting density of states" [80,81]) of quantum many-body Hamiltonian in terms of the original localized spin operators vs. Hamiltonian using bosonic operators generated by truncated HP transformation. The magnonic spectral function and related excited eigenstates are studied in Sec. III E. Since both ground and excited states of electron-magnon interacting system are many-body entangled [82], we compute their entanglement entropy in Sec. III F which makes it possible to quantify how far they are from the eigenstates of a system where the interaction between electrons and localized spins is turned off. Finally, Sec. III G studies time evolution of diagonal and off-diagonal elements of time-dependent lesser electronic and magnonic GFs which demonstrates that often employed "local self-energy approximation" [83][84][85][86] for electron-boson interacting systems, neglecting the offdiagonal elements, is generally not justified. We conclude in Sec. IV.

A. Quantum many-body Hamiltonian of electrons interacting with localized spins
The quantum many-body Hamiltonian of 1D chain composed of N sites (with open boundary conditions assumed), each of which hosts spin-S localized spin which interacts with conduction electron spins [as illustrated in Fig. 1(b)], is given by [87] H =Ĥ e +Ĥ lspins +Ĥ e−lspins .
It acts in the total space F e ⊗ H lspins which is a tensor product of the Fock space of electrons, F e , and the Hilbert space of all localized spins The Fock space of electrons [23] F e = C ⊕ H e ⊕Â(H e ⊗ H e ) ⊕Â(H e ⊗ H e ⊗ H e ) ⊕ · · ·, (6) is induced by the one-electron Hilbert space H e as the completion (indicated by overline) of the direct sum of antisymmetrized n-fold tensor products of H e . The op-eratorÂ antisymmetrizes tensors for fermionic particles. In the sector of F e ⊗ H lspins with N e = 0 electrons, we have a chain hosting only spin-S localized spins [as illustrated in Fig. 1(a)], which is described solely byĤ lspins termĤ chosen as the quantum Heisenberg Hamiltonian with the nearest-neighbor (NN) exchange interaction (as signified by ij notation) of strength J H = 1 eV. When electrons are present, they are described byĤ e term chosen as the tight-binding Hamiltonian with the NN hopping γ = 1 eV between single s-orbitals residing on each site. The Hamiltonian describing sd exchange interaction of strength J sd = 0.2 eV [67] between conduction electron spin and localized spins is given bŷ The row vector operatorψ † i = (ψ † i↑ , ψ † i↓ ) consists of operatorsψ † iσ which create an electron of spin σ =↑, ↓ on site i;ψ i is a column vector operator that contains the corresponding annihilation operators; andσ = (σ x ,σ y ,σ z ) is the vector of the 2 × 2 Pauli spin matrices as matrix representation of spin-1 2 operator of electronic spin. Using notation {Ô 1 ,Ô 2 } for the anticommutator and [Ô 1 ,Ô 2 ] for the commutator of two operatorsÔ 1 andÔ 2 , fermionic operators of electrons satisfy where 1 is 4 N × 4 N unit matrix in the antisymmetrized 2N -particle subspace of the Fock space F e . The localized spin operatorsŜ α i (α = x, y, z) on site i satisfy the angular momentum algebra The square of the localized spin operator, For computational convenience in calculations of electronic GFs, we change [88] the basis of one-particle electronic states from site basis to eigenenergy basis to obtain where c † i = (c † i↑ , c † i↓ ) is a row vector consisting of c † iσ operators which create an electron with spin σ =↑, ↓ in one-particle electronic eigenstate | i with the discrete eigenenergy i , so thatĤ e | i = i | i . These eigenenergies and eigenstates are evaluated by diagonalizing the one-particle tight-binding Hamiltonian where |i denotes s-orbital [whose coordinate representation is r|i = φ(r − R i )] of an electron centered on site i. Using change of basis transformation rules for operators in second-quantization formalism and substituting this into Eq. (9) we get Since eachĉ † i orĉ i operator is represented by 4 × 4 matrix (see Sec. II C), and eachŜ i operator is represented by a (2S + 1) × (2S + 1) matrix, the quantum many-body Hamiltonian in Eq. (4) for the chain of N sites in Fig. 1 Although systems containing larger than our choice N = 3 (when electrons are present) or N ≤ 7 (when electrons are absent) sites could be diagonalized with state-of-the-art numerical algorithms [65], we restrict our analysis to such smaller number of sites in order to make the analysis transparent and pedagogical-see, e.g., easy-to-follow visualization of ground and excited quantum many-body states depicting population of a small number of energy levels i in Figs. 7 and 9, respectively.

B. Symmetries of quantum many-body Hamiltonian
The exact diagonalization of quantum many-body Hamiltonian in Eq. (4) yields its many-body eigenenergies E k and many-body eigenstates |Ψ k . The total electron number operator is given byN The operator of total spin in the z-direction is the sum of electronic spin operators (first term) and localized spin operators (second term) along the z-axis at each site i. The many-body Hamiltonian in Eq. (4) has two symmetries encoded by the commutation relations: which is due to conservation of the number of electrons N e ; and which is due to conservation of total z-spin (electronic + localized spin) S z tot . Therefore, N e and S z tot , as eigenvalues ofN e andŜ z tot , respectively, serve as "good quantum numbers" for labeling quantum many-body eigenstates together with many-body eigenenergy E k . The effect of two symmetries in Eq. (20) and Eq. (21) can also be visualized in the matrix representation (see Secs. II C and II D) of quantum many-body Hamiltonian H [Eq. (4)]. For example, in Fig. 2(a) the matrix elements ofĤ in the original site basis are visually scattered throughout the whole matrix. However, whenĤ is represented in the basis of eigenstates ofN e , Fig. 2(b) shows that its matrix becomes block-diagonal where each block contains the nonzero matrix elements associated with states with fixed number of electrons N e . Finally, in Fig. 2(c)Ĥ is represented in the basis composed of eigenstates ofN e andŜ z tot simultaneously, which isolates additional submatrices with fixed S z tot within blocks associated to fixed N e .

C. Matrix representation of electronic creation and annihilation operators
A fermionic operator creating or annihilating electrons on a single site operate within the natural basis of kets |0 , |↑ , |↓ and |↑↓ which denote the empty state; state with one spin-↑ electron; state with one spin-↓ electron; and the state with one spin-↑ and one spin-↓ electron. Thus, these basis states are represented by column vectors In the same basis, creation and annihilation operators that act in the 1-site or 2-particle subspace which satisfy the fermionic commutation relations in Eq. (10).
If we consider two sites, then electronic creation (annihilation) operators,ψ † 1 (ψ 1 ) and ψ † 2 (ψ 2 ), act in the 2-site or 4-particle subspace H N =2 =Â(H e ⊗ H e ⊗ H e ⊗ H e ) of the Fock space F e and are represented by matrices of size 4 2 × 4 2 . For example, the action ofψ † 1σ in H N =2 is given by where 1 is 4 × 4 unit matrix. However, the action of requires [88] the permutation matrix P = diag(1, −1, −1, 1), instead of naïvely using only the unit matrix 1, in order to preserve the correct anticommutation relations of fermionic operators at different sites in Eq. (10). The next step is to construct the matrix representation of electronic creation and annihilation operators for three sites, which is done in a similar fashion [88] to furnish where each operator on the left hand side (LHS) is a 4 3 × 4 3 matrix. Equations (27)

E. Truncated Holstein-Primakoff transformation
The HP transformation shown in Eq. (1) expresses localized spin operators in terms of bosonic operators. However, to make MBPT for such bosonic operators tractable [26,27], one typically expands the square root of Eq. (1) in a power series in x =n i /2S which is further truncated [9,[11][12][13][14][15] to a finite number of terms N T . Inserting this result in Eq. (1), and using thus obtainedŜ α i in Eq. (7), we can re-writê as the sum of two terms. Herê (32) is one-particle Hamiltonian of noninteracting HP bosons covered in textbook literature [8,9], whereaŝ is composed of many-particle interacting terms that we write explicitly for truncation number N T = 1 to emphasize how nontrivial multi-boson interactions arise even in this lowest order truncated HP transformation. The bosonic operatorâ † is represented by an infinite matrix for a single site, so that matrix representation ofâ † i in the case of N sites is given bŷ where 1 is the unit matrix of the same size asâ † . The matrix representation of operatorâ i is the Hermitian conjugate ofâ † i .

F. Resummed Holstein-Primakoff transformation
In numerical calculations,â † orâ are first truncated to a finite N B ×N B matrices, so that the matrix representation of localized spin operators is then composed of matrix blocks associated with physical states {|n } n=0,··· ,2S and unphysical states {|n } n=2S+1,··· ,N B where ∆ c is the coupling between physical and unphysical states. The numerically exact computation of the square root of an operator in Eq. (1) ensures ∆ c = 0, but Taylor expansion of square root in Eq. (30) leads to ∆ c = 0 which, therefore, couples the physical and unphysical states. This feature reveals the trouble with the truncated HP transformation. Alternatively, Refs. [17,18] have recently proposed a resummed HP transformation that furnishes a polynomial expansion for the square root in Eq. (1) where the iterative relation for coefficients C n was derived in Ref. [17] by using flow-equations, whereas an equivalent closed-form expression was derived in Ref. [18] by using Newton-series expansion. Equation (37) ensures that for N max = 2S the matrix-block S ±,z phys associated with the physical states is exact, whereas coupling between the physical and unphysical states is ∆ c = 0, which makes nonzero subma-trixŜ ±,z unphys irrelevant for all practical purposes.

G. Relationship between localized spin operators and their mapping to Holstein-Primakoff bosons
For physically transparent understanding of the relationship between localized spin operators and their mapping to HP bosons, let us consider an example of 1D chain of N = 7 sites hosting spin-5 2 localized spins. We use arrows of different length to denote eigenvalues S z i of localized spin operatorŜ z i [Eq. (28c)] with m = −5/2, −3/2, −1/2, 1/2, 3/2, 5/2, respectively, as illustrated in Fig. 1(c). The ferromagnetic ground state of this system is identical to HP bosonic vacuum state |0 with zero HP bosons on each site n i = 0 and, therefore, total number of HP bosons also being zero, N mag = 0. Inside the ket vector on the left hand side (LHS) of Eq. (41), we indicate eigenstate with eigenvalue m = 5/2 of the localized spin operator S z i for all sites i = 1 to i = 4. Equation (41) is proved by noting that |0 on the right hand side (RHS) of Eq. (41) and ket on the LHS of Eq. (41) are both eigenstates of the same operatorŜ z i with eigenvalue m = 5/2 i.e., so they must be identical. Thus, creating n i = 1 or n i = 2 HP bosons on site i = 1, which we depict bŷ respectively, corresponds to reducing the size of localized spin on site i = 1 by 1 or 2 units , i.e., m = 5/2 → 3/2 in Eq. (44) and m = 5/2 → 1/2 in Eq. (45). Similarly, the state with a total of N mag = 2 HP bosons created on different sites i = 1 and i = 3 is depicted bŷ Thus, creating a total of N mag ∈ 0, 1, 2, . . . [Eq. (42)] HP bosons is interpreted physically as the reduction of the total localized z-spin by N mag units. Since in quantum state (â † i ) n |0 the expectation value of the z-component of localized spin operator is Ŝ z i = S − n, the constraint 0 ≤ n ≤ 2S (i.e., at a given site i one cannot create more than 2S HP bosons) must be obeyed in order to remain in the subspace of physical states [Eq. (36)].

H. Numerically exact time evolution of quantum many-body states
The solution of time-dependent Schrödinger equation for quantum many-body state |Ψ(t) is formally given by where T is the time-ordering operator. While many numerical algorithms are available to propagate Eq. (48), including direct computation of matrix exponential when H is time independent, in general by using sufficiently small δt and by consideringĤ(t) to be constant over such δt the Crank-Nicolson algorithm we employ offers propagation scheme that is unitary, accurate to second order in δt, and unconditionally stable [89].
Using thus obtained |Ψ(t) , the time evolution of the expectation value of the α-component of localized spin operator on site i is given by When localized spin operators are represented directly by finite size matrices in Eqs. (28) and (29), the corresponding expectation values S α i are numerically exact and, therefore, serve as a benchmark for alternative computation of the same expectation value whenŜ α i are represented by polynomial expressions in bosonic operators introduced in Secs. II E and II F.

I. From Holstein-Primakoff bosons to one-or two-magnon Fock states
In contrast to HP bosons created on a given site,â † i |0 , which are not the eigenstates ofĤ 0 in Eq. (32), onemagnon states are linear combinations ofâ † i |0 which diagonalize HamiltonianĤ 0 (but with periodic boundary conditions included) Thus, they can be visualized bosonic quasiparticle which carries momentum q (assuming 1D chains we use as examples) and angular momentum and is completely "delocalized" over all sites. Here E 0 = −2J H S 2 N is the ground state energy of a ferromagnetic spin chain. To find explicit expression for excited eigenstate |q , we consider 1D chain [ Fig. 1(a)] composed of N sites each of which is hosting spin-1 localized spin and with periodic boundary conditions so that its first and last site are coupled by J H inĤ 0 in Eq. (32). For the clarity of notation, we use ↓, , and ↑, to denote eigenstates of localized spin operatorŜ z i with eigenvalues [Eq. (28c)] m = −1, m = 0, and m = 1, respectively. The onemagnon state is then given by where x n = na 0 is the real-space position of the localized spin on site n + 1 and a 0 is the lattice spacing. The which indicates that creation of magnon with wavevector q removes one unit of total z-spin from the ferromagnetic ground state. Because of this feature, presence of one HP boson or one HP magnon is labeled by the same N mag = 1 throughout the paper. In addition, the expectation value of the localized z-spin operator at arbitrary site i shows that excitation of one HP magnon reduces the zcomponent of each localized spin by 1/N . This rigorous quantum-mechanical result justifies the LLG picture [4] of spin wave in which classical vectors of localized spins precess with frequency ω and with some small cone angle around the z-axis, while the phase of the precession of adjacent vectors varies harmonically in space over the wavelength λ.
In second-quantization description produced by HP transformation, |q =b † q |0 is one-magnon Fock state [90] where the creation operator of HP magnon is given bŷ Note that such one-magnon Fock state has been realized experimentally only very recently in a millimeter-sized ferrimagnetic crystal and detected by superconducting qubit as quantum sensor [90], thereby representing a counterpart in quantum magnonics of a single-photon detection from quantum optics. Note that in spintronics and magnonics literature [45] one also findsâ † i |0 denoted as "one magnon created in real space at position i" whileb † q |0 is "one magnon created in the reciprocal space with momentum q". However, the former is not an eigenstate of Hamiltonian in Eq. (51), while the later is, so we differentiate between them by using "HP boson" for the former and "HP magnon" for the latter. As already highlighted, for both situations we use label N mag = 1 for simplicity of notation because in both cases one unit of total z-spin is removed from the ferromagnetic ground state [Eq. (53)].
Nevertheless, we illustrate the distinction between HP boson and HP magnon by initializing N = 7 site chain [ Fig. 1(a)] in quantum state |Ψ 1 (t = 0) = | ↑↑↑↑↑↑ in Figs. 3(a) and 3(b); or in quantum state |Ψ 2 (t = 0) = |↓↑↑↑↑↑↑ in Figs. 3(c) and 3(d). This means that N mag = 1 HP boson is created on site i = 1 at t = 0 in the former case; while "full spin flip" of localized spin on site i = 1 in the latter case means that N mag = 2 HP bosons are created on site i = 1. Besides pedagogical value, such initial states and one or two magnon propagation including magnon bound states, has also been probed experimentally in ultracold atoms in an optical lattice where tracking of the localized spin expectation values is possible with single-spin and single-site resolution [91].
Since |Ψ 1 (t = 0) is not an eigenstate, it evolves in time to produce spatio-temporal profile of the expectation value Ŝ z i (t) [ Fig. 3(a)] in quantum state |Ψ 1 (t) . For quantum time evolution we use the scheme explained in Sec. II H where interacting HamiltonianĤ lspins from Eq. (5) is plugged in, but since only one HP boson is excited this is equivalent to using noninteractingĤ 0 in Eq. (32). The "white trace" in Fig. 3(a) visualizes how HP boson moves from the left to the right edge of the chain while undergoing reflection on site i = 7 at t ≈ 2.5 fs, as indicated by horizontal dashed line, followed by multiple back-and-forth reflections. Note that since 1D chain in Fig. 1(a) has open boundary condi-tions, its low-energy excited eigenstates differ [93] from textbook [8,9] HP magnons |q in Eq. (52) as eigenstates of interacting localized spin systems with translational invariance (which is, therefore, either infinite or finite length but with periodic boundary conditions). Figure 3(b) visualizes the overlap, n q (t) = | q|Ψ 1 (t) | 2 , between many-body quantum state |Ψ 1 (t) with one HP boson and one-magnon Fock state |q . Large values of n q (t) are observed in the region where q > 0 and t 2.5 fs, coinciding with left-to-right motion of HP boson in Fig. 3(a), which signifies excitation of N mag = 1 magnon with positive momentum. On the other hand, after reflection of the HP boson at the boundary (i.e., site i = 7) and t ≈ 2.5 fs, a rapid rise of n q (t) in the q < 0 region is observed which indicates excitation of N mag = 1 magnon with negative momentum. This is consistent with the intuitive picture of HP boson reflecting back-and-forth between the hard walls of our 1D chain with open boundary conditions. The Fock states of N mag = 2 magnons carrying momentum q and q are defined by [92] |q Here N is the normalization constant, and f mn (q, q ) is symmetric under exchange q ↔ q ensuring |q, q = |q , q in order to satisfy the symmetrization postulate of quantum mechanics for bosonic particles-as manifestly encoded by second-quantization formalism, Figure 3(c) plots spatio-temporal profile of Ŝ z i (t) in quantum state |Ψ 2 (t) starting from |Ψ 2 (t = 0) = |↓↑↑↑↑↑↑ .
For quantum time evolution we use the scheme explained in Sec. II H where interacting HamiltonianĤ lspins from Eq. (5) is plugged in, so that two HP bosons are correlated by: (i) bosonic statistics; (ii) interactions inĤ int [Eq. (33)] where N T → ∞. The two HP bosons propagate immediately from the left to the right for t > 0, as shown by "white traces" in Fig. 3(c). The corresponding overlap, n q,q (t) = | q, q |Ψ 2 (t) | 2 , in Fig. 3(d) is n q,q (t) = 0 for t 1 fs which is explained by Eq. (56) where two-magnon Fock state is composed of terms containing two HP bosons on different sites m = n. Since at t = 0 the two HP bosons are on the same site i = 1, we find n q,q (t = 0) = 0. However, this holds until t 1 fs (indicated by horizontal dashed line), after which the two HP bosons are physically separated in real space, as confirmed by the emergence of nonzero values of n q,q (t) thereafter. We also note that for 1 t 4 fs (indicated by horizontal dashed line) the region near q = 0.5 π/a 0 shows large values of n q,q (t), and since q = 0.5 π/a 0 is fixed for all values of q in Fig. 3(d), we can conclude that two HP bosons posses nearly the same velocity. Beyond t ≈ 4 fs, nonzero values of n q,q (t) in the region with q < 0 and q > 0 coexist, which indicates that one HP boson moves toward the right while the other moves toward the left edge of the chain.

J. Retarded and lesser one-particle Green functions
The fundamental quantities of nonequilibrium GF formalism [22,23] for fermions are the one-particle retarded GF and the one-particle lesser GF which describe the density of available quantum states and how electrons occupy those states, respectively. Here Θ(t−t ) is the Heaviside-function;ĉ iσ (t) indicates Heisenberg picture time evolution ofĉ iσ ; and · · · = Tr(ρ · · · ) is the quantum statistical average, whereρ is the density operator of the system at t = 0. Analogously, the bosonic one-particle retarded GF is defined by and the lesser GF is defined by In equilibrium or in steady-state nonequilibrium, these GFs depend solely on τ = t − t and can be Fourier transformed to energy domain [36], such as for electrons; and for bosons.

K. Spectral function for electrons and magnons
The electronic spectral function A(E), or the "interacting density of states" [80,81], is computed using the retarded GF in Eq. (62) as where ∆ k = (E k − E 0 ) and E k are the eigenenergies of quantum many-body HamiltonianĤ [Eq. (4)]. The prefactors of δ-function in A(E) define the "weight" of the many-body eigenstate |Ψ k within A(E). Since the ground state |Ψ 0 is an eigenstate of the electron number operatorN e [Eq. (18)], it has a well-defined number of electrons N e . Thus, the action ofĉ † iσ andĉ iσ on |Ψ 0 in Eq. (65) reveals that the δ-function peaks at E = ±∆ k can only be contributed by those quantum many-body eigenstates |Ψ k which describe systems containing N e ± 1 electrons.
Similarly, the bosonic spectral function D(E) is evaluated using the bosonic retarded GF in Eq. (63) where define the "weight" of many-body eigenstate |Ψ k within D(E). The δ-function peaks in Eq. (66) at E = ±∆ k come from many-body eigenstates |Ψ k . However, unlike the electronic case, they do not have a well-defined total magnon number N mag as they are not eigenstates of the total magnon number operatorN mag = N i=1n i . This is illustrated by Fig. 9(e) with the structure of one selected many-body eigenstate |Ψ k which is a linear combination of many-body states with total magnon number N mag = 0, N mag = 1 and N mag = 2.
Both A(E) and D(E) must satisfy the sum rule This feature allows for physical interpretation where A(E)dE/2N or D(E)dE/N can be viewed as probabilities to find fermion or boson within energy window dE around E in a general quantum many-body system where fermions interact with other fermions and bosons interact with other bosons, as well as with each other. Note that our fermion-boson interacting system, as illustrated in Fig. 1

(b) and described by Hamiltonian in Eq. (4), includes HP bosons interacting [Eq (33)] with other HP bosons when N mag > 1 and electrons interacting with HP bosons while electron-electron interactions are excluded.
Since the sum rule is an exact result, in practical GF calculations it can be employed to test the quality of a verity of analytical and numerical approximations schemes [36].

III. RESULTS AND DISCUSSION
A. Range of validity of truncated HP transformation for nonequilibrium interacting system of magnons Figure 4 compares Ŝ z 1 (t) for 1D chain [ Fig. 1(a)], hosting S = 5/2 or S = 1 localized spins in the absence of electrons (i.e., N e = 0), computed using the original localized spin operators vs. their mapping to bosonic operators via the truncated HP transformation. In the ferromagnetic ground state |Ψ 0 , the expectation value Ŝ z i (t = 0) = Ψ 0 |Ŝ z i |Ψ 0 = 5/2 for all sites i at t = 0 in Fig. 4(a)-(c). To initiate nonequilibrium dynamics for times t > 0, we choose an initial state |Ψ(0) such that the expectation value of the localized spin on site i = 1 is reduced by N mag units, , while on other sites it remains S z i =1 (t = 0) = 5/2. This is equivalent to introducing N mag HP bosons on site i = 1 at t = 0, so that the initial quantum many-body state of HP bosons is given by When N mag = 1, Fig. 4(a) shows that Ŝ z 1 (t), evaluated by truncated HP transformation (green dashed line) solely containing single-particle HamiltonianĤ 0 of noninteracting HP bosons in Eq. (32), accurately tracks the exact time dependence (black lines) of S z 1 (t) evaluated using the localized spin operators, as trivially expected. That is, because there is only one HP boson in the system, magnon-magnon interaction terms active within H int part of the Hamiltonian [Eq. (33)] cannot influence the dynamics of localized spins.
To understand the significance of magnon-magnon interaction terms withinĤ int on the dynamics of localized spins, we next introduce N mag = 2 HP bosons on site i = 1. The time dependence of Ŝ z 1 (t) in Fig. 4(b), evaluated via the truncated HP transformation with truncation number N T = 1 (Sec. II E), matches the exact time evolution obtained using the original localized spin operators only for short enough times (0 < t < 10 fs). At longer times (30 < t < 35 fs), discrepancy emerges due to missing effects from N T > 1 magnon-magnon interaction terms withinĤ int . Thus, to recover the agreement between two types of calculations at longer times requires increasing N T , such as by using N T = 3 (orange dashed lines) in Fig. 4(b). However, progressively larger N T must be employed (Fig. 5) to increase the "breakdown-time" t break [marked in Fig. 4(d),(e)] at which disagreement between two types of calculations emerges. We define t break as the time when the deviation ∆ = Ŝ z 1 (t)| HP − Ŝ z 1 (t)| exact , between Ŝ z 1 (t)| HP (evaluated from truncated HP transformation with a truncation number N T ) and the exact Ŝ z 1 (t)| exact becomes larger than the chosen tolerance |∆| > 10 −4 . As demonstrated by Figs. 4(d),(e) and Fig. 5, t break sensitively depends on the density of HP bosons, N mag /N , whose increase makes magnon-magnon interaction terms withinĤ int more relevant, thereby requiring larger N T in Fig. 5. Figure 4(d),(e) explicitly confirms this conclusion by showing the effect of reduced density of HP bosons on the range of validity of truncated HP transformation, where we employ spin-1 localized spins allowing us to exactly diagonalize larger chains [than those composed of N = 4 sites in Fig. 4(a)-(c) with spin-5 2 on each site]. In Fig. 4(d), at t = 0 we flip the localized spins on sites i = 1 and i = 2 [see inset in Fig. 4(d)], thereby introducing two HP bosons on each of these sites. Thus, the total number of HP bosons within the system in Fig. 4(d) is N mag = 4, whereas the HP boson density is N mag /N = 1. For such parameters, t break ≈ 11 fs (|∆| = 0.05 is chosen solely for visualization of t break at fs time scales). On the other hand, in Fig. 4(e), where HP boson density is reduced to N mag /N = 0.57 by making 1D chain longer from N = 4 sites to N = 7 sites, we find that t break for truncated HP transformation increases to t break ≈ 23 fs. This observation is easily explained since in longer 1D chains the probability for magnon-magnon scattering events is reduced, which makes inclusion of high-order magnon-magnon interaction terms less important and thus the breakdowntime for truncated HP transformation increases. Figure 5 demonstrates how for a given breakdowntime t = t break , the horizontal distance between consecutive curves from left to right increases nonlinearly. This means that N T needed to accurately track Ŝ z 1 (t) via the truncated HP transformation increases nonlinearly with the number of HP bosons N mag excited in the system. On the other hand, if we consider the roughly constant slope 'p' of each curve in Fig. 5 (for the part before a sudden jump), and note the logarithmic scale for the ordinate axis of Fig. 5, we can conclude that t break ∝ exp(pN T ). At first sight, the exponential dependence of t break on N T appears to be favorable i.e., by using larger values of N T (and hence including more and more multi-magnon terms), we can increase t break exponentially and yield accurate dynamics for longer times. However, to obtain a practically tractable MBPT for electron-boson interacting systems [94] a small N T is required but Fig. 5 shows that using small N T = 1-5 allows one to track dynamics of localized spins only up to time t break ≈ 1.5 /J H ≈ 15.0 fs (for J H = 0.1 eV). This is insufficient to model even ultrafast optical manipulation of magnetism requiring simulation times ∼ 10 fs [95], and it is much further away from current-driven magnetization dynamics via spin torque which occurs on ∼ 1 ns time scales [48,50].

B. Range of validity of truncated HP transformation for nonequilibrium interacting system of electrons and magnons
In this Section, we repeat the same analysis as in Sec. III A-but with electron-localized-spin or, equivalently electron-magnon-interaction turned on within 1D quantum many-body system composed of N = 3 sites [ Fig. 1(b)]. These sites host spin-5 2 localized spins interacting with half-filled (N e = 3) tight-binding electrons via the sd exchange interaction [67] of strength J sd as encoded by Eq. (9).
At t = 0, we fully flip the localized spin-5 2 on site i = 1 to initiate nonequilibrium dynamics. From the viewpoint of HP transformation, this is equivalent to introducing N mag = 5 HP bosons on site i = 1, and thus, the initial quantum many-body state is given by in the notation of second-quantization formalism. Here |0 is the vacuum state of electrons and HP bosons combined. Figure 6(a) with J sd = 0 serves as a reference. When electron-magnon interaction is turned on  Fig. 1(b)] comprised of N = 3 sites hosting spin- 5 2 localized spins which interact with conduction electrons. The electron-localized-spin sd exchange interaction is turned off (J sd = 0 eV) in (a) for reference, and turned on in (b) using J sd = 1 eV. The insets in both panels depict the state of localized spins at t = 0, where the localized spin on site i = 1 is completely reversed [i.e., Nmag = 5 HP bosons are created on site i = 1 via Eq. (70)] to initiate nonequilibrium dynamics. When sd interaction is turned on, the disagreement between two types of calculations is actually alleviated when moving from panel (a) to panel (b), which is explained by panel (c) as being due to a rapid loss of total number of magnons, , from the localized spin subsystem to the electronic subsystem. Panels (d) and (e) are counterparts of panels (a) and (b), respectively, for Ŝ z i (t) evaluated using the resummed HP transformation in Eq. (37), which is approximate for Nmax = 3 but it becomes exact for Nmax = 2S = 5.
(J sd = 1 eV) in Fig. 6(b), Ŝ z 1 (t) computed by truncated HP transformation follows the exact result for longer times t 15 fs than in Fig. 6(a). This is explained by Fig. 6(c) which shows that the total number of magnons as a function of time, N mag (t) = N i=1 Ψ(t)|n i |Ψ(t) , is reduced in the course of quantum time evolution. Therefore, this leads to fewer magnon-magnon scattering events which facilitates accurate tracking over longer time intervals of nonequilibrium dynamics of localized spins by truncated HP transformation, in accord with Fig. 5. The lost magnons in Fig. 6(c) are absorbed by the electronic subsystem and mediate transfer of spin angular momentum between the subsystems of electrons and localized spins, while the total z-spin remains conserved [Eq. (21)].
Furthermore, Figs. 6(d) and 6(e), as the counterpart of Figs. 6(a) and 6(b), respectively, demonstrate that electron-boson interacting Hamiltonian can track exact time evolution if truncated HP transformation is replaced by resummed HP transformation in Eq. (37). That is, when N max = 3 is used in Eq. (37), there is disagreement between the two calculations of Ŝ z 1 (t)-compare resummed HP transformation (magenta solid line) vs. the exact one (black solid line)-but increasing N max = 5 in Eq. (37) ensures that both methods match perfectly.
Thus, Figs. 6(d) and 6(e) with properly chosen N max motivate us to derive electron-boson Hamiltonian as the exact mapping of the original electron-localizedspin Hamiltonian in Eq. (4). The former is required for equilibrium or nonequilibrium MBPT [22,23] which can handle [36] systems in two-or three-dimensions composed of large number of sites N 1-that is, the problems where exact diagonalization [65] or (time-dependent) density matrix renormalization group [96][97][98][99] (suitable for N 1 but only in quasi-1D [100]) are inapplicable. Here the terms in Eq. (71) are given by Their physical meaning is transparent:Ĥ e is the tightbinding Hamiltonian of noninteracting electrons;Ĥ mag is the Hamiltonian of noninteracting HP bosons;Ĥ mag−mag describes various interactions between two (first term in H mag−mag ) or more HP bosons; andĤ e−mag describes electron-boson interactions, such as absorption or emission of HP bosons accompanied by electron spin flip as the spin angular momentum is transferred. We note that Eq. (72) is much more complex that what is typically used in spintronics literature [36,39,[41][42][43]45]. Most importantly, it shows that accurate MBPT or diagrammatic Monte Carlo calculations [101,102] in the future for interacting electron-magnon system will have to deal with nonlinear [94] electron-boson interactions.

C. Ground state of interacting system of electrons and magnons
It is also instructive to explore the structure of the exact quantum many-body ground state (GS) of conduction electrons plus localized spins in equilibrium (GS implies zero temperature as well), as described by Hamiltonian in Eq. (4) in terms of the original localized spin operators; as well as to find out how many terms of truncated HP transformation (Sec.II E) need to be retained in order to obtain the same ground state by exact diagonalization of the Hamiltonian of the same system but expressed in terms of electronic and bosonic operators. In this Section, comparison of ground states in two methods is performed for the system depicted in Fig. 1(b) composed of N = 3 sites hosting spin-1 localized spins, while in Secs. III D and III E we also perform comparison of electronic and magnonic spectral functions, respectively, which require additional information about the excited quantum many-body states.
The exact GS is obtained in three steps: (i)Ĥ is represented as a matrix in the basis of eigenstates ofN e andŜ z tot to render a block-diagonal matrix as shown in Fig. 2(c); (ii) to ensure half-filling for electrons, the matrix block corresponding to N e = 3 electrons is isolated; (iii) This matrix block is diagonalized and the eigenstate with the lowest eigenenergy E 0 is identified as the GS. Obviously, if in step (i)Ĥ is expressed directly in terms of localized spin operators [Eq. (28)], then step (iii) yields the numerically exact GS |Ψ 0 . On the other hand, ifĤ is expressed using the truncated HP transformation with a truncation number N T [Eq. (31)], then thus obtained GS |Ψ 0 HP N T is not guaranteed to be the same as |Ψ 0 . In particular, we are interested to know what value of N T ensures that |Ψ 0 HP N T = |Ψ 0 . Figure 7(a) depicts the numerically exact GS |Ψ 0 as a linear combination [88] of many-body kets where red arrows denote quantum state of localized spins (using the same notation as introduced in Sec. II I for spin-1 localized spins) while blue arrows denote spin-↑ or spin-↓ electrons filling three single-particle energy levels 0 = − √ 2 eV, 1 = 0, and 2 = √ 2 eV of noninteracting tight-binding HamiltonianĤ e [Eq. (8)]. In contrast, we find in Fig. 7(b) that GS |Ψ 0 HP N T =6 evaluated using truncated HP transformation with N T = 6 is entirely different from |Ψ 0 shown in Fig. 7(a). Only when the truncation number is increased to N T = 20 in Fig. 7(c) we find |Ψ 0 HP N T =20 ≡ |Ψ 0 . It is worth examining further the structure of the exact GS |Ψ 0 in Fig. 7(a). Its many-body eigenenergy is E 0 = −3.33 eV while the other quantum numbers [Eq. (22)] are N e = 3 and S z tot = 3.5. The largest contribution (greater than 99%) to |Ψ 0 comes from the first term on the RHS in Fig. 7(a) where N e = 3 electrons fill up the single-particle energy levels 0 , 1 , and 2 of noninteracting HamiltonianĤ e [Eq. (8)] in accord with the Pauli exclusion principle while the localized spins are in the ferromagnetic configuration. In the absence of electron-localized-spin interaction, the first term on the RHS would be the only one. Thus, interactions give rise to three states [indicated by horizontal overline in Fig. 7(a)] where N mag = 1 HP boson is created on one of the three sites. This HP boson is actually emitted when the spin-↓ electron in eigenenergy level 0 undergoes a spin-flip process and emerges as a spin-↑ electron in eigenenergy level 2 . This process respects conservation of total z-spin encoded by Eq. (21). The remaining four kets on the RHS of Fig. 7(a) are purely electronic excitations where electrons are excited among eigenenergy levels 0 , 1 , and 2 but are not accompanied by any spin-flip process of localized spins.

D. Electronic spectral function in interacting system of electrons and magnons
For the same system considered in Sec. III C, Fig. 8(a)-(c) compares the electronic spectral function A(E) [Eq. (64)] evaluated from truncated HP transformation vs. the exact one evaluated using localized spin operators. To set a reference point, in Fig. 8(a) we first consider the noninteracting electronic spectral function (black line) when electron-localized-spin interaction is turned off (J sd = 0). For such a case, the available single-particle states are simply the eigenstates of the noninteracting tight-binding HamiltonianĤ e [Eq. (8)] with single-particle energy levels 0 , 1 , and 2 , so that A(E) consists of sharp peaks centered at 0 , 1 , and 2 . Upon turning on electron-localized-spin interaction (J sd = 0.2 eV), the exact A(E) (red line) evaluated using localized spin operators is modified to exhibit peak splitting [with respect to black line reference result when J sd = 0] at energies 0 , 1 , and 2 . Also, few additional peaks around single-particle energy levels 0 and 2 emerge.
In Fig. 8(b), we compute A(E) using truncated HP transformation with a truncation number N T = 6. Although it reproduces the peak-splitting near 0 , 1 and 2 , it exhibits several additional peaks that are absent in the exact result in Fig. 8(a). This discrepancy can be understood as follows. The function A(E) depends on the exact GS |Ψ 0 through Eqs. (64) and (65). However, for N T = 6 Fig. 7 At first sight, it appears that the same argument should produce exact A(E) in Fig. 8(c) using N T = 20 because |Ψ 0 HP N T =20 ≡ |Ψ 0 in Fig. 7(c). However, the discrepancy between exact A(E) (red line) in Fig. 8(a) and blue line in Fig. 8(c) is explained by Eqs. (64) and (65) where A(E) depends both on GS |Ψ 0 and excited many-body states |Ψ k [Eq. (65)] for which truncation number N T = 20 appears to be insufficient. The repeated analysis from Fig. 8(a)-(c) for spin-1 localized spins, but by using spin-5 2 localized spins in Fig. 8(d)-(f), shows that requirement of large N T = 20 cannot be bypassed by increasing the value of localized spins which make them more "classicallike" [3,57,59].  64) for 1D quantum many-body system, illustrated in Fig. 1(b) and comprised of N = 3 sites hosting spin-1 localized spins, is evaluated using localized spin operators in Hamiltonian in Eq. (4). This is compared to approximate A(E) in panels (b) and (c) evaluated by mapping localized spin operators to bosonic operators via truncated HP transformation (Sec. II E) with truncation numbers NT = 6 and NT = 20, respectively. Panels (d)-(f) show counterpart information to panels (a)-(c), but for spin- 5 2 localized spins on each of N = 3 sites. In all panels we set JH = J sd = 0.2 eV, except for black curves in panels (a) and (d) for which electronlocalized-spin interaction is turned off i.e., J sd = 0.

E. Magnonic spectral function in interacting system of electrons and magnons
The exact HP transformation in Eq. (1) makes it possible to define magnonic spectral function D(E) in Eq. (66) and compute it without any approximations by numerically evaluating square root of matrices in Eq. (1). Using the same systems of spin-1 or spin- 5 2 localized spins that are studied in Fig. 8, we first establish a reference magnonic spectral function by computing D(E) in Figs. 9(a) and 9(b) with electron-localized-spin interaction turned off (J sd = 0). Such reference D(E) [ Fig. 9(a)] exhibits three peaks at energies E = 0 eV, E = 0.2 eV, and E = 0.6 eV which correspond to available states in the presence of solely localized-spin-localized-spin (or equivalently magnon-magnon) interactions. Conversely, when we turn on J sd = 0.2 eV in Figs. 9(c) and 9(d), we find: (i) the original noninteracting peaks remain largely intact, except for the one near E = 0.2 eV which undergoes a tiny splitting; (ii) far away the original peaks, D(E) exhibits new additional peaks (marked by dotted circles) near energies E = 1.6 eV and E = 2.1 eV. Analogous features are observed for spin- 5 2 localized spins when switching from J sd = 0 in Fig. 9(b) to J sd = 0 in Fig. 9(d).
We note that similar additional peaks in magnonic spectral function, generated by turning on electronmagnon interaction, were previously observed in MBPT calculations [36] despite being based on resummation of an infinite class of selected diagrams-in contrast, calculations in Figs. 9(c) and 9(d) are nonperturbative and, therefore, correspond to all diagrams being summed to infinite order. These additional peaks in D(E) computed by MBPT were interpreted in Feynmann diagrammatic language as quasibound states of magnons dressed by the cloud of electron-hole pair excitations. Also, MBPT calculations of Refs. [36,87] find much smaller modification of electronic A(E) upon tuning on electron-magnon interaction. This is explained by magnons being in the strongly interacting regime vs. electrons being in the weakly interacting regime due to [36] J sd divided by the bandwidth of noninteracting magnons being much larger than J sd divided by the bandwidth of noninteracting electrons.
To clarify the origin of these peaks further in the context of our exact nonperturbative calculations in Figs. 9(c) and 9(d), we focus on the peak near E = 1.6 eV in Fig. 9(c). This peak is due to manybody excited state |Ψ k whose composition is given explicitly in Fig. 9(e). This state has a nonzero "weight" Q + k = 0.006 in Eq. (67). Although, the value of Q + k appears to be small, it contributes about 2% in the sum rule in Eq. (68b) and thus it cannot be ignored. Interestingly, Fig. 9(e) reveals that this specific |Ψ k is a linear superposition of states with N mag = 0, 1, or 2 HP bosons.

F. Entanglement entropy of ground and excited states of interacting system of electrons and magnons
All three different version of the GS |Ψ 0 in Fig. 7, as well as selected excited state |Ψ k shown in Fig. 9(e), are examples of pure but entangled quantum many-body states [82]. In particular, these states encodes entanglement between electronic and localized-spins subsystems. The von Neumann entanglement entropy [82] for electronic or localized-spins subsystems of the total bipartite system are identical, S e = S lspins , and can be computed from the reduced density matrixρ e where the (improper) mixed quantum state of the electronic subsystem is described by reduced density matrix ρ e = Tr lspins |Ψ Ψ|, which is obtained by partial trace of the pure state density matrix, |Ψ Ψ|, over the basis of states in H lspins . For example, S 0 e = 5.6 × 10 −3 for the exact GS in Fig. 7(a), which means that this many-body entangled state is quite close to separable (characterized by S e ≡ 0) noninteracting (i.e., for J sd = 0) GS as the first term depicted in  66). (e) Composition of specific excited quantum many-body state |Ψ k for spin-1 case which is responsible for peak in D(E) near E = 1.6 eV in panel (c)-this state is a superposition of kets with zero magnons (Nmag = 0) and solely electronic excitations; followed by kets with one magnon excitation (Nmag = 1); and then states with two magnons excitations (Nmag = 2). Numbers on the top of each ket are coefficients in their linear superposition leading to |Ψ k . Fig. 7(a). On the other hand, S 0 e = 0.604 for the GS in Fig. 7(b) which is incorrect [unlike the correct GS in Fig. 7(c) which matches the exact GS in Fig. 7(a)] due to too small N T employed in truncated HP transformation [Sec. II E]. Note that selected excited many-body entangled state |Ψ k analyzed in Fig. 9(e) has much larger S k e = 0.467.

G. Diagonal and off-diagonal elements of time-dependent electronic and magnonic lesser Green functions
The time-dependent electronic lesser GF G < iσ,jσ (t, t ) in Eq. (59) generally depends [22,23] on two time arguments, t and t . At equal times t = t, it yields electronic one-particle nonequilibrium density matrix [22,52,103,104] Its diagonal elements in, e.g., coordinate (or site for discrete lattice) representation contain information about the time-dependent electronic charge and spin density [52], whereas the off-diagonal elements encode quantum-mechanical interference effects [104] and measure the degree of quantum coherence [105]. To illustrate their time evolution, we use the same 1D quantum manybody system employed in Fig. 6 where localized spin-  Conversely, Figs. 10(c) and 10(d) use J sd = 1 eV which leads to nontrivial time dependence of both diagonal and off-diagonal elements of −i G < (t, t). Interestingly, the diagonal elements, −i G < iσ,iσ (t, t) in Fig. 10(c) satisfy −i [G < i↑,i↑ (t, t)+G < i↓,i↓ (t, t)] = Q i with Q i being the total electronic density on site i, are time-independent. This means that no charge currents flows between sites i and j. Instead, population of electrons with spin σ =↑, ↓ on site i exchanges solely between spin σ =↑, ↓ states at that site. This is also accompanied by time evolution of the off-diagonal elements −i G < iσ,jσ (t, t) in Fig. 10(d). The off-diagonal elements of the lesser GF are also required to calculate many-body lesser self-energy Σ < (t 1 , t 2 ) [22,23], which is connected to lesser GF in a self-consistent fashion via the Keldysh equation  tities. Here G a (t, t ) = [G r (t , t)] † is the advanced GF. A self-consistent solution to Eq. (76) can yield exact manybody lesser self-energy. Alternatively, one can systematically approximate it [36] using the so-called "conserving approximations" [106] in MBPT. One such "conserving approximation" for the lesser self-energy of electronboson interacting systems is the so-called self-consistent Born approximation (SCBA) [36,[106][107][108]. The SCBA ensures charge conservation in nonequilibrium [88], and in steady-state nonequilibrium one can Fourier transform Σ < ij (t 1 − t 2 ) to energy domain and operate with Σ < ij (E). To reduce computational complexity [107] of calculations of Σ < ij (E) and enable simulations of devices containing large number of atoms, the "local self-energy" approximation is often employed [83][84][85][86] when modeling inelastic scattering of electrons and bosons. In this approximation, one assumes |Σ < ii (E)| |Σ < ij (E)|, i.e., the off-diagonal elements of self-energy are minuscule when compared to the diagonal ones, and thus, one can set them to zero. This is done in conjunction with discarding the off-diagonal elements of the electronic lesser GF i.e., G < ij (E) ≈ G < ii (E)δ ij . Using our numerically exact electronic lesser GF in Fig. 10(c),(d), we can explicitly check if the "local selfenergy" approximation [83][84][85][86] is warranted for electronmagnon realization of electron-boson quantum manybody system. The off-diagonal elements of the lesser GF in Fig. 10(d) are not minuscule, but are instead approximately one-fifth of the diagonal elements in Fig. 10(c). Therefore, "local self-energy" approximation cannot be justified in the case of many-body electron-magnon in-teracting systems. Figure 10(e)-(h) shows the counterpart of Fig. 10(a)-(d) but for the lesser GF of HP bosons. Here the off-diagonal elements of the bosonic lesser GF, D < (t, t) in Eq. (61), at equal times (i.e., of bosonic oneparticle nonequilibrium density matrix) are always comparable to the diagonal ones independently of whether the electron-localized-spin interaction J sd is turned off

IV. CONCLUSIONS
By applying numerically exact diagonalization techniques to two versions of the Hamiltonian of quantum many-body system of conduction electrons interacting with localized spins that are widely used in spintronics and magnonics, we compare predictions from these two Hamiltonians for: ground state and spectral functions extracted from the retarded GF in equilibrium; and time evolution of the expectation values of localized spin operators and lesser GF in nonequilibrium. The two Hamiltonians, describing systems illustrated in Fig. 1 chosen as 1D and small in order to make calculations tractable, differ in their treatment of localized quantum spins-they are described by either finite-size matrices of the original spin operators or infinite matrices of bosonic operators after the original localized spin operators are mapped to bosonic ones using the popular HP transformation. The truncation [Sec. II E] of HP transformation is always done to make diagrammatic MBPT [22,36] or Monte Carlo [101,102] calculations possible, but mapping of fi-nite size to infinite matrices necessarily requires some approximations which can lead to spurious effects in equilibrium (Fig. 8) or incorrect time evolution (Figs. 4-6) out of equilibrium. Our conclusions are summarized as follows: 1. For quantum many-body systems composed of localized spins alone, Fig. 4 shows that as more interacting HP bosons are introduced into the system, progressively larger number of terms N T is required in truncated HP transformation to incorporate multi-magnon interactions and accurately track the nonequilibrium dynamics of localized spins. Figure 5 shows that the breakdowntime t break for truncated HP transformation follows t break ∝ exp(pN T ). Although, the exponential dependence of t break on the truncation number N T is favorable, the reasonable value of N T = 1-5 typically used in practical calculations does not allow one to track dynamics beyond ∼ 15 fs time scale which is insufficient for ultrafast [95] or spin torque applications [48,50].
2. When electrons are introduced and electronlocalized-spin interaction is turned on, Fig. 6(a)-(c) shows that N T required to accurately track nonequilibrium dynamics of localized spins is actually reduced due to the transfer of spin angular momentum between the two subsystems, which effectively reduces the total number of interacting magnons within the localized spin subsystem. Furthermore, Figs. 6(d) and 6(e) show that the recently introduced [17] resummed HP transformation [Sec. II F] makes it possible to completely evade artifacts of the usual truncated HP transformation. However, the electron-magnon Hamiltonian furnished by it in Eq. (72) is much more complex for MBPT and diagrammatic Monte Carlo calculations than previously used electron-magnon Hamiltonians [36] based on low-order truncated HP transformation.
3. Figure 7 reveals how truncated HP transformation with a small truncation number [such as N T = 6 in Fig. 7(b)] produces an incorrect GS of the interacting electron-magnon system. Only when truncation number is increased [such as to N T = 20 in Fig. 7(c)], exact diagonalization of electron-boson Hamiltonian reproduces the exact GS obtained by diagonalizing the original electron-localized-spinoperators Hamiltonian [ Fig. 7(a)]. However, even large truncation number [such as N T = 20 in Fig. 8(c),(f)] does not ensure that correct electronic spectral function can be obtained from electronboson Hamiltonian due to the fact that spectral functions depends [Eq. (64)] on both the GS and excited quantum many-body states.
4. The magnonic spectral function can be substantially modified [ Fig. 9(c),(d)] upon introduction of conduction electrons and their interaction with localized spins, even when such interaction appears small for electrons, due to much smaller bandwidth of magnons. That is, magnons are effectively pushed into strongly interacting regime, and the new peaks in their spectral function (or "interacting density of states" [80,81]) can be directly related to specific excited quantum many-body states. The structure of excited states [ Fig. 9(a)] reveals superpositions of many-body states in which holes in electronic single particle levels are formed and accompanied by flips of localized spins or, equivalently, creation of one or more virtual HP bosons. 5. The time evolution of the matrix elements of the lesser GF (electronic or magnonic) at equal times in real-space representation, which yields the one-particle nonequilibrium density matrix in realspace representation, shows that the magnitude of the off-diagonal elements is always comparable to the magnitude of the diagonal ones (Fig. 10). Thus, "local self-energy" approximation neglecting the off-diagonal elements, as often employed [83][84][85][86] to enable MBPT modeling of electron-boson systems with large number of atoms, is not warranted.