text:
\section{Introduction}\label{sec:intro}

The concept of spin waves was introduced by Bloch~\cite{Bloch1930} 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 $\lambda$. The quanta of energy of spin waves behave as a quasiparticle termed magnon which carries energy $\hbar \omega$ and spin $\hbar$. As regards terminology, we note that in spintronics and magnonics~\cite{Chumak2015} it is common to use ``spin wave'' for excitations described by the classical Landau-Lifshitz-Gilbert (LLG) equation~\cite{Wieser2015} within numerical micromagnetics~\cite{Kim2010} or atomistic spin dynamics~\cite{Evans2014}, while ``magnon'' is used for quantized version of the same excitation. In other fields, terms ``spin waves'' and ``magnons'' are sometimes used to distinguish between long- and short-wavelength excitations, respectively, or both names are used interchangeably~\cite{Zhitomirsky2013}.

The second quantized description of magnons was introduced by Holstein and Primakoff (HP)~\cite{Holstein1940} by mapping the localized spin operator $\hat{\mathbf{S}}_i$ at site $i$ of the lattice to bosonic operators
\begin{subequations}\label{eq:hp}
	\begin{eqnarray}
	\hat{S}^+_i &=& \hat{S}^x_i + i\hat{S}^y_i
	=  \sqrt{2S}\bigg(1 - \frac{\hat{n}_i}{2S}\bigg)^{1/2}\hat{a}_i, \label{eq:hp_sp}
	\\
	\hat{S}^-_i &=& \hat{S}^x_i - i\hat{S}^y_i
	=  \sqrt{2S} \hat{a}_i^\dagger \bigg(1 - \frac{\hat{n}_i}{2S}\bigg)^{1/2}, \label{eq:hp_sm}
	\\
	\hat{S}^z_i &=&  (S - \hat{n}_i). \label{eq:hp_sz}
	\end{eqnarray}
\end{subequations}
Here $\hat{a}_i^\dagger$ ($\hat{a}_i$) creates (annihilates) HP boson at site $i$ and satisfies the bosonic commutation relations
\begin{equation}\label{eq:hp_comm}
[\hat{a}_i,\hat{a}_j^\dagger] = \id \delta_{ij},
\end{equation}
where $\id$ is the unit operator in the infinite dimensional space of bosons. The HP boson number operator,  \mbox{$\hat{n}_i =\hat{a}_i^\dagger \hat{a}_i$} whose eigenvalues and eigenstates are defined by $\hat{n}_i |n\rangle = n |n\rangle$, measures how much the localized spin deviates away from the ground state, where we assume the ferromagnetic ground state with the $z$-axis as the easy axis. Thus, the creation of one HP boson is equivalent to the removal of one unit of spin angular momentum from the ground state [see Fig.~\ref{fig:fig1}(c) for illustration and Sec.~\ref{sec:magnons} for details].

Since bosonic operators $\hat{a}_i^\dagger$, $\hat{a}_i$  act on an {\em infinite}-dimensional Hilbert space, but the physical Hilbert space corresponding to
a single site is spanned by only $2S + 1$ states, the extra unphysical states are decoupled from the physical ones by the square root in Eq.~\eqref{eq:hp}. However, the  textbook~\cite{Mahan2011,Chudnovsky2006} and research literature are typically focused on band structure of noninteracting magnons (which can also be topologically nontrivial~\cite{Kim2016}), so they employ the lowest-order truncation
\begin{subequations}\label{eq:approxhp}
\begin{eqnarray}
\hat{S}^+_i  & \approx & \sqrt{2S}  \hat{a}_i, \\
\hat{S}^-_i & \approx  & \sqrt{2S} \hat{a}_i^\dagger,
\end{eqnarray}
\end{subequations}
of original HP transformation in Eq.~\eqref{eq:hp}, and then retaining only the terms in the Hamiltonian that are up to the quadratic order in the bosonic operators. This assumes low-density limit, \mbox{$\langle \hat{n}_i \rangle/2S \ll 1$}, achieved at low temperatures and/or large $S \gg 1$ in which HP bosons can be treated as noninteracting.

Taking into account higher order terms in the Taylor expansion of the square root in Eq.~\eqref{eq:hp} produces higher-than-quadratic terms in the Hamiltonian  which describe {\em boson-boson interactions}~\cite{Chudnovsky2006,Tupitsyn2008,Yuan2020,Takei2019}. Note the exact HP transformation in Eq.~\eqref{eq:hp} splits the infinite dimensional Hilbert space spanned by boson number states $\{ |n\rangle \}_{n \in \mathbb{N}}$ into two sectors---physical states $\{ |n\rangle \}_{n=0,...,2S}$; and all the unphysical ones $\{ |n\rangle \}_{n>2S}$ [see also Eq.~\eqref{eq:phys_states}]. Those sectors cannot be connected by $\hat{S}^+_i$ and $\hat{S}^-_i$ operators. However, when the square root in Eq.~\eqref{eq:hp} is expanded in power series and then truncated (see Secs.~\ref{sec:hptruncated} and~\ref{sec:hpresum}) to any {\em finite} order $N_T$, the physical and unphysical subspaces become coupled. In addition, canonical commutation relations for the spin operators are only satisfied approximately, resulting in artificial breaking of rotational symmetries that may be present in the original Hamiltonian~\cite{Vogl2020,Konig2020}.

These higher order terms of the Taylor series expansion  are required to study equilibrium properties at increasing temperature~\cite{Dyson1956,Hofmann2011} or  in nonequilibrium~\cite{Tupitsyn2008,Yuan2020,Schuckert2018,Takei2019}. They describe magnon-magnon interaction~\cite{Radosevic2015}, which  renormalize magnon energy and lead to magnon decay over a finite lifetime~\cite{Zhitomirsky2013}. For example, Dyson~\cite{Dyson1956} 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 + \mathcal{O}(T^{9/2})$, where $T^{3/2}$ is the so-called Bloch law for noninteracting magnons with parabolic energy-momentum dispersion; second and third term also stem from noninteracting magnons but with nonparabolic dispersion on a discrete lattice; and magnon-magnon interaction starts manifesting at order $T^4$. Such calculations require many-body perturbation theory (MBPT)~\cite{Stefanucci2013,Schlunzen2020} which is most easily formulated in terms of bosonic or fermionic operators. For such  operators, the Wick theorem~\cite{Leeuwen2012} 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 $\mathcal{O}(g^n)$ ($g$ is the parameter of many-body interaction strength). 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.~\eqref{eq:hp}, conjectured by Kubo~\cite{Kubo1953} to be  only an asymptotic series, lead to cumbersome MBPT~\cite{Harris1971,Hamer1992}, a plethora of other mappings of original localized spin operators to bosonic or fermionic operators has been proposed. This includes mappings to Dyson-Maleev bosons~\cite{Dyson1956,Maleev1958}, Schwinger bosons~\cite{Schuckert2018}, fermions~\cite{Jordan1928,Affleck1998}, Majorana fermions~\cite{Tsvelik1992}, sypersymmetric operators~\cite{Coleman2000} and even exotic particles called semions~\cite{Kiselev2000}. We note that mapping of localized spin operators to standard bosonic or fermionic operators could be evaded in path integral formulation by using spin coherent states, but that leads to topological terms associated with the Berry phase, so that even in path integral formalism mapping to bosonic operators is preferred~\cite{Schuckert2018}.
\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{fig1.pdf}
	\caption{(a) Schematic view of a finite-size quantum many-body system comprised of 1D chain of $N=4$ sites hosting spin-$S$ localized spins ($S=1$ or $\frac{5}{2}$) interacting with each other via nearest-neighbor exchange interaction of strength $J_\mathrm{H}$; at $t=0$ nonequilibrium dynamics is initiated by flipping the localized spin at site $i=1$. The quantum many-body system in (b) is same as (a) but with $N=3$ sites hosting \mbox{spin-$S$} localized spins that can also interact with conduction electrons spin (blue arrows) via the $sd$ exchange interaction of strength $J_\mathrm{sd}$ that turns-on electron-magnon interaction. The electrons can hop along the 1D chain with a hopping potential \mbox{$\gamma = 1$ eV} where a ``half-filling" i.e., $N_\mathrm{e} = 3$ is used. (c) Cartoon showing that reducing the expectation value $S^z_i = \braket{\hat{S}^z_i}$ of a spin-$\frac{5}{2}$ localized spin at site $i=1$ by $N_\mathrm{mag}$ units is equal to creating $N_\mathrm{mag}$ number of Holstein-Primakoff bosons at site $i=1$.}
	\label{fig:fig1}
\end{figure}

However, although the Dyson-Maleev transformation evades usage of the square root of operators in Eq.~\eqref{eq:hp}, it generates Hamiltonian that is {\em no longer manifestly Hermitian}. The Schwinger transformation require the introduction of auxiliary fields. The Jordan-Wigner transformation~\cite{Jordan1928} to fermionic operators is exact, but it works only for spin-$\frac{1}{2}$ operators. These disadvantaged have prompted a very recent revival of the HP transformation by considering possible nonperturbative replacement of the Taylor expansion of the square root in Eq.~\eqref{eq:hp} which can be written as the polynomial expression in bosonic operators that does not couple physical and unphysical subspaces while generating {\em manifestly Hermitian} Hamiltonian~\cite{Vogl2020,Konig2020}.  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~\cite{Vogl2020}.

The MBPT calculations of equilibrium magnon properties based on Dyson-Maleev vs. truncated Taylor series expansion of the HP transformation have been carefully compared in the literature over many decades~\cite{Harris1971,Hamer1992}. On the other hand, much less is know about the range of validity of truncated HP transformation when describing nonequilibrium dynamics of localized spins~\cite{Marcuzzi2016}, {\em especially in the presence of additional interactions with conduction electrons}. This is precisely the situation frequently encountered in transport phenomena within spintronics. The nonequilibrium MBPT~\cite{Stefanucci2013,Schlunzen2020} for such problems is {\em virtually always} conducted using HP transformation, but {\em truncated}, as exemplified by theoretical and computational modeling of inelastic electron tunneling spectroscopy in magnetic tunnel junctions~\cite{Mahfouzi2014}; spin-transfer~\cite{Tay2013,Cheng2019,Bender2019} and spin-orbit torques~\cite{Okuma2017}; ultrafast demagnetization~\cite{Tveten2015}; and conversion of  magnonic spin currents into electronic spin current or vice versa at magnetic-insulator/normal-metal interfaces~\cite{Zheng2017,Troncoso2019}. Similarly, truncated HP transformation is chosen mapping from localized spin operators to bosonic operators for problems of quantum magnonics focused on nonequilibrium dynamics of localized spins within magnetic insulators which interact with each other~\cite{Tupitsyn2008,Yuan2020,Takei2019,Kamra2016} or with external electromagnetic fields~\cite{Parvini2020}. One can expect that any truncated HP transformation will break down at sufficiently long times where higher-order terms in the expansion of square root in HP transformation become important, introducing interactions among bosons.  The breakdown then precludes 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 torques~\cite{Tay2013}. Surprisingly, {\em rigorous analysis of such breakdown is lacking}.

Instead, current-driven magnetization reversal due to spin torque~\cite{Ralph2008} is standardly modeled by the LLG equation~\cite{Berkov2008}, which can be combined in  multiscale  fashion with steady-state~\cite{Ellis2017} or time-dependent quantum transport calculations~\cite{Petrovic2018,Bajpai2019a,Suresh2020,Suresh2021,Bajpai2020,Stahl2017,Bostrom2019} considering single-particle quantum Hamiltonians for electrons. However, such hybrid quantum-classical theories, which treat electrons quantum-mechanically and employ classical dynamics of localized spins via the LLG equation, are justified only in the limit of very large localized spins $S \rightarrow \infty$ and $\hbar \rightarrow 0$ (while $S \times \hbar  \rightarrow 1$)~\cite{Wieser2015,Stahl2017}, as well as in the absence of entanglement~\cite{Wieser2015,Mondal2019,Petrovic2020} between quantum states of localized spins. For example, in the emerging concept of quantum spin torque~\cite{Mondal2019,Petrovic2020,Petrovic2020a,Mitrofanov2020,Mitrofanov2021}, describing transfer of angular momentum between spins of flowing electrons and localized spins in situations~\cite{Zholud2017} where the latter {\em must} be described by quantum-mechanical operators, the whole system of electrons and localized spins must be described by a quantum many-body Hamiltonian [as exemplified by Eqs.~\eqref{eq:ham_sys} and ~\eqref{eq:em_ham}].

In this study, we use exact diagonalization techniques~\cite{Wang2019} applied to quantum many-body Hamiltonians defined on chains of  $N \le 7$ sites hosting fermionic (for electrons) and localized spin operators, or fermionic and bosonic (obtained by mapping the original localized spin operators) operators. Comparing their many-body quantum states and Green functions (GFs) in equilibrium and nonequilibrium makes it possible to precisely delineate {\em the range of validity of truncated HP transformation}. We consider a one-dimensional (1D) quantum Heisenberg ferromagnetic spin-$S$ XXX chain hosting localized spins which interact via the nearest-neighbor exchange interaction of strength $J_\mathrm{H}$, as illustrated in Fig.~\ref{fig:fig1}(a). We use both spin $S=1$ as the ultraquantum limit and $S=5/2$. Naively, the eigenvalue of $\hat{\mathbf{S}}_i^2$ being $S^2(1+1/S)$ suggests that quantum effects become progressively less important for $S>1$, but they exist for all $S < \infty$ vanishing as $1/2S$ in the classical limit~\cite{Parkinson1985}. The nonzero electron hopping $\gamma$ between chain of  $N=3$ sites, chosen when electrons are present and illustrated in Fig.~\ref{fig:fig1}(b), means that such 1D chain models a ferromagnetic metal (FM) where $N_\mathrm{e}$  conduction electron [we consider half filled lattice, so $N_\mathrm{e}=3$ for systems in Fig.~\ref{fig:fig1}(b)] spins interact with localized spins via $sd$ exchange interaction~\cite{Cooper1967} usually considered in spintronics. From the viewpoint of the physics of strongly correlated electrons, the model illustrated in Fig.~\ref{fig:fig1}(b) is interpreted as the Kondo-Heisenberg chain~\cite{Tsvelik2017}.

For technical reasons (i.e., exponential increase of the Hilbert space size), we consider 1D chains of $N \le 7$ sites while concentrating on {\em 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~\cite{Spinelli2014} or antiferromagnetically~\cite{Loth2012} coupled 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~\cite{Spinelli2014}. 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~\cite{Schumann2010,Carrascal2015,Hermanns2014} described by the ``pure'' and extended Hubbard models~\cite{Schumann2010}, respectively; or electrons interacting with bosons~\cite{Sakkinen2015,Sakkinen2015a,Dimitrov2017}---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-molecular-optical physics.

The paper is organized as follows. In Sec.~\ref{sec:theory} we introduce different versions of quantum many-body Hamiltonian describing systems in Fig.~\ref{fig:fig1}  and their matrix representations in the Fock space, as well as a procedure to obtain the exact one-particle retarded and lesser GFs that  depend on two times in general. In particular, its subsection~\ref{sec:hptruncated}  introduces an infinite power series expansion of the HP transformation from Eq.~\eqref{eq:hp} and its truncation, while its subsection~\ref{sec:hpresum} provides a brief summary of the very recently proposed~\cite{Vogl2020}  resummed HP transformation. The time evolution of quantum many-body states of a spin chain with no electrons ($N_\mathrm{e}=0$) is employed in Sec.~\ref{sec:validity_mm} to examine the range of validity of truncated HP transformation when tracking time-evolution of localized spins in the presence of \emph{magnon-magnon} interaction and different magnon number $N_\mathrm{mag}$ as a fraction of the maximum possible number. Then in Sec.~\ref{sec:validity_em} we introduce electrons into the 1D chain to examine the range of validity of truncated HP transformation when tracking time evolution of localized spins in the presence of both \emph{electron-magnon} and \emph{magnon-magnon} interactions. In the same Sec.~\ref{sec:validity_em}, we additionally employ the resummed HP transformation~\cite{Vogl2020} to derive quantum many-body Hamiltonian [Eq.~\eqref{eq:em_ham}] for electron-magnon systems in terms of fermionic and bosonic operators whose usage reproduces numerically exact result from calculations based on original localized spin operators. In Secs.~\ref{sec:groundstate} and ~\ref{sec:spectralelectron} we compare ground state and electronic spectral function (or ``interacting density of states''~\cite{Balzer2011,Nocera2018}) of many-body Hamiltonian using the original localized spin operators vs. Hamiltonian using bosons from truncated HP transformation. The excited states and magnonic spectral function are  studied in Sec.~\ref{sec:spectralmagnon}. Finally, Sec.~\ref{sec:diag_offdiag} 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''~\cite{Luiser2009,Rhyner2014,Cavassilas2016,Bescond2018} for electron-boson interacting systems, neglecting the off-diagonal elements, is generally not justified. We conclude in Sec.~\ref{sec:conclusions}.

 \section{Models and Methods}\label{sec:theory}
\subsection{Quantum many-body Hamiltonian of electrons interacting with localized spins}\label{sec:hamilt}
The many-body Hamiltonian of 1D chain composed of $N$ sites, hosting spin-$S$ localized spins and interacting with conduction electrons as illustrated in Fig.~\ref{fig:fig1}(b), is given by~\cite{Woolsey1970}
\begin{equation}\label{eq:ham_sys}
\hat{H} = \hat{H}_\mathrm{e} + \hat{H}_\mathrm{lspins} + \hat{H}_\mathrm{e-lspins}.
\end{equation}
It acts in the total space $\mathcal{F}_\mathrm{e} \otimes \mathcal{H}_\mathrm{lspins}$ which is tensor product of the Fock space of electrons, $\mathcal{F}_\mathrm{e}$, and the Hilbert space of all localized spins
\begin{equation}\label{eq:hlspins}
\mathcal{H}_\mathrm{lspins} = \mathcal{H}_1 \otimes \cdots \otimes \mathcal{H}_N.
\end{equation}
The Fock space of electrons~\cite{Schlunzen2020}
\begin{equation}\label{eq:fock}
   \mathcal{F}_\mathrm{e}= \overline{\mathbb{C} \oplus  \mathcal{H}_\mathrm{e} \oplus \hat{\mathcal{A}}(\mathcal{H}_\mathrm{e} \otimes \mathcal{H}_\mathrm{e})  \oplus  \hat{\mathcal{A}}(\mathcal{H}_\mathrm{e} \otimes \mathcal{H}_\mathrm{e}  \otimes \mathcal{H}_\mathrm{e}) \oplus \cdots},
\end{equation}
is induced by the one-electron Hilbert space $\mathcal{H}_\mathrm{e}$ as the completion (indicated by overline) of the direct sum of
antisymmetrized $n$-fold tensor products of ${\mathcal{H}_\mathrm{e}}$. The operator $\hat{\mathcal{A}}$ antisymmetrizes tensors for fermionic particles. In the sector of $\mathcal{F}_\mathrm{e} \otimes \mathcal{H}_\mathrm{lspins}$ with $N_\mathrm{e}=0$ electrons, we have a chain hosting only  spin-$S$ localized spins, as illustrated in Fig.~\ref{fig:fig1}(a), which is described solely by $\hat{H}_\mathrm{lspins}$ term
\begin{equation}\label{eq:spin_ham}
\hat{H}_\mathrm{lspins} = -J_\mathrm{H} \sum_{\braket{ij}} \hat{\bold{S}}_i\cdot \hat{\bold{S}}_j,
\end{equation}
chosen as the quantum Heisenberg Hamiltonian with nearest-neighbor (NN) exchange interaction (as signified by $\braket{ij}$ notation) of strength \mbox{$J_\mathrm{H}=1$ eV}. When electrons are present, they are described by $\hat{H}_\mathrm{e}$ term
\begin{equation} \label{eq:elec_ham}
	\hat{H}_\mathrm{e} = -\gamma \sum_{\braket{ij}} \hat{\psi}_i^\dagger \hat{\psi}_j,
\end{equation}
chosen as the tight-binding Hamiltonian with the NN hopping \mbox{$\gamma=1$ eV} between single $s$-orbitals residing on each site. The Hamiltonian describing $sd$ exchange interaction of strength \mbox{$J_\mathrm{sd}=0.2$ eV}~\cite{Cooper1967} between conduction electron spin and localized spins is given by
\begin{equation}\label{eq:sd_ham}
\hat{H}_\mathrm{e-lspins} = -J_\mathrm{sd} \sum_{i=1}^{N} \hat{\psi}_i^\dagger \hat{\bm \sigma} \hat{\psi}_i  \cdot \hat{\bold{S}}_i.
\end{equation}
The row vector operator $\hat{\psi}_i^\dagger = (\psi_{i\uparrow}^\dagger,\psi_{i\downarrow}^\dagger )$ consists of operators $\hat{\psi}_{i\sigma}^\dagger$ which create an electron of spin $\sigma = \uparrow, \downarrow$ at site $i$; $\hat{\psi}_i$ is a column vector operator that contains the corresponding annihilation operators; and  $\hat{\bm \sigma} = (\hat{\sigma}^x,\hat{\sigma}^y,\hat{\sigma}^z)$ is the vector of the Pauli spin matrices describing spin-$\frac{1}{2}$ of conduction electron.

Using notation $\{\hat{O}_1,\hat{O}_2\}$ for the anticommutator and $[\hat{O}_1,\hat{O}_2]$ for the commutator of two operators $\hat{O}_1$ and $\hat{O}_2$, fermionic operators for electrons satisfy
\begin{equation}\label{eq:fermi_comm}
\{ \hat{\psi}_{i\sigma},\hat{\psi}^\dagger_{j\sigma'}\} = \id \delta_{ij}\delta_{\sigma\sigma'},
\end{equation}
where $\id$ is a $4^N \times 4^N$ identity matrix in the antisymmetrized $2N$-particle subspace of $\mathcal{F}_\mathrm{e}$ Fock-space of electrons. Additionally, the localized spin operators, $\hat{S}^\alpha_i$ ($\alpha = x,y,z$), at site $i$ satisfy the angular momentum algebra
\begin{subequations}\label{eq:ang_al}
\begin{eqnarray}
&[\hat{S}^x_a,\hat{S}^y_b] = i \hat{S}_b^z\delta_{ab}, \label{eq:ang_al_a}
\\
&[\hat{S}^y_a,\hat{S}^z_b] = i \hat{S}_b^x\delta_{ab}, \label{eq:ang_al_b}
\\
&[\hat{S}^z_a,\hat{S}^x_b] = i \hat{S}_b^y\delta_{ab}. \label{eq:ang_al_c}
\end{eqnarray}
\end{subequations}
The square of the localized spin operator,  $\hat{S}^2_i = (\hat{S}^x_i)^2 + (\hat{S}^y_i)^2 + (\hat{S}^z_i)^2$, commutes with each component
\begin{equation}\label{eq:tot_ang}
[\hat{S}^2_i,\hat{S}^\alpha_j] = 0.
\end{equation}
For computational convenience in calculations of electronic GFs, we change the basis of one-particle electronic states from site basis to eigenenergy basis to obtain
\begin{equation}\label{eq:eham_diag}
\hat{H}_\mathrm{e} = \sum_{i=1}^{N} \epsilon_{i}\hat{c}^\dagger_{i} \hat{c}_{i},
\end{equation}
where $c_i^\dagger = (c_{i\uparrow}^\dagger,c_{i\downarrow}^\dagger )$ is a row vector consisting of $c^\dagger_{i\sigma}$ operators which create an electron in the one-particle electronic eigenstate $\ket{\epsilon_i}$ of eigenenergy $\epsilon_{i}$, such that \mbox{$\hat{H}_\mathrm{e} \ket{\epsilon_{i}} = \epsilon_{i}\ket{\epsilon_{i}}$}, with a spin $\sigma = \uparrow, \downarrow$. These eigenenergies and eigestates are evaluated by diagonalizing one-particle tight-binding Hamiltonian
\begin{equation}\label{eq:ham_1ele}
\hat{H}_\mathrm{e} = \sum_{\braket{ij}}-\gamma \ket{i}\bra{j},
\end{equation}
where $\ket{i}$ denotes $s$-orbital [whose coordinate representation is $\langle \mathbf{r}| i \rangle = \phi(\mathbf{r}-\mathbf{R}_i)$] of an electron centered at site $i$. Using change of basis transformation rules for operators in second quantization
\begin{equation}\label{eq:cc_dag_trans}
\hat{\psi}_{i\sigma} = \sum_{j=1}^N \braket{i|\epsilon_j}\hat{c}_{j\sigma},
\end{equation}
and substituting this into Eq.~\eqref{eq:sd_ham} we get
\begin{equation}\label{eq:sd_new_bas}
\hat{H}_\mathrm{e-lspins} = -J_\mathrm{sd}
\sum_{i}^N\sum_{j}^N\sum_{j'}^N
\braket{\epsilon_{j'}|i}\braket{i|\epsilon_j}
\hat{c}_{j'}^\dagger \hat{\bm \sigma} \hat{c}_j \cdot \hat{\mathbf{S}}_i.
\end{equation}
Since each $\hat{c}_i^\dagger$ or $\hat{c}_i$ operator is represented by $4 \times 4$ matrix (see Sec.~\ref{sec:matrixofccdagger}), and each $\hat{\mathbf{S}}_i$ operator is represented by a $(2S+1) \times (2S+1)$ matrix, the matrix representation of many-body Hamiltonian in Eq.~\eqref{eq:ham_sys} for the chain of $N$ sites in Fig.~\ref{fig:fig1}(b) has size $4^N \times (2S+1)^N$. Although larger than our choice of $N=3$ or $N=4$ sites could be diagonalized with state-of-the-art techniques, we restrict our analysis to such small number of sites in order to make the population of energy levels $\epsilon_i$ for a chosen quantum many-body state easy to analyze.

\subsection{Symmetries of many-body Hamiltonian}
\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{fig2.pdf}
    \caption{Visualization of the absolute value of matrix elements of quantum many-body Hamiltonian $\hat{H}$ [Eq.~\eqref{eq:ham_sys}] for 1D chain in Fig.~\ref{fig:fig1}(b) of $N=3$ sites hosting spin-$1$ localized spins and variable number of electrons $N_\mathrm{e}$ hopping between these   sites while interacting with localized spins. The matrix representation of $\hat{H}$ is in: (a) site basis; (b) basis composed of eigenstates of electron number operator $\hat{N}_\mathrm{e}$ [Eq.~\eqref{eq:ele_num}], where it is block-diagonal and each block consists of states with fixed number of electrons $N_\mathrm{e}$ ($N_\mathrm{e} = 2, 3, 4$ is marked explicitly); and (c) basis composed of eigenstates of \emph{both} $\hat{N}_\mathrm{e}$  [Eq.~\eqref{eq:ele_num}] and total $z$-spin operator $\hat{S}^z_\mathrm{tot}$  [Eq.~\eqref{eq:total_sz}].}
    \label{fig:fig2}
\end{figure*}
The exact diagonalization of many-body Hamiltonian in Eq.~\eqref{eq:ham_sys}
\begin{equation}\label{eq:ham_spec}
\hat{H} \ket{\Psi_k} = E_k \ket{\Psi_k},
\end{equation}
yields its many-body eigenenergies $E_k$ and many-body eigenstates $\ket{\Psi_k}$. The total electron number operator is given by
\begin{equation}\label{eq:ele_num}
\hat{N}_\mathrm{e} = \sum_{i=1}^{N} (\hat{c}_{i\uparrow}^\dagger \hat{c}_{i\uparrow} + \hat{c}_{i\downarrow}^\dagger \hat{c}_{i\downarrow}).
\end{equation}
The operator of total spin in the $z$-direction
\begin{equation}\label{eq:total_sz}
\hat{S}^z_\mathrm{tot} = \frac{1}{2}\sum_{i=1}^N(\hat{c}_{i\uparrow}^\dagger \hat{c}_{i\uparrow} - \hat{c}_{i\downarrow}^\dagger \hat{c}_{i\downarrow}) + \sum_{i=1}^N \hat{S}^z_i,
\end{equation}
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.~\eqref{eq:ham_sys} has two symmetries encoded in the commutation relations
\begin{equation}\label{eq:num_sym}
[\hat{H},\hat{N}_\mathrm{e}] = 0,
\end{equation}
which is due to conservation of number of electrons $N_\mathrm{e}$, and
\begin{equation}\label{eq:sz_sym}
[\hat{H},\hat{S}_\mathrm{tot}^z] = 0,
\end{equation}
which is due to conservation of total $z$-spin (electronic + localized spin) $S^z_\mathrm{tot}$. Therefore, $N_\mathrm{e}$ and $S^z_\mathrm{tot}$ serve as ``good quantum numbers" for labeling many-body eigenstates
\begin{equation}\label{eq:qnumbers}
\ket{\Psi_k} = \ket{E_k, N_\mathrm{e}, S^z_\mathrm{tot}},
\end{equation}
together with many-body eigenenergy $E_k$. The two symmetries in Eq.~\eqref{eq:num_sym} and Eq.~\eqref{eq:sz_sym} can be visualized in the matrix-representation (see Secs.~\ref{sec:matrixofccdagger} and~\ref{sec:exact}) of many-body Hamiltonian $\hat{H}$~[Eq.~\eqref{eq:ham_sys}].

For example, in Fig.~\ref{fig:fig2}(a)  the matrix elements of $\hat{H}$ in the original site basis are visually scattered throughout the whole matrix. However,  when $\hat{H}$ is represented in the basis of eigenstates of $\hat{N}_\mathrm{e}$, Fig.~\ref{fig:fig2}(b) shows that its matrix becomes block-diagonal where each block containing the nonzero matrix elements is associated with states with fixed number of electrons $N_\mathrm{e}$. Finally, in Fig.~\ref{fig:fig2}(c) $\hat{H}$ is represented in the basis of eigenstates $\ket{E_k, N_\mathrm{e}, S^z_\mathrm{tot}}$ of \emph{both} $\hat{N}_\mathrm{e}$ and $\hat{S}^z_\mathrm{tot}$ which isolates additional submatrices with fixed $S^z_\mathrm{tot}$  within blocks associated with fixed $N_\mathrm{e}$.

\subsection{Matrix representation of electronic creation and annihilation operators}\label{sec:matrixofccdagger}
A fermionic operator creating or annihilating electrons on a single site operate within the natural basis of kets  $\ket{0}$, $\ket{\uparrow}$, $\ket{\downarrow}$ and $\ket{\uparrow\downarrow}$ which denote the empty state; state with one \mbox{spin-$\uparrow$} electron; state with one \mbox{spin-$\downarrow$} electron; and the state with one \mbox{spin-$\uparrow$} and one \mbox{spin-$\downarrow$} electron. Thus, these basis states are represented by column vectors
\begin{equation}\label{eq:electron_basis}
\ket{0} =
\begin{pmatrix}
1 \\
0 \\
0 \\
0 \\
\end{pmatrix} \hspace{0.1cm},\hspace{0.1cm}
\ket{\uparrow} =
\begin{pmatrix}
0 \\
1 \\
0 \\
0 \\
\end{pmatrix} \hspace{0.1cm},\hspace{0.1cm}
\ket{\downarrow} =
\begin{pmatrix}
0 \\
0 \\
1 \\
0 \\
\end{pmatrix} \hspace{0.1cm},\hspace{0.1cm}
\ket{\uparrow \downarrow} =
\begin{pmatrix}
0 \\
0 \\
0 \\
1 \\
\end{pmatrix}.
\end{equation}
In the same basis, creation and annihilation operators that act in the 1-site or 2-particle subspace \mbox{$\mathscr{H}_{N=1} = \hat{\mathcal{A}}(\mathcal{H}_\mathrm{e} \otimes \mathcal{H}_\mathrm{e})$} of the Fock space $\mathcal{F}_\mathrm{e}$ are represented by $4 \times 4$ matrices
\begin{equation}
\hat{\psi}^\dagger_{\uparrow} = \begin{pmatrix}
0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0
\end{pmatrix} \hspace{0.1cm},\hspace{0.1cm}
\hat{\psi}^\dagger_{\downarrow} = \begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 \\
0 & -1 & 0 & 0
\end{pmatrix} \hspace{0.1cm},\hspace{0.1cm}
\hat{\psi}_{\sigma} = [\hat{\psi}^\dagger_{1\sigma}]^T,
\end{equation}
which satisfy the fermionic commutation relations in Eq.~\eqref{eq:fermi_comm}.
If we consider two sites, then electronic creation (annihilation) operators,  $\hat{\psi}_1^\dagger$ ($\hat{\psi}_1$) and $\hat{\psi}_2^\dagger$ ($\hat{\psi}_2$), act in the 2-site or 4-particle subspace \mbox{$\mathscr{H}_{N=2} = \hat{\mathcal{A}}(\mathcal{H}_\mathrm{e} \otimes \mathcal{H}_\mathrm{e} \otimes \mathcal{H}_\mathrm{e} \otimes \mathcal{H}_\mathrm{e})$} of the Fock space $\mathcal{F}_\mathrm{e}$  and are represented by matrices of size $4^2 \times 4^2$. For example, the action of $\hat{\psi}^\dagger_{1\sigma}$ in $\mathscr{H}_{N=2}$ is given by
\begin{equation}
(\hat{\psi}^\dagger_{1\sigma})_{\mathscr{H}_{N=2}} = \hat{\psi}^\dagger_{\sigma}\otimes \id,
\end{equation}
where $\id$ is a 4 $\times$ 4 unit matrix. However, the action of $(\psi^\dagger_{2\sigma})_{\mathscr{H}_{N=2}}$ is given by~\cite{Frederiksen2004}
\begin{equation}\label{eq:pmatrix}
(\psi^\dagger_{2\sigma})_{\mathscr{H}_{N=2}} = \hat{P} \otimes \hat{\psi}^\dagger_\sigma.
\end{equation}
Here the permutation matrix $\hat{P}=\text{diag}(1, -1, -1, 1)$ must be included, instead of na\"ive usage of the unit matrix $\id$, in order to preserve the correct anticommutation relations of fermionic operators at different sites in Eq.~\eqref{eq:fermi_comm}. The next step is to construct the matrix representations of  electronic creation and annihilation operators for three sites, which is done in a similar fashion
\begin{subequations}\label{eq:ele_create}
\begin{eqnarray}
(\hat{\psi}_{1\sigma}^\dagger)_{\mathscr{H}_{N=3}} &=& \psi^\dagger_\sigma \otimes \id \otimes \id, \\
(\hat{\psi}_{2\sigma}^\dagger)_{\mathscr{H}_{N=3}} &=& \hat{P} \otimes \hat{\psi}^\dagger_\sigma \otimes \id, \\
(\hat{\psi}_{3\sigma}^\dagger)_{\mathscr{H}_{N=3}}  &=&  \hat{P} \otimes \hat{P} \otimes \hat{\psi}^\dagger_\sigma,
\end{eqnarray}
\end{subequations}
where each operator above is a matrix of size $4^3 \times 4^3$. This also makes it clear how to construct inductively matrix representations of electronic creation and annihilation operators for arbitrary number of sites $N$, where these operators act in  $\mathscr{H}_N = \hat{\mathcal{A}}(\underbrace{\mathcal{H}_\mathrm{e} \otimes \mathcal{H}_\mathrm{e}  \cdots \mathcal{H}_\mathrm{e} \otimes \mathcal{H}_\mathrm{e}}_{2N~\text{times}})$ subspace of the Fock space $\mathcal{F}_\mathrm{e}$.

\subsection{Localized spin operators}\label{sec:exact}
The matrix representation of the localized spin operator $\hat{\bold{S}}=(\hat{S}^x,\hat{S}^y,\hat{S}^z)$ is given by
\begin{subequations}\label{eq:spin_mat}
\begin{eqnarray}
\bra{m\pm1}\hat{S}^x\ket{m} &=& \frac{1}{2}\sqrt{S(S+1) - m(m\pm1)}, \label{eq:sx_mat}
\\
\bra{m\pm1}\hat{S}^y\ket{m} &=& \frac{1}{2i} \sqrt{S(S+1) - m(m\pm1)}, \label{eq:sy_mat}
\\
\bra{m}\hat{S}^z\ket{m} &=& m , \label{eq:sz_mat}
\end{eqnarray}
\end{subequations}
where $\ket{m}$ is an eigenstate of $\hat{S}^z$; $m \in \{-S,-S+1,\cdots,S-1,S\}$; and $\hat{S}^\alpha$ is a $(2S+1)\times(2S+1)$ matrix that acts in the single-site subspace $\mathcal{H}_i$ of $\mathcal{H}_\mathrm{lspins}$. For the chain in Fig.~\ref{fig:fig1}(a) composed of $N$ sites hosting spin-$S$ localized spins, their operators act in the total space of all localized spins, $\mathcal{H}_\mathrm{lspins}$  in Eq.~\eqref{eq:hlspins}, as
\begin{equation}\label{eq:spinop_N}
\hat{S}^\alpha_i = \underbrace{\id \otimes \id \cdots \id \otimes \id}_{i-1~{\rm times}} \otimes~\hat{S}^\alpha \otimes \underbrace{\id \otimes \id \cdots \id \otimes \id}_{N-i~{\rm times}},
\end{equation}
where $\id$ is a $(2S+1)\times(2S+1)$ identity matrix.

\subsection{Truncated Holstein-Primakoff transformation}\label{sec:hptruncated}
The HP transformation shown in Eq.~\eqref{eq:hp} expresses localized spin operators in terms of bosonic operators. However, to make use of MBPT, one typically expands the square root of Eq.~\eqref{eq:hp} in a power series in $x = \hat{n}_i/2S$
\begin{equation}\label{eq:tay}
(1 - x)^{1/2} = \sum_{n=0}^{\infty} \frac{2n!}{4^n(1-2n)n!^2}x^n \approx \sum_{n=0}^{N_T} \frac{2n!}{4^n(1-2n)n!^2}x^n,
\end{equation}
which is further truncated~\cite{Chudnovsky2006,Tupitsyn2008,Yuan2020,Takei2019} to a finite number of terms $N_T$. Inserting this result in Eq.~\eqref{eq:hp}, and using thus obtained $\hat{S}^\alpha_i$ in Eq.~\eqref{eq:spin_ham}, we can re-write
\begin{equation}\label{eq:split_ham}
\hat{H}_\mathrm{lspins} = \hat{H}_{0} + \hat{H}_\mathrm{int},
\end{equation}
as the sum of two terms. Here
\begin{equation}\label{eq:ni_ham}
\hat{H}_0 = -J_\mathrm{H}(N-1)S^2 + J_\mathrm{H}S \sum_{i=1}^{N} \hat{n}_i - J_\mathrm{H}S \sum_{\braket{ij}} (\hat{a}_i^\dagger a_j + a_ia_j^\dagger),
\end{equation}
represents a one-particle Hamiltonian of noninteracting HP bosons, whereas
\begin{eqnarray}\label{eq:mag_int_ham}
\hat{H}_\mathrm{int} = -J_\mathrm{H}\sum_{\braket{ij}} &\bigg[&
  \hat{n}_i\hat{n}_j
- \frac{\hat{a}_i\hat{a}_j^\dagger \hat{n}_j}{4S}
- \frac{a_i^\dagger \hat{n}_j \hat{a}_j}{4S}
- \frac{\hat{n}_i\hat{a}_i \hat{a}_j^\dagger}{4S} \nonumber \\
&-& \frac{\hat{a}_i^\dagger\hat{n}_i\hat{a}_j}{4S}
+ \frac{\hat{n}_i\hat{a}_i\hat{a}_j^\dagger \hat{n}_j}{16S^2}
+ \frac{\hat{a}_i^\dagger \hat{n}_i\hat{n}_j\hat{a}_j}{16S^2}
\bigg] + \cdots \nonumber \\
\end{eqnarray}
is composed of many-body interaction terms that we write explicitly for truncation number $N_T=1$ to emphasize the nontrivial multi-magnon interactions present even in this lowest order of truncated HP transformation. Using the bosonic operator $\hat{a}^\dagger$ represented by an infinite matrix
\begin{equation}\label{eq:adag_mat}
\hat{a}^\dagger = \begin{pmatrix}
0 & 0& 0  & \cdots &0& \cdots \\
\sqrt{1} & 0  & 0 & \cdots &0 &\cdots \\
0 & \sqrt{2}  & 0 & \cdots &0& \cdots \\
0 & 0 & \sqrt{3}  & \ddots& \vdots & \cdots \\
0 & 0 & 0 & \cdots& \sqrt{n} & \cdots \\
\vdots & \vdots & \vdots & \vdots& \vdots & \ddots
\end{pmatrix},
\end{equation}
for a single site, we construct its matrix representation
\begin{equation}\label{eq:adag_N}
\hat{a}^\dagger_i = \underbrace{\id \otimes \id \cdots \id \otimes \id}_{i-1~{\rm times}} \otimes~\hat{a}^\dagger \otimes \underbrace{\id \otimes \id \cdots \id \otimes \id}_{N-i~{\rm times}},
\end{equation}
for the case of $N$ sites where $\id$ stands for an identity matrix of size same as of $\hat{a}^\dagger$. The matrix representation of operator $\hat{a}_i$ is the Hermitian conjugate of $\hat{a}^\dagger_i$.

\subsection{Resummed Holstein-Primakoff transformation}\label{sec:hpresum}
In numerical calculations,  $\hat{a}^\dagger$ or $\hat{a}$ are first truncated to a finite size of $N_B \times N_B$, so that the matrix representation of localized spin operators is then composed of matrix blocks associated with physical states $\left\{ \ket{n} \right\}_{n=0, \cdots, 2S}$ and unphysical states $\left\{ \ket{n}\right\}_{n=2S+1, \cdots, N_B}$, such that
\begin{equation}\label{eq:phys_states}
    \hat{S}^{\pm, z} = \begin{pmatrix}
    \hat{S}^{\pm, z}_\mathrm{phys} & \Delta_\mathrm{c} \\
    \Delta^\dagger_\mathrm{c} & \hat{S}^{\pm, z}_\mathrm{unphys}
    \end{pmatrix},
\end{equation}
where $\Delta_\mathrm{c}$ is the coupling between physical and unphysical states. The numerically exact computation of the square root of an operator in Eq.~\eqref{eq:hp} ensures $\Delta_\mathrm{c} = 0$, but Taylor expansion of square root in Eq.~\eqref{eq:tay} leads to $\Delta_\mathrm{c} \neq 0$ which, therefore, couples the physical and unphysical states. This reveals the trouble with the truncated HP transformation.

Alternatively, Ref.~\cite{Vogl2020} has very recently proposed \emph{resummed HP transformation} that furnishes a polynomial expansion for the square root in Eq.~\eqref{eq:hp}
\begin{subequations}\label{eq:hp_resum}
\begin{eqnarray}
    \hat{S}^+_i
    &\approx& \sqrt{2S} \left[ \sum_{n=0}^{N_\mathrm{max}}C_n (\hat{a}^\dagger_i)^n(\hat{a}_i)^n\right] \hat{a}_i, \label{eq:hp_resum_sp} \\
    C_n &=& \frac{1}{n!}\bigg( 1 - \frac{n}{2S}\bigg)^{1/2} - \sum_{m=0}^{n - 1} \frac{C_m}{(n-m)!}. \label{eq:hp_resum_qn}
\end{eqnarray}
\end{subequations}
This expansion ensures that when $N_\mathrm{max} = 2S$, the matrix-block $S^{\pm, z}_\mathrm{phys}$ associated with the physical states is \emph{exact}, whereas coupling between the physical and unphysical states is $\Delta_\mathrm{c} = 0$, which makes $\hat{S}^{\pm, z}_\mathrm{unphys}$ irrelevant for all practical purposes.

\subsection{Relationship between localized spin operators and their mapping to Holstein-Primakoff bosons}\label{sec:exacthp}
To interpret in physically transparent way the relationship between localized spin operators and their mapping to HP bosons, we consider 1D chain of $N=7$ sites hosting spin-$\frac{5}{2}$ localized spins. We use arrows of different length
\begin{equation}
    \bigg\downarrow, \big\downarrow, \downarrow, \uparrow, \big\uparrow, \bigg\uparrow,
\end{equation}
to denote the eigenstates  of localized spin operator $\hat{S}^z_i$ [Eq.~\eqref{eq:sz_mat}] with $m=-5/2$, $-3/2$, $-1/2$, $1/2$, $3/2$, $5/2$, respectively, as illustrated in Figs.~\ref{fig:fig1}(a) and ~\ref{fig:fig}(c). The ferromagnetic ground state of this system
\begin{equation}\label{eq:hp_vac_state}
\bigg\vert\bigg\uparrow~\bigg\uparrow~\bigg\uparrow~\bigg\uparrow~\bigg\uparrow \bigg\rangle \equiv \ket{0}
\end{equation}
is identical to HP bosonic vacuum state~$\ket{0}$ with zero HP bosons on each site, $n_i=0$, and therefore \mbox{$N_\mathrm{mag} = \sum_i n_i =0$}. Inside the ket vector on the left hand side (LHS) of  Eq.~\eqref{eq:hp_vac_state}, we indicate eigenstate with eigenvalue $m=5/2$ of the localized spin operator $\hat{S}^z_i$ for all sites $i=1$ to $i=4$. Equation~\eqref{eq:hp_vac_state} is proved by noting that $\ket{0}$ on right hand side (RHS) of Eq.~\eqref{eq:hp_vac_state} and ket on the LHS of Eq.~\eqref{eq:hp_vac_state}  are both eigenstates of the same operator $\hat{S}^z_i$ with eigenvalue $m=5/2$
\begin{subequations}\label{eq:hp_vac}
\begin{eqnarray}
    \hat{S}^z_i \bigg\vert\bigg\uparrow~\bigg\uparrow~\bigg\uparrow~\bigg\uparrow~\bigg\uparrow \bigg\rangle
    &=& \frac{5}{2} \bigg\vert\bigg\uparrow~\bigg\uparrow~\bigg\uparrow~\bigg\uparrow~\bigg\uparrow \bigg\rangle, \\
    \hat{S}^z_i\ket{0} = (5/2 - \hat{n}_i)\ket{0}
    &=&  5/2\ket{0},
\end{eqnarray}
\end{subequations}
and, therefore, they must be identical. Thus, creating $n_i=1$ or $n_i = 2$ HP bosons on site $i=1$ corresponds to reducing the size of localized spin at site $i=1$ by $1$ or $2$ units $\hbar$, i.e.,  $m=5/2 \mapsto 3/2$ or $m=5/2 \mapsto 1/2$, respectively. We depict this graphically by using
\begin{eqnarray}\label{eq:hp_one_mag}
\hat{a}_1^\dagger \ket{0}
&=& \bigg\vert\big\uparrow~\bigg\uparrow~\bigg\uparrow~\bigg\uparrow~ \bigg\uparrow \bigg\rangle, \\
(\hat{a}_1^\dagger)^2 \ket{0} &=& \bigg\vert \hspace{-0.04cm}\uparrow~ \bigg\uparrow~\bigg\uparrow~\bigg\uparrow~\bigg\uparrow \bigg\rangle.
\end{eqnarray}
In addition, for the total of $N_\mathrm{mag} = 2$ HP bosons created on different sites, such as sites $i=1$ and $i=3$,
\begin{equation}\label{eq:hp_two_mag2}
\hat{a}^\dagger_1 \hat{a}^\dagger_3 \ket{0} =  \bigg\vert\big\uparrow~ \bigg\uparrow~\big\uparrow~\bigg\uparrow~\bigg\uparrow~\bigg\rangle.
\end{equation}
Thus, creating the total of $N_\mathrm{mag}= n$ HP bosons ($n \in 0, 1, 2, \ldots$) is interpreted physically as the reduction of the size of localized spins by $n$ units $\hbar$. However, in the quantum state $(\hat{a}_i^\dagger)^n|0\rangle$, the expectation value of localized spin operator is $\braket{\hat{S}^z_i} = S - n$, so that the constraint $0 \le n \le 2S$ must be obeyed in order to remain in the subspace of physical states shown in Eq.~\eqref{eq:phys_states}. In other words, on a given site $i$ one {\em cannot} create $N_\mathrm{mag} > 2S$ HP bosons.

\subsection{Numerically exact time evolution of quantum many-body states}\label{sec:cn}

Using Eq.~\eqref{eq:ele_create} and Eq.~\eqref{eq:spinop_N} in Eq.~\eqref{eq:ham_sys} yields the matrix representation of many-body Hamiltonian $\hat{H}$. The solution of time-dependent Schr\"{o}dinger equation for quantum many-body state
\begin{equation}\label{eq:psievolve}
	i\hbar \frac{\ket{\Psi(t)}}{dt} = \hat{H}(t) \ket{\Psi(t)},
\end{equation}
is formally given by
\begin{equation}\label{eq:evolutionop}
	\ket{\Psi(t+\delta t)} = \mathcal{T} \exp \left( -\frac{i}{\hbar}\int\limits_t^{t+\delta t} \hat{H}(t) \right) \ket{\Psi(t+\delta t)},
\end{equation}
where $\mathcal{T}$ is the time-ordering operator. While many numerical algorithms are available to propagate Eq.~\eqref{eq:evolutionop}, including direct computation of matrix exponential when $\hat{H}$ is time independent, in general by using small $\delta t$ and by considering $\hat{H}(t)$ constant over such time step the Crank-Nicolson algorithm
\begin{equation}\label{eq:cn}
	\left( 1 + \frac{i\delta t}{2\hbar}\hat{H}(t) \right) \ket{\Psi(t+\delta t)} = \left( 1 - \frac{i\delta t}{2\hbar}\hat{H}(t) \right) \ket{\Psi(t)},
\end{equation}
we employ offers propagation scheme which is unitary, accurate to second-order in $\delta t$ and unconditionally stable~\cite{Wells2019}.

Using thus obtained $\ket{\Psi(t)}$, the time-dependent expectation value of a localized spin operator at site $i$ is given by
\begin{equation}\label{eq:expectationspin}
	S^\alpha_i(t) = \bra{\Psi(t)}\hat{S}^\alpha_i\ket{\Psi(t)}.
\end{equation}
When the localized spin operators are represented directly by finite size matrices in Eqs.~\eqref{eq:spin_mat} and ~\eqref{eq:spinop_N}, the corresponding expectation values $S^\alpha_i(t)$ are ~\emph{numerically~exact} and, therefore, serve as a benchmark for alternative computation of the same expectation value when $\hat{S}^\alpha_i$ are represented by polynomial expressions in bosonic operators introduced in Secs.~\ref{sec:hptruncated} and ~\ref{sec:hpresum}.

\subsection{From Holstein-Primakoff bosons to one- or multi-magnon Fock states}\label{sec:magnons}

In contrast to HP bosons created on a given site, $\hat{a}_i^\dagger|0\rangle$, which are not the eigenstates of $\hat{H}_0$ in Eq.~\eqref{eq:ni_ham}, one-magnon states are linear combinations of $\hat{a}_i^\dagger|0\rangle$ which diagonalize Hamiltonian $\hat{H}_0$
\begin{equation}\label{eq:eigenmagnon}
	\hat{H}_0 |q\rangle = [E_0 + \hbar \omega(q)]|q\rangle.
\end{equation}
Thus, they can be visualized as a completely ``delocalized'' over all sites bosonic quasiparticle which carries momentum $\hbar q$ (in 1D we use for examples) and angular momentum $1\hbar$. Here $E_0=-2J_\mathrm{H}S^2 N$ is the ground state energy of ferromagnetic spin chain.

To find explicit expression for $|q\rangle$, we consider 1D chain [Fig.~\ref{fig:fig1}(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_\mathrm{H}$ in $\hat{H}_0$ in Eq.~\eqref{eq:ni_ham}. For the clarity  of notation, we use $\uparrow$, $\odot$, and $\downarrow$  to denote eigenstates of localized spin operator $\hat{S}^z_i$ with eigenvalues [Eq.~\eqref{eq:sz_mat}] $m=1$, $m=0$, and $m=-1$, respectively. The one-magnon state is then given by
\begin{equation}\label{eq:locs_one_mag}
    \ket{q} \equiv \frac{1}{\sqrt{N}}
    \sum_{n=0}^{N-1}e^{i q x_n}\ket{ \underbrace{\uparrow \cdots\uparrow}_{n~\text{times}} \odot  \underbrace{\uparrow\cdots \uparrow}_{N-n-1~\text{times}}},
\end{equation}
where $x_n = na_0$ is the real-space position of the localized spin at site $n+1$ and $a_0$ is the lattice spacing. This then determines magnon energy-momentum dispersion in Eq.~\eqref{eq:eigenmagnon},  $\hbar \omega(q) = 2JS[1-\cos(qa_0)]$. The expectation value of the total $z$-spin operator of localized spins is given by
\begin{equation}\label{eq:sam_reduction}
   \bra{q}\sum_{i=1}^N \hat{S}^z_i\ket{q} = (NS - 1),
\end{equation}
which indicates that a unit of total $z$-spin is removed from the ferromagnetic ground state. Thus, both one HP boson state and one HP magnon state are labeled by  $N_\mathrm{mag} = 1$ throughout the paper. In addition, the expectation value of the localized $z$-spin operator at arbitrary site $i$
\begin{equation}\label{eq:loc_sam_reduction}
    \bra{q} \hat{S}^z_i \ket{q} = \left( S - \frac{1}{N}\right),
\end{equation}
shows that excitation of a magnon reduces the $z$-component of each spin by $1/N$. This also underlies the LLG picture of spin wave comprised of precessing classical vectors with frequency $\omega$ and some small cone angle, where the phase of the precession of adjacent vectors varies harmonically~\cite{Kim2010}.

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{fig3.pdf}
	\caption{ (a) Spatio-temporal profile of $S^z_i(t)$ within 1D chain in Fig.~\ref{fig:fig1}(a) composed of $N=7$ sites hosting spin-$1$ localized spins where at $t=0$ $N_\mathrm{mag} = 1$ HP boson is introduced at site $i=1$, $\ket{\Psi(t=0)} = \ket{\odot \uparrow \uparrow\uparrow \uparrow \uparrow \uparrow}$. (b) The corresponding overlap, $n_q(t) = |\braket{q|\Psi(t)}|^2$, with {\em one-magnon Fock state} in $q > 0$ region as the HP boson moves from the left to the right edge of the chain [white shadow in (a)] which switches towards $q < 0$ region when reflection occurs near the boundary at site $i=7$ and $t \gtrsim 2.5$ fs (indicated by dashed horizontal lines). (c) Same as (a) but at $t=0$, $N_\mathrm{mag} = 2$ HP bosons are introduced at site $i=1$, $\ket{\Psi(t=0)} = \ket{\downarrow \uparrow \uparrow \uparrow \uparrow \uparrow \uparrow}$. (d) The corresponding overlap, $n_{q, q'}(t) = |\braket{q, q'|\Psi(t)}|^2$, with {\em two-magnon Fock state}. Dashed horizontal lines refer to times \mbox{$t = 1$ fs} and \mbox{$t = 4$ fs}. In all panels $J_H = 1$ eV.}
	\label{fig:fig3}
\end{figure}

In the second quantized description produced by HP transformation, $\ket{q} = \hat{b}_q^\dagger\ket{0}$ is {\em one-magnon Fock state}~\cite{Quirion2020} where the creation operator of HP magnon is given by
\begin{eqnarray}\label{eq:fockmagnon}
        \hat{b}^\dagger_q = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} e^{ikx_n}\hat{a}^\dagger_{n+1}.
\end{eqnarray}
The one-magnon Fock state has been realized experimentally~\cite{Quirion2020} only very recently in a millimeter-sized ferrimagnetic crystal and detected by superconducting qubit as quantum sensor, thereby representing  a single-photon detection counterpart in quantum magnonics.

Note that in spintronics and magnonics literature~\cite{Kamra2016} one sometimes finds $\hat{a}_i^\dagger|0\rangle$ denoted as ``one magnon created in real space at position $i$'' while $\hat{b}^\dagger_q |0\rangle$ is ``one magnon created in the reciprocal space with momentum q''. However, the former is not an eigenstate of Hamiltonian in Eq.~\eqref{eq:eigenmagnon}, while the later is, so we differentiate between them by using ``HP boson'' for the former and  ``HP magnon'' for the latter. Nevertheless, for both situations we use label $N_\mathrm{mag} =1$ for simplicity of notation because in both cases a unit of total $z$-spin is removed from the ferromagnetic ground state, as shown in Eq.~\eqref{eq:sam_reduction}.

The Fock states of $N_\mathrm{mag}=2$ magnons carrying momentum $\hbar q$ and $\hbar q'$ are defined by~\cite{Morimae2005}
\begin{equation}\label{eq:two_mag_bs}
	\ket{q, q'} = \sum_{n>m} f_{mn}(q, q') \ket{\underbrace{\uparrow\cdots\uparrow}_{m-1~\text{times}} \odot \uparrow\cdots\uparrow \odot \underbrace{\uparrow \cdots \uparrow}_{N-n~\text{times}}},
\end{equation}
where
\begin{equation}
	f_{mn}(q, q') = \frac{1}{\mathcal{N}} \bigg[ e^{i q x_m}e^{iq'x_n} + e^{iq'x_n}e^{iqx_m}\bigg].
\end{equation}
Here $\mathcal{N}$ is normalization constant~\cite{Morimae2005}, and $f_{mn}(q, q')$ is symmetric under exchange $q \leftrightarrow q'$ leading to $\ket{q, q'} = \ket{q', q}$ which  reflects the bosonic nature of localized spin operators that commute on different sites $i \neq j$. The is manifestly encoded in the second quantization formalism since $\ket{q, q'} = \hat{b}^\dagger_q\hat{b}^\dagger_{q'}\ket{0} = \hat{b}^\dagger_{q'}\hat{b}^\dagger_{q}\ket{0} = \ket{q', q}$.

To illustrate the distinction between HP boson and HP magnon, we initialize $N=7$ site chain [Fig.~\ref{fig:fig1}(a)] in quantum state $\ket{\Psi_1(t=0)} = \ket{\odot \uparrow \uparrow \uparrow \uparrow \uparrow \uparrow}$ in Figs.~\ref{fig:fig3}(a) and ~\ref{fig:fig3}(b); or in quantum state $\ket{\Psi_2(t=0)}=\ket{\downarrow \uparrow \uparrow \uparrow \uparrow \uparrow \uparrow}$ in Figs.~\ref{fig:fig3}(c) and ~\ref{fig:fig3}(d). This means that $N_\mathrm{mag}=1$ HP boson is created at site $i=1$ at $t=0$ in the former case, while full flip of localized spin at site $i=1$ in the latter case means that $N_\mathrm{mag} = 2$ HP bosons are created at site $i=1$. Besides pedagogical value, such initial states and one or two magnon propagation including magnon bound states, has also been observed 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~\cite{Fukuhara2013}.

Since $\ket{\Psi_1(t=0)}$ is not an eigenstate, it evolves in time to produce spatio-temporal profile of $S^z_i(t)$ in $\ket{\Psi_1(t)}$ plotted in Fig.~\ref{fig:fig3}(a). For quantum time evolution we use the scheme explained in Sec.~\ref{sec:cn} where interacting Hamiltonian $\hat{H}_\mathrm{lspins}$ from Eq.~\eqref{eq:hlspins} is plugged in, but since only one HP boson is excited this is equivalent to using noninteracting $\hat{H}_0$ in Eq.~\eqref{eq:ni_ham}. The white trace  Fig.~\ref{fig:fig3}(a) indicates how HP boson moves from the left to the right edge of the chain while undergoing reflection  at site $i=7$ at \mbox{$t \approx 2.5 $ fs}, as indicated by horizontal dashed line, followed by multiple back-and-forth reflections. Note that since 1D chain in Fig.~\ref{fig:fig1}(a) has open boundary conditions, its low-energy excited eigenstates differ~\cite{Haque2010} from textbook~\cite{Mahan2011,Chudnovsky2006}  HP magnons $|q\rangle$ in Eq.~\eqref{eq:locs_one_mag} as eigenstates of interacting localized spin systems with translational invariance (that is, either infinite or finite length but with periodic boundary conditions). Nevertheless, it is insightful to visualize in Fig.~\ref{fig:fig3}(b) the  overlap, $n_q(t) = |\braket{q|\Psi(t)}|^2$, between many-body quantum state  $\ket{\Psi_1(t)}$  with one HP boson and one-magnon Fock state $|q\rangle$. Large values of $n_q(t)$ are observed in the region where $q>0$ and \mbox{$t \lesssim 2.5$ fs}, coinciding with left-to-right motion of HP boson in Fig.~\ref{fig:fig3}(a), which signifies excitation of $N_\mathrm{mag}=1$ magnon with positive momentum. On the other hand, after reflection of the HP boson near the boundary at site $i=7$ and \mbox{$t \approx  2.5$ fs}, a rapid rise of $n_{q}(t)$ in the $q < 0$ region is observed which indicates excitation of $N_\mathrm{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.

Figure~\ref{fig:fig3}(c) plots spatio-temporal profile of $S^z_i(t)$ in state $\ket{\Psi_2(t)}$ starting from $\ket{\Psi_2(t=0)} = \ket{\downarrow \uparrow \uparrow \uparrow \uparrow \uparrow \uparrow}$. For quantum time evolution we use the scheme explained in Sec.~\ref{sec:cn} where interacting Hamiltonian $\hat{H}_\mathrm{lspins}$ from Eq.~\eqref{eq:hlspins} is plugged in, so that two HP bosons are correlated both by the Pauli exclusion principle and by interactions in $\hat{H}_\mathrm{int}$ [Eq.~\ref{eq:mag_int_ham}] where $N_T \rightarrow \infty$. The two HP bosons move immediately from the left to the right for $t>0$, as shown by white traces in Fig.~\ref{fig:fig3}(c). The corresponding overlap, $n_{q, q'}(t) = |\braket{q, q'|\Psi(t)}|^2$, in Fig.~\ref{fig:fig3}(d) is $n_{q, q'}(t) = 0$ for $t \lesssim 1$ fs which is explained by Eq.~\eqref{eq:two_mag_bs} where two-magnon Fock state is composed  of terms containing two HP bosons on different sites $m\neq n$. Since at $t=0$ the two HP bosons are on the same site $i=1$, we find $n_{q, q'}(t)=0$. However, this holds until $t \lesssim 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\lesssim t \lesssim 4$ fs (indicated by horizontal dahsed line) the region near \mbox{$q = 0.5$ $\pi/a_0$} shows large values of $n_{q, q'}(t)$, and since \mbox{$q'=0.5$ $\pi/a_0$} is fixed for all values of $q$ in Fig.~\ref{fig:fig3}(d), we can conclude that two HP bosons posses nearly the same velocity. Beyond $t~\approx~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 to the right while the other moves to the left edge of the chain.

\subsection{Retarded and lesser one-particle Green functions}\label{ref:gf_sf}
The fundamental quantities of nonequilibirum GF formalism~\cite{Stefanucci2013,Schlunzen2020} for fermions are the one-particle retarded GF

\begin{equation}\label{eq:ele_GF_ret}
G^r_{i\sigma, j\sigma'}(t,t') =
-i\hbar^{-1}\Theta(t-t')\braket{\{\hat{c}_{i\sigma}(t),\hat{c}_{j\sigma'}^\dagger(t') \}},
\end{equation}
and the one-particle lesser GF
\begin{equation}\label{eq:ele_GF_less}
G^<_{i\sigma,j\sigma'}(t,t') =
i\hbar^{-1} \braket{\hat{c}_{j\sigma'}^\dagger(t')\hat{c}_{i\sigma}(t)},
\end{equation}
which describe the density of available quantum states and how electrons occupy those states, respectively. Here $\Theta(t-t')$ is the Heaviside-function; $\hat{c}_{i\sigma}(t)$ indicates Heisenberg time-evolution of $\hat{c}_{i\sigma}$; and $\braket{\cdots} = \mathrm{Tr}(\hat{\rho} \cdots)$ is the quantum statistical average, where $\hat{\rho}$ is the density operator of the system at $t=0$.

Analogously, the bosonic one-particle retarded GF is given by
\begin{equation}\label{eq:mag_GF_ret}
D^r_{ij}(t,t') =
-i\hbar^{-1}\Theta(t-t')\braket{[\hat{a}_{i}(t),\hat{a}_{j}^\dagger(t')]},
\end{equation}
and the lesser one is given by
\begin{equation}\label{eq:mag_GF_less}
D^<_{ij}(t,t') =
-i \hbar^{-1}\braket{\hat{a}_{j}^\dagger(t')\hat{a}_{i}(t)}.
\end{equation}
At zero-temperature and in equilibrium, where \mbox{$\hat{\rho} = \ket{\Psi_0}\bra{\Psi_0}$} with $\ket{\Psi_0}$ being the ground state of the many-body Hamiltonian $\hat{H}$ in Eq.~\eqref{eq:ham_sys}, these GFs  depend solely on $t-t'$ and can be Fourier transformed to energy. Such energy-dependent electronic retarded and lesser GFs are given by
\begin{equation}\label{eq:ele_FTGF}
G^{r,<}_{i\sigma,j\sigma'}(E) = \frac{1}{2\pi\hbar}
\int\displaylimits_{-\infty}^{+\infty} dE\,
G^{r,<}_{i\sigma,j\sigma'}(t-t') e^{iE(t-t')/\hbar},
\end{equation}
whereas the bosonic ones are given by
\begin{equation}\label{eq:mag_FTGF}
D^{r,<}_{ij}(E) =  \frac{1}{2\pi\hbar}
\int\displaylimits_{-\infty}^{+\infty} dE\,
D^{r,<}_{ij}(t-t') e^{iE(t-t')/\hbar}.
\end{equation}

\subsection{Spectral function for electrons and magnons}
The  electronic spectral function $A(E)$, or the ``interacting density of states''~\cite{Balzer2011,Nocera2018}, is computed using the retarded GF in Eq.~\eqref{eq:ele_FTGF} as
\begin{eqnarray}\label{eq:ele_spec_f}
A(E) &=&
-2\sum_{i=1}^{N}\sum_{\sigma=\uparrow, \downarrow} \mathrm{Img}[G^r_{i\sigma,i\sigma}(E)] \nonumber \\
&=& \sum_{k} W_k^+ \delta(E-\Delta_k) + W_k^-\delta(E+\Delta_k),
\end{eqnarray}
where $\Delta_k = (E_k - E_0)$ and $E_k$ are the eigenenergies of many-body Hamiltonian $\hat{H}$ in Eq.~\eqref{eq:ham_sys}. The pre-factor of $\delta$-function peaks in $A(E)$
\begin{subequations}\label{eq:ele_weight}
\begin{eqnarray}
W_k^+ &=& \sum_{i=1}^N\sum_{\sigma=\uparrow, \downarrow}|\bra{\Psi_k} \hat{c}_{i\sigma}^\dagger\ket{\Psi_0}|^2, \\
W_k^- &=& \sum_{i=1}^N\sum_{\sigma=\uparrow, \downarrow}|\bra{\Psi_k} \hat{c}_{i\sigma}\ket{\Psi_0}|^2,
\end{eqnarray}
\end{subequations}
define the ``weight'' of the many-body eigenstate $\ket{\Psi_k}$ within $A(E)$. Since the ground state $\ket{\Psi_0}$ is an eigenstate of the electron number operator $\hat{N}_\mathrm{e}$ [Eq.~\eqref{eq:ele_num}], it has a well-defined number  of electrons $N_\mathrm{e}$. Thus, the action of $\hat{c}_{i\sigma}^\dagger$ and $\hat{c}_{i\sigma}$ on $\ket{\Psi_0}$ in Eq.~\eqref{eq:ele_weight} reveals that the $\delta$-function peaks at $E = \pm \Delta_k$ can only come from many-body eigenstates $\ket{\Psi_k}$ that posses  $N_\mathrm{e}\pm 1$ electrons.

Similarly, the bosonic spectral function $D(E)$ is evaluated using the bosonic retarded GF  in Eq.~\eqref{eq:mag_FTGF}
\begin{eqnarray}\label{eq:mag_spec_f}
D(E) &=& -2\sum_{i=1}^{N}\mathrm{Img}[D^r_{ii}(E)] \nonumber \\
&=& \sum_{k} Q_k^+ \delta(E-\Delta_k) - Q_k^-\delta(E+\Delta_k),
\end{eqnarray}
where
\begin{subequations}\label{eq:mag_weight}
\begin{eqnarray}
Q_k^+ &=& \sum_{i}|\bra{\Psi_k} \hat{a}_{i}^\dagger\ket{\Psi_0}|^2,
\\
Q_k^- &=& \sum_{i}|\bra{\Psi_k} \hat{a}_{i}\ket{\Psi_0}|^2.
\end{eqnarray}
\end{subequations}
define the ``weight'' of many-body eigenstate $\ket{\Psi_k}$ within $D(E)$. The $\delta$-function peaks in Eq.~\eqref{eq:mag_spec_f} at $E = \pm \Delta_k$ come from many-body eigenstates $\ket{\Psi_k}$. However, unlike the electronic case, they do not have a  well-defined total magnon number $N_\mathrm{mag}$ as they are not eigenstates of the total magnon number operator $\hat{N}_\mathrm{mag} = \sum_{i=1}^N \hat{n}_i$. This can be easily seen in Fig.~\ref{fig:fig9}(e) where the structure of one such many-body eigenstate $\ket{\Psi_k}$ is shown to be a linear combination of many-body states with total magnon number $N_\mathrm{mag} = 0$, $N_\mathrm{mag} = 1$ and $N_\mathrm{mag} = 2$.

Both $A(E)$ and $D(E)$ satisfy the sum rule
\begin{subequations}
    \begin{eqnarray}
        \int\displaylimits_{-\infty}^{+\infty} \frac{dE}{2\pi} A(E)   &=& 2N, \label{eq:ele_sum_rule}\\
        \int\displaylimits_{-\infty}^{+\infty} \frac{dE}{2\pi} D(E)  &=& N.
        \label{eq:mag_sum_rule}
    \end{eqnarray}
\end{subequations}
The sum rule allows us to invoke 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 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.~\ref{fig:fig1}(b) and described by Hamiltonian in Eq.~\eqref{eq:ham_sys}, includes HP bosons interacting [Eq~\eqref{eq:mag_int_ham}] with other HP bosons when $N_\mathrm{mag} > 1$ and electrons interacting with HP bosons while electron-electron interactions are excluded. Since the sum rule is an exact result, in actual GF calculations it is often employed to test the quality of the approximations made~\cite{Mahfouzi2014}.

\section{Results and Discussion}\label{sec:randd}
\subsection{Range of validity of truncated HP transformation for nonequilibrium magnon-magnon interacting systems}\label{sec:validity_mm}

Figure~\ref{fig:fig4} compares $S_1^z(t)$ for 1D chain [Fig.~\ref{fig:fig1}(a)], hosting $S=5/2$ or $S=1$ localized spins and without any electrons, computed using original localized spin operators vs. their mapping to bosonic operators via  truncated HP transformation. In the ferromagnetic ground state $\ket{\Psi_0}$, the expectation value $S^z_i = \braket{\Psi_0|\hat{S}^z_i|\Psi_0} = 5/2$ for all sites $i$ at $t=0$ in Fig.~\ref{fig:fig4}(a)--(c). To initiate nonequilibrium dynamics for $t>0$, we choose an initial state $\ket{\Psi(0)}$ such that the expectation value of the localized spin at site $i=1$ is reduced by $N_\mathrm{mag}$ units, i.e., $S^z_{i=1}(0) = \braket{\Psi(0)|\hat{S}^z_{i=1}|\Psi(0)} = (5/2 - N_\mathrm{mag})$ while  $S^z_{i\neq 1}(0) = 5/2$.  From the viewpoint of HP transformation, this is equivalent to introducing $N_\mathrm{mag}$ HP bosons at site $i=1$ at $t=0$ and the initial state is given by
\begin{equation}
	\ket{\Psi(0)} = (\hat{a}_1^\dagger)^{N_\mathrm{mag}} \ket{0}.
\end{equation}
When $N_\mathrm{mag} = 1$, Fig.~\ref{fig:fig4}(a) shows that $S^z_1(t)$, evaluated by truncated HP transformation (green dashed line) solely containing single-particle Hamiltonian $\hat{H}_0$ of noninteracting HP bosons in Eq.~\eqref{eq:ni_ham}, accurately matches the exact time dependence  (black lines) of $S^z_1(t)$ evaluated using the localized spin operators. This result is intuitively expected---since there is only \emph{one} HP boson in the system,  magnon-magnon interaction terms active within $\hat{H}_\mathrm{int}$ part of the Hamiltonian [Eq.~\eqref{eq:mag_int_ham}] cannot influence the dynamics of localized spins.

To understand the significance of magnon-magnon interaction terms within $\hat{H}_\mathrm{int}$ on the dynamics of localized spins, in the next step we  introduce $N_\mathrm{mag} = 2$ HP bosons at site $i=1$. The time dependence of $S^z_1(t)$ in Fig.~\ref{fig:fig4}(b), evaluated via the truncated HP transformation with a truncation number $N_T = 1$ (Sec.~\ref{sec:hptruncated}), matches the exact one only for short enough times ($0 < t < 10$~fs). At longer times ($30 <t < 35$~fs), disagreement emerges due to missing effects from $N_T>1$ magnon-magnon interaction terms within $\hat{H}_\mathrm{int}$. Thus, to recover the agreement between two types of calculations at longer times requires increasing $N_T$, such as $N_T=3$ (orange dashed lines) in Fig.~\ref{fig:fig4}(b).

However, progressively larger $N_T$ must be used [Fig.~\ref{fig:fig5}] to extend breakdown time $t_\mathrm{break}$, marked in Fig.~\ref{fig:fig4}(d),(e), at which disagreement between two types of calculations emerges. We define $t_\mathrm{break}$ as the time when the deviation $\Delta~=~S^z_1(t)|_\mathrm{HP} - S^z_1(t)|_\mathrm{exact}$, between $S^z_1(t)|_\mathrm{HP}$ (evaluated from truncated HP transformation  with a truncation number $N_T$) and the exact $S^z_1(t)|_\mathrm{exact}$ becomes larger than the chosen tolerance, $|\Delta| > 10^{-4}$. As demonstrated by Figs.~\ref{fig:fig4}(d),(e) and Fig.~\ref{fig:fig5}, $t_\mathrm{break}$ sensitively depends on the density of HP bosons, $N_\mathrm{mag}/N$, whose increase makes magnon-magnon interaction terms within $\hat{H}_\mathrm{int}$ more relevant, thereby requiring larger $N_T$ in Fig.~\ref{fig:fig5}. Figure~\ref{fig:fig4}(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 longer quantum spin chains [than $N=4$ sites used in Fig.~\ref{fig:fig4}(a)--(c) with spin-$5/2$ on each site]. In Fig.~\ref{fig:fig4}(d), at $t=0$ we flip the localized spins at sites $i=1$ and $i=2$ (see inset), thereby introducing two HP bosons on each of these sites. Thus, the total number of HP bosons within the system is $N_\mathrm{mag} = 4$, whereas the HP boson density is $N_\mathrm{mag}/N = 1$. For such parameters, \mbox{$t_\mathrm{break} \approx 11 $ fs} ($|\Delta| = 0.05$ is chosen solely for visualization of $t_\mathrm{break}$ at fs time scales).  On the other hand, in Fig.~\ref{fig:fig4}(e), where HP boson density is reduced to $N_\mathrm{mag}/N = 0.57$ by making 1D chain longer from $N= 4$ sites to $N=7$ sites, we find that $t_\mathrm{break}$ for truncated HP transformation increases to \mbox{$t_\mathrm{break} \approx 23$ fs}.

\begin{figure}
	\centering
	\includegraphics[width=8cm]{fig4.pdf}
	\caption{Comparison of exact trajectory of $S^z_1(t)$ (black lines) evaluated using localized spin operators vs. the one obtained by truncated Holstein-Primakoff (HP) transformation power-series [Eq.\eqref{eq:tay}] with a truncation number $N_T$, for 1D quantum spin chain illustrated in Fig.~\ref{fig:fig1}(a) composed of $N=4$ sites hosting spin-$\frac{5}{2}$ [panels (a)--(c)] or spin-1 [panels (d)--(e)] localized spins. At $t=0$, $N_\mathrm{mag}$ HP bosons are introduced at site $i=1$ or $2$, indicated by the reduced size of or flipped localized spin in the inset. In panel~(a), $N_\mathrm{mag} = 1$ and the green dotted line evaluated \emph{solely} from noninteracting HP boson Hamiltonian [Eq.~\eqref{eq:ni_ham}] is identical to the exact trajectory of $S_1^z(t)$. In panel~(b), $N_\mathrm{mag}=2$ and $S_1^z(t)$ evaluated from truncated HP transformation power-series with $N_T=1$ disagrees, while that for $N_T=3$ matches the exact one. In panel~(c), more HP magnons are introduced ($N_\mathrm{mag} = 10$) by \emph{flipping two localized spins} while retaining $N_T=3$ to show emergence of disagreement with exact. In panels (d) and (e), $N=4$ and $7$ (see inset), respectively, where orange solid lines indicate $S_1^z(t)$ evaluated from truncated HP power-series with $N_T=2$. The vertical dot-dash lines indicate $t=t_\mathrm{break}$ at which truncated HP transformation power-series breaks down (see Fig.~\ref{fig:fig4}). Here \mbox{$J_\mathrm{H} = 1$} eV.}
	\label{fig:fig4}
\end{figure}
\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{fig5.pdf}
	\caption{Breakdown time for truncated Holstein-Primakoff (HP) transformation power-series [Eq.~\eqref{eq:tay}] with a truncation number $N_T$ for 1D quantum spin chains illustrated in Fig.~\ref{fig:fig1}(a) composed of $N=4$ sites hosting \mbox{spin-$\frac{5}{2}$} localized spins. The inset near each line indicates the number of HP bosons $N_\mathrm{mag}$, introduced in the chain at $t=0$ depicted by the reduced size of or flipped localized spin at site $i=1$. Here \mbox{$J_\mathrm{H} = 1$ eV}.}
	\label{fig:fig5}
\end{figure}

Figure~\ref{fig:fig5} demonstrates how for a given time $t=t_\mathrm{break}$, the horizontal distance between consecutive curves from left to right increases nonlinearly. This means that $N_T$ needed to accurately track $S_1^z(t)$ via the truncated HP transformation increases nonlinearly with the number of HP bosons $N_\mathrm{mag}$ present in the system. Since we use logarithmic scale for the ordinate axis of Fig.~\ref{fig:fig5}, a roughly constant slope of all curves suggests \mbox{$t_\mathrm{break} \propto e^{N_T}$}. This indicates that with increasing time more and more magnon-magnon interaction terms become relevant within $\Hat{H}_\mathrm{int}$, so they have to be included by increasing $N_T$. While small $N_T$ greatly simplifies MBPT for electron-boson interacting systems~\cite{Marini2018}, using \mbox{$N_T=1$--$5$} allows one to track dynamics of localized spins only up to time \mbox{$t_\mathrm{break} \approx 1.5\hbar/ J_\mathrm{H} \approx 15.0$ fs} (for \mbox{$J_\mathrm{H} = 0.1$ eV}). This is insufficient to model even ultrafast optical manipulation of magnetism requiring simulation times of few tens of fs~\cite{Siegrist2019}, and it is much further away  from current-driven magnetization dynamics via spin torque which occurs on $\sim 1$ ns time scales~\cite{Ralph2008,Berkov2008}.

\begin{figure*}
	\centering
	\includegraphics[width=\linewidth]{fig6.pdf}
	\caption{(a) Comparison of the exact (black line) time dependence $S_1^z(t)$ evaluated using localized spin operators vs. the one obtained by truncated HP transformation [Eq.~\eqref{eq:tay}] with a truncation number $N_T = 3$ (orange line) for 1D quantum many-body system  illustrated in Fig.~\ref{fig:fig1}(b) comprised of $N=3$ sites hosting spin-$\frac{5}{2}$ localized spins. The electron--localized-spin $sd$ exchange interaction is turned off (\mbox{$J_\mathrm{sd} = 0$ eV}) in (a) for reference, and turned on in (b) using \mbox{$J_\mathrm{sd} = 1$ eV}. The insets in both panels depict the state of localized spins at $t=0$, where the localized spin at site $i=1$ is completely flipped [i.e., $N_\mathrm{mag} = 5$ HP bosons are introduced at site $i=1$ via Eq.~\eqref{eq:fivehpbosons}] to initiate nonequilibrium dynamics.  With the $sd$ interaction turned on,  the disagreement between the 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 interacting magnons, $\sum_{i=1}^N \bra{\Psi(t)} \hat{n}_i \ket{\Psi(t)}$, from the localized spin subsystem to the electronic subsystem. Panels (d) and (e) are counterparts of panels (a) and (b), respectively, for $S_1^z(t)$ evaluated using the resummed HP transformation in Eq.~\eqref{eq:hp_resum} which is approximate when $N_\mathrm{max} = 3$ but \emph{exact} for $N_\mathrm{max} = 2S = 5$.
	}
	\label{fig:fig6}
\end{figure*}

Figure~\ref{fig:fig4}(d),(e) shows explicitly 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 longer quantum spin chains [than $N=4$ sites used in Fig.~\ref{fig:fig4}(a)--(c) with spin-$5/2$ on each site]. In Fig.~\ref{fig:fig4}(d), at $t=0$ we flip the localized spins at sites $i=1$ and $i=2$ (see inset), thereby introducing two HP bosons on each of these sites. Thus, the total number of HP bosons within the system is $N_\mathrm{mag} = 4$, whereas the HP boson density is $N_\mathrm{mag}/N = 1$. For such parameters, \mbox{$t_\mathrm{break} \approx 11 $ fs} ($|\Delta| = 0.05$ is chosen solely for visualization of $t_\mathrm{break}$ at fs time scales).  On the other hand, in Fig.~\ref{fig:fig4}(e), where HP boson density is reduced to $N_\mathrm{mag}/N = 0.57$ by making 1D chain longer from $N= 4$ sites to $N=7$ sites, we find that $t_\mathrm{break}$ for truncated HP transformation increases to \mbox{$t_\mathrm{break} \approx 23$ fs}. Qualitatively, this observation is easily explained. For a longer 1D chain, the probability for magnon-magnon scattering events is low which makes inclusion of high-order magnon-magnon interaction terms less important and thus the breakdown-time for truncated HP transformation is larger.

\subsection{Range of validity of truncated HP transformation for  nonequilibrium electron-magnon interacting systems}\label{sec:validity_em}

In this Section, we repeat the same analysis as in Sec.~\ref{sec:validity_mm} but with electron--localized-spin interaction turned on within 1D quantum many-body system illustrated in Fig.~\ref{fig:fig1}(b) composed of $N=3$ sites. These sites host spin-$\frac{5}{2}$ localized spins interacting with ``half-filled'' ($N_\mathrm{e} = 3$) tight-binding electrons via the $sd$ exchange interaction~\cite{Cooper1967} of strength $J_\mathrm{sd}$ as encoded by Eq.~\eqref{eq:sd_ham}.

At $t=0$, we completely flip the localized spin at site $i=1$ to initiate nonequilibrium dynamics. From the viewpoint of HP transformation, this is equivalent to introducing $N_\mathrm{mag} = 5$ HP bosons at site $i=1$, and thus, the initial  many-body quantum state is given by
\begin{equation}\label{eq:fivehpbosons}
    \ket{\Psi(0)} = (\hat{a}_1^\dagger)^5\hat{c}^\dagger_{2\uparrow} \hat{c}^\dagger_{1\uparrow}\hat{c}^\dagger_{1\downarrow}\ket{0}.
\end{equation}
Here $\ket{0}$ is the vacuum state of electrons and HP bosons combined. Figure~\ref{fig:fig6}(a) with electron--localized-spin interaction turned-off ($J_\mathrm{sd}$ = 0) serves as a reference. When electron--localized-spin interaction is  turned on (\mbox{$J_\mathrm{sd} = 1$ eV}) in  Fig.~\ref{fig:fig6}(b),  $S_1^z(t)$ computed exactly vs. using truncated HP transformation  follow each other for longer times \mbox{$t \approx 15$ fs} than
in Fig.~\ref{fig:fig6}(a). This is explained by  Fig.~\ref{fig:fig6}(c) which shows that the total number of interacting magnons, given by $\sum_{i=1}^N \bra{\Psi(t)}\hat{n}_i\ket{\Psi(t)}$, is  reduced during the time evolution. Therefore, this leads to fewer magnon-magnon scattering events which facilitates accurate tracking of nonequilibrium dynamics of localized spin for longer times by truncated HP transformation. The lost magnons in  Fig.~\ref{fig:fig6}(c) are absorbed by the electronic subsystem and mediate transfer of spin angular momentum between the two subsystems, in accord with the conservation of total $z$-spin shown in Eq.~\eqref{eq:sz_sym}.

Figures~\ref{fig:fig6}(d) and ~\ref{fig:fig6}(e), as the counterpart of Figs.~\ref{fig:fig6}(a) and ~\ref{fig:fig6}(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.~\eqref{eq:hp_resum}. That is, when $N_\mathrm{max} = 3$ is used in Eq.~\eqref{eq:hp_resum},  there is disagreement between two calculations of
$S_1^z(t)$---resummed HP transformation (magenta solid line) vs. the exact one (black solid line)---but increasing $N_\mathrm{max} = 5$ ensures that both methods match exactly.

Thus, Figures~\ref{fig:fig6}(d) and ~\ref{fig:fig6}(e) with proper $N_\mathrm{max} = 5$ provide motivation to derive electron-boson Hamiltonian that
{\em exactly} maps the original electron--localized-spin Hamiltonian in Eq.~\eqref{eq:ham_sys}. The former is required for equilibrium or nonequilibrium MBPT~\cite{Stefanucci2013,Schlunzen2020} which can handle~\cite{Mahfouzi2014} systems in two- or three-dimensions composed of large number of sites  $N \gg 1$ which otherwise renders exact diagonalization~\cite{Wang2019} (or methods like time-dependent density matrix renormalization group~\cite{White2004,Schmitteckert2004,Daley2004,Feiguin2011} tailored for 1D systems) inapplicable. The {\em exact} mapping of the Hamiltonian in Eq.~\eqref{eq:ham_sys} to its electron-boson version
\begin{equation}\label{eq:exacteb}
\hat{H} = \hat{H}_e + \hat{H}_m + \hat{H}_{m-m} + \hat{H}_{e-m}
\end{equation}
using the resummed HP transformation is given by
\begin{widetext}
\begin{subequations}\label{eq:em_ham}
\begin{eqnarray}
    \hat{H}_e
    &=&
    -\gamma\sum_{\braket{ij}}\hat{\psi}^\dagger_i\hat{\psi}_j, \\
    \hat{H}_m
    &=& 2SJ_\mathrm{H}\sum_{i=1}^N\hat{a}_i^\dagger \hat{a}_i  -J_\mathrm{H}S \sum_{\braket{ij}} \big(\hat{a}_i^\dagger \hat{a}_j + \hat{a}_j^\dagger \hat{a}_i \big), \\
    \hat{H}_{m-m}
    &=& -J_\mathrm{H} \sum_{\braket{ij}}(\hat{a}_i^\dagger \hat{a}_i)(\hat{a}_j^\dagger \hat{a}_j)
    \nonumber \\
    &&- J_\mathrm{H}S\sum_{\braket{ij}}\sum_{n=0}^{N_\mathrm{max}}\sum_{m=0}^{N_\mathrm{max}}C_nC_m(1 - \delta_{n0}\delta_{m0})\big[
    (\hat{a}_i^\dagger)^n(\hat{a}_i)^n\hat{a}_i\hat{a}_j^\dagger (\hat{a}_j^\dagger)^m (\hat{a}_j)^m + \hat{a}_i^\dagger (\hat{a}_i^\dagger)^n(\hat{a}_i)^n(\hat{a}_j^\dagger)^m(\hat{a}_j)^m\hat{a}_j \big], \\
    \hat{H}_{e-m} &=& \sqrt{2}S J_\mathrm{sd}\sum_{i=1}^N \sum_{n=0}^{N_\mathrm{max}}C_n\big[\hat{\psi}_{i\uparrow}^\dagger \hat{\psi}_{i\downarrow} \hat{a}_i^\dagger(\hat{a}_i^\dagger)^n(\hat{a}_i)^n +  \hat{\psi}_{i\downarrow}^\dagger \hat{\psi}_{i\uparrow} (\hat{a}_i^\dagger)^n(\hat{a}_i)^n \hat{a}_i\big].
\end{eqnarray}
\end{subequations}
\end{widetext}
The physical meaning of the terms in Eq.~\eqref{eq:em_ham} quite transparent: $\hat{H}_e$ is the tight-binding Hamiltonian of noninteracting electrons; $\hat{H}_m$ is the  Hamiltonian of noninteracting HP bosons; $\hat{H}_{m-m}$ describes various many-body magnon-magnon interaction terms where the first term resembles a nearest-neighbor Hubbard interaction whereas the second term describes other multi-boson interactions captured by the resummed HP transformation; and $\hat{H}_{e-m}$ describes the electron-boson interactions where absorption and emission of HP bosons accompanies electron spin flips as the spin angular momentum is transferred. We note that Eq.~\eqref{eq:em_ham} is much more complex that what is typically used in spintronics literature~\cite{Mahfouzi2014, Tveten2015, Zheng2017, Bender2019, Troncoso2019, Kamra2016}.  Most importantly, it shows that MBPT or diagrammatic Monte Carlo simulations~\cite{Bertrand2019,Bertrand2019a} for interacting electron-magnon system must deal with nonlinear~\cite{Marini2018} electron-boson interactions. Nevertheless, resummed HP transformation makes such MBPT easier and more accurate than truncated HP transformation requiring very large number of terms $N_T$ in order to keep track of localized spin dynamics [Fig.~\ref{fig:fig5}] accurately over sufficiently long times encountered in spintronics phenomena.

\subsection{Ground state of electron-magnon interacting system}\label{sec:groundstate}
\begin{figure*}
	\centering
	\includegraphics[width=\linewidth]{fig7.pdf}
	\caption{(a) The composition of {\em exact} GS $\ket{\Psi_0}$ of 1D quantum many-body system, illustrated in Fig.~\ref{fig:fig1}(b) and composed of \mbox{$N=3$} sites hosting spin-1 localized spins, is computed by exact diagonalization using localized spin operators in Hamiltonian in Eq.~\eqref{eq:ham_sys}. In every ket in the sum, red arrows depict the 	quantum state of localized spins whereas blue arrows indicate the quantum state of $N_\mathrm{e} = 3$ electrons. Here, \mbox{$\epsilon_0 = -\sqrt{2}$ eV}, \mbox{$\epsilon_1 = 0$ eV}, and   \mbox{$\epsilon_2 = \sqrt{2}$ eV} are eigenenergies of the electronic noninteracting single-particle Hamiltonian [Eq.~\eqref{eq:eham_diag}] and \mbox{$J_\mathrm{H} = J_\mathrm{sd} = 0.2$ eV}. (b) The composition of GS evaluated by the truncated HP transformation  (Sec.\ref{sec:hptruncated}) with a truncation number $N_T = 6$ for the same physical parameters as in panel (a). Panel (c) shows that to match the composition of exact GS from panel (a) requires $N_T = 20$ in truncated HP transformation. Numbers on the top of each ket are coefficients in their linear superposition to obtain $|\Psi_0\rangle$.}
	\label{fig:fig7}
\end{figure*}

It is also interesting to understand the structure of exact quantum many-body ground state (GS) of conduction electrons plus localized spins in {\em equilibrium} (GS implies zero temperature as well), described by Hamiltonian in Eq.~\eqref{eq:ham_sys} in terms of original localized spin operators; as well as how many terms of truncated HP transformation (Sec.\ref{sec:hptruncated}) is required to obtain the same ground state by exactly diagonalizing Hamiltonian of the same system expressed in terms of electronic and bosonic operators. In this Section, such analysis is performed for  the system depicted in Fig.~\ref{fig:fig1}(b) which is composed of $N=3$ sites hosting spin-$1$ localized spins, while in Secs.~\ref{sec:spectralelectron} and  \ref{sec:spectralmagnon} we also perform the same comparison for electronic and magnonic spectral functions, respectively, which require additional
information of {\em excited} quantum many-body states.

The exact GS is obtained in three steps: (\emph{i})  $\hat{H}$ is represented as a matrix in the basis of eigenstates of $\hat{N}_\mathrm{e}$ and $\hat{S}^z_\mathrm{tot}$ to render a block-diagonal matrix as shown in Fig.~\ref{fig:fig2}(c); (\emph{ii}) to ensure half-filling for electrons, the matrix block corresponding to $N_\mathrm{e} = 3$ electrons is isolated; (\emph{iii}) This matrix block is diagonalized and the eigenstate with the lowest eigenenergy $E_0$ is identified as the GS. Obviously, if in step (\emph{i}) $\hat{H}$ is expressed directly in terms of localized spin operators [Eq.~\eqref{eq:spin_mat}], then step-(\emph{iii}) yields the \emph{numerically exact}  GS $\ket{\Psi_0}$. On the other hand, if $\hat{H}$ is expressed using the truncated HP transformation with a truncation number $N_T$ [Eq.~\eqref{eq:split_ham}], then thus obtained GS $\ket{\Psi_0}^\mathrm{HP}_\mathrm{N_T}$ is not guaranteed to be the same as $\ket{\Psi_0}$. In particular, we are interested to know what value of $N_T$ ensures that $\ket{\Psi_0}^\mathrm{HP}_\mathrm{N_T} = \ket{\Psi_0}$.

Figure~\ref{fig:fig7}(a) depicts the numerically exact GS $\ket{\Psi_0}$ as a linear combination~\cite{Frederiksen2004} of many-body kets  where red arrows  denotes quantum state of localized spins (using the same notation as introduced in Sec.~\ref{sec:magnons} for spin-1 localized spins) while blue arrows denote  \mbox{spin-$\uparrow$} or \mbox{spin-$\downarrow$} electrons filling the single-particle energy levels \mbox{$\epsilon_0 = -\sqrt{2}$ eV}, $\epsilon_1 = 0$, and \mbox{$\epsilon_2 = \sqrt{2}$ eV} of noninteracting tight-binding Hamiltonian $\hat{H}_e$ [Eq.~\eqref{eq:elec_ham}]. In contrast, we immediately notice in Fig.~\ref{fig:fig7}(b) that the interacting many-body GS $\ket{\Psi_0}^\mathrm{HP}_\mathrm{N_T=6}$ evaluated using truncated HP transformation with $N_T = 6$ is entirely different from $\ket{\Psi_0}$ shown in Fig.~\ref{fig:fig7}(a).  Only when the truncation number is increased to $N_T = 20$ in Fig.~\ref{fig:fig7}(c) we find $\ket{\Psi_0}^\mathrm{HP}_\mathrm{N_T=20} = \ket{\Psi_0}$.

It is worth examining further the structure of the exact GS $\ket{\Psi_0}$ in Fig.~\ref{fig:fig7}(a). Its many-body eigenenergy is \mbox{$E_0 = -3.33$ eV} while the other quantum numbers [Eq.~\eqref{eq:qnumbers}] are $N_\mathrm{e} = 3$ and $S^z_\mathrm{tot} = 3.5$. The largest contribution (greater than $99\%$) to $\ket{\Psi_0}$ comes from the first term on the RHS in Fig.~\ref{fig:fig7}(a) where $N_\mathrm{e} = 3$ electrons fill up the single-particle energy levels  $\epsilon_0$, $\epsilon_1$, and $\epsilon_2$ of noninteracting Hamiltonian $\hat{H}_\mathrm{e}$ [Eq.~\eqref{eq:elec_ham}] 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.~\ref{fig:fig7}(a)] where $N_\mathrm{mag} = 1$ HP boson is created on one of the three sites. This HP boson is actually emitted when the spin-$\downarrow$ electron in eigenenergy level $\epsilon_0$  undergoes a spin-flip process and emerges as a spin-$\uparrow$ electron in eigenenergy level $\epsilon_2$. This process respects conservation of total \mbox{$z$-spin} as encoded in Eq.~\eqref{eq:sz_sym}. The final four kets on the RHS of Fig.~\ref{fig:fig7}(a) are purely electronic excitations where electrons  are excited
among eigenenergy levels $\epsilon_0$, $\epsilon_1$, and $\epsilon_2$ but are not accompanied by any spin-flip process of localized spins.

\subsection{Electronic spectral function for electron-magnon interacting system} \label{sec:spectralelectron}

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{fig8.pdf}
	\caption{(a) The exact electronic spectral function $A(E)$ (red lines) in Eq.~\eqref{eq:ele_spec_f} for 1D quantum many-body system,  illustrated in Fig.~\ref{fig:fig1}(b) and  comprised  of $N=3$ sites hosting spin-1 localized spins, is evaluated using localized spin operators in Hamiltonian in Eq.~\eqref{eq:ham_sys}. This is compared to $A(E)$ in panels (b) and (c) evaluated by mapping localized spin operators to HP bosonic operators via truncated HP transformation  [Eq.~\eqref{eq:tay}] with truncation numbers $N_T= 6$ and $N_T=20$, respectively. Panels (d)--(f) show counterpart information to panels (a)--(c), but for spin-$\frac{5}{2}$ localized spins on each of $N=3$ sites. In all panels \mbox{$J_\mathrm{H} = J_\mathrm{sd} =0.2$ eV}, except for black curves in panels (a) and (d) for which electron--localized-spin interaction is turned off ($J_\mathrm{sd} =0$ while keeping \mbox{$J_\mathrm{H} =0.2$ eV}).}
	\label{fig:fig8}
\end{figure}
\begin{figure*}
	\centering
	\includegraphics[width=\linewidth]{fig9.pdf}
	\caption{The exact magnonic spectral function $D(E)$ in Eq.~\eqref{eq:mag_spec_f} evaluated by mapping localized spin operators to HP bosons in a numerically exact fashion via Eq.~\eqref{eq:hp} for 1D quantum many-body system  illustrated in Fig.~\ref{fig:fig1}(b) comprised  of $N=3$ sites
	hosting (a),(c) spin-$1$ or (b),(d) spin-$5/2$ localized spins. The electron--localized-spin interaction is turned off  (\mbox{$J_\mathrm{sd} = 0$ eV}) in
	panels (a) and (b), or turned on (\mbox{$J_\mathrm{sd} = 0.2$ eV}) in panels (c) and (d). In panels (c) and (d), dotted circle mark additional peaks in $D(E)$ originating excited states $|\Psi_k\rangle$ entering into Eq.~\eqref{eq:mag_spec_f}. (e) Composition of selected excited quantum many-body state $|\Psi_k\rangle$ for spin-$1$ case which is responsible for peak in $D(E)$ near \mbox{$E=1.6$ eV} in panel (c).  The state $|\Psi_k\rangle$ is superposition of kets with zero magnons and solely electronic excitations; followed by kets with one magnon; and then states with two magnons. Numbers on the top of each ket are coefficients in their linear superposition to obtain $|\Psi_k\rangle$.}
	\label{fig:fig9}
\end{figure*}

For the same system considered in Sec.~\ref{sec:groundstate}, Fig.~\ref{fig:fig8}(a)--(c) compares the electronic spectral function $A(E)$ [Eq.~\eqref{eq:ele_spec_f}] evaluated from truncated HP transformation vs. the exact one evaluated using localized spin operators. To set a reference point,
in Fig.~\ref{fig:fig8}(a) we first consider the noninteracting electronic spectral function (black line) for electron--localized-spin interaction
turned off ($J_\mathrm{sd} = 0$). For such a case, the available single-particle states are simply the eigenstates of the noninteracting tight-binding  Hamiltonian  $\hat{H}_\mathrm{e}$ [Eq.~\eqref{eq:elec_ham}] with single-particle energy levels $\epsilon_0$, $\epsilon_1$, and $\epsilon_2$, so that  $A(E)$  consists of sharp peaks centered at $\epsilon_0$, $\epsilon_1$, and $\epsilon_2$. For electron--localized-spin interaction  turned on  (\mbox{$J_\mathrm{sd} = 0.2$ eV}), the exact $A(E)$  (red line) evaluated using localized spin operators  is modified by exhibiting peak splitting [with respect to black line reference for $J_\mathrm{sd} = 0$] at energies $\epsilon_0$, $\epsilon_1$, and $\epsilon_2$. Also, few additional peaks around single-particle energy levels $\epsilon_0$ and $\epsilon_2$ emerge.

In Fig.~\ref{fig:fig8}(b), we compute $A(E)$ using truncated HP transformation with a truncation number $N_T=6$. Although it reproduces the peak-splitting near $\epsilon_0$,  $\epsilon_1$ and  $\epsilon_2$, it exhibits several additional peaks that are absent in the exact result in Fig.~\ref{fig:fig8}(a). This discrepancy can be understood as follows. A correct computation of $A(E)$ requires exact interacting many-body GS $\ket{\Psi_0}$ in Eq.~\eqref{eq:ele_spec_f} and Eq.~\eqref{eq:ele_weight}. However, for  $N_T=6$ the interacting many-body GS $\ket{\Psi_0}^\mathrm{HP}_\mathrm{N_T = 6} \neq \ket{\Psi_0}$ as shown by Fig.~\ref{fig:fig7}(b). At first sight, it appears that the same argument should produce exact $A(E)$ in Fig.~\ref{fig:fig8}(c) using $N_T=20$ because
$\ket{\Psi_0}^\mathrm{HP}_\mathrm{N_T = 20} = \ket{\Psi_0}$ in Fig.~\ref{fig:fig7}(c). However, the discrepancy between exact $A(E)$ (red line) in Fig.~\ref{fig:fig8}(a) and blue line in Fig.~\ref{fig:fig8}(c) can be explained by inspecting  Eqs.~\eqref{eq:ele_spec_f} and ~\eqref{eq:ele_weight}. A correct evaluation of $A(E)$ not only requires usage of exact GS $\ket{\Psi_0}$, but it also requires correct interacting many-body excited states $\ket{\Psi_k}$ as shown by Eq.~\eqref{eq:ele_weight} for which a truncation number $N_T = 20$ becomes  {\em insufficient}. We also repeated analysis from Fig.~\ref{fig:fig8}(a)--(c) for spin-$1$ localized spins by using spin-$\frac{5}{2}$ localized spins in Fig.~\ref{fig:fig8}(d)--(f), which shows that our conclusions are not limited to specific value of chosen localized spins.

\subsection{Magnonic spectral function for electron-magnon interacting system} \label{sec:spectralmagnon}

The {\em exact} HP transformation in Eq.~\eqref{eq:hp} makes it possible to define magnonic spectral function $D(E)$ in Eq.~\eqref{eq:mag_spec_f} and compute it without any approximations numerically by evaluating square roots of matrices in Eq.~\eqref{eq:hp}. Using the same systems of spin-$1$ or spin-$5/2$ localized spins studied in Fig.~\ref{fig:fig8}, we first establish a reference point by computing $D(E)$ in Figs.~\ref{fig:fig9}(a) and \ref{fig:fig9}(b)
with electron--localized-spin interaction turned off ($J_\mathrm{sd} = 0$). Such $D(E)$ in Fig.~\ref{fig:fig9}(a) exhibits three peaks at energies $E=0$ eV, $E = 0.2$ eV, and $E=0.6$ eV which correspond to the available states in the presence of solely localized-spin--localized-spin (or equivalently magnon-magnon) interactions. Conversely, when we turn on \mbox{$J_\mathrm{sd} = 0.2$ eV} in Figs.~\ref{fig:fig9}(c) and ~\ref{fig:fig9}(d), we find: (\emph{i}) the original noninteracting peaks are largely the same, except for the one near \mbox{$E=0.2$ eV} which undergoes a tiny splitting; (\emph{ii}) far away from these peaks, $D(E)$ acquires additional peaks (marked by dotted circles) near energies \mbox{$E=1.6$ eV} and \mbox{$E=2.1$ eV}. Analogous features are observed for spin-$5/2$ localized spins when switching from \mbox{$J_\mathrm{sd} = 0$} in Fig.~\ref{fig:fig9}(b) to \mbox{$J_\mathrm{sd} \neq 0$} in Fig.~\ref{fig:fig9}(d).

Similar additional peaks in magnonic spectral function were previously found in MBPT calculations~\cite{Mahfouzi2014} using summation of infinite class of selected diagrams and were interpreted in that language as quasibound states of magnons dressed by the cloud of electron-hole pair excitations. Also, MBPT calculations of Refs.~\cite{Mahfouzi2014,Woolsey1970} find small modification of electronic $A(E)$ vs. large modification of magnonic $D(E)$ upon tuning electron-magnon interactions on, which is in accord with modification in Fig.~\ref{fig:fig8} vs. that in Fig.~\ref{fig:fig9} when switching from  \mbox{$J_\mathrm{sd} = 0$} to  \mbox{$J_\mathrm{sd} \neq 0$}. This is explained by magnons being in the strongly interacting regime vs. electrons being in the weakly interacting regime due to~\cite{Mahfouzi2014} $J_\mathrm{sd}$ divided by the bandwidth of noninteracting magnons being much larger than $J_\mathrm{sd}$ divided by the bandwidth of noninteracting electrons.

To clarify the origin of these peaks further in the context of our nonperturbative calculations in Figs.~\ref{fig:fig9}(c) and ~\ref{fig:fig9}(d), we focus on the peak near \mbox{$E=1.6$ eV} in Fig.~\ref{fig:fig9}(c). This peak is due to quantum many-body excited state $\ket{\Psi_k}$ whose composition is given explicitly in Fig.~\ref{fig:fig9}(e). This state has a nonzero ``weight''  $Q_k^+ = 0.006$ in Eq.~\eqref{eq:mag_weight}. Although, the value of $Q_k^+$ appears to be small, and it contributes only about $2\%$ in the ``sum rule'' in Eq.~\eqref{eq:mag_sum_rule}, it cannot be ignored.  Interestingly, Fig.~\ref{fig:fig9}(e) reveals that this specific $\ket{\Psi_k}$ is a linear superposition of states with $N_\mathrm{mag} = 0$, $1$, or $2$ HP bosons.

\subsection{Diagonal and off-diagonal elements of time-dependent electronic and bosonic lesser Green functions}\label{sec:diag_offdiag}

The time-dependent electronic lesser GF $G^<_{i\sigma, j\sigma'}(t, t')$ in Eq.~\eqref{eq:ele_GF_less} generally depends on two times $t$ and $t'$~\cite{Stefanucci2013,Schlunzen2020}. At equal times $t'=t$, it yields one-particle electronic density matrix~\cite{Stefanucci2013,Gaury2014,Bajpai2019b}
\begin{equation}
    \bm{\rho}(t) = -i\hbar \bold{G}^<(t, t')\big|_{t'=t}.
\end{equation}
Its diagonal elements in, e.g., coordinate (or site for discrete lattice) representation contain information about the time-dependent electronic charge and spin density~\cite{Petrovic2018}, whereas the off-diagonal elements encode quantum-mechanical interference effects~\cite{Bajpai2019b} and measure the degree of quantum coherence~\cite{Schlosshauer2005}. The illustrate their time evolution, we use the same 1D quantum many-body system employed in Fig.~\ref{fig:fig6} where localized spin-$5/2$ at site  $i=1$ is completely flipped [i.e., $N_\mathrm{mag} = 5$ HP bosons are introduced at site $i=1$ via Eq.~\eqref{eq:fivehpbosons}] to initiate nonequilibrium dynamics. Fig.~\ref{fig:fig10}(a)--(d) shows the ensuing time-evolution of diagonal elements, $-i\hbar G^<_{i\sigma, i \sigma}(t, t)$; and off-diagonal elements,  $-i\hbar G^<_{i\sigma, j \sigma'}(t, t)$, for the same electron--localized-spin 1D quantum many-body system whose equilibrium properties are analyzed Fig.~\ref{fig:fig8}. For reference, we turn the electron--localized-spin interaction off ($J_\mathrm{sd} = 0$) in Figs.~\ref{fig:fig10}(a) and ~\ref{fig:fig10}(b), which trivially leads to all elements being time independent because for $J_\mathrm{sd} = 0$ the quantum state of the electronic subsystem is an eigenstate of the electronic Hamiltonian $\hat{H}_\mathrm{e}$ [Eq.~\eqref{eq:elec_ham}].
\begin{figure*}
	\centering
	\includegraphics[width=\linewidth]{fig10.pdf}
	\caption{Panels (a) and (b) show time-evolution of diagonal and off-diagonal elements (absolute value of latter), respectively, of the electronic lesser GF $G^<_{i\sigma, j \sigma'}(t, t)$ [Eq.~\eqref{eq:ele_GF_less}] in the site basis for the system depicted in Fig.~\ref{fig:fig1}(b) composed of $N=3$ sites hosting spin-$\frac{5}{2}$ localized spins but electron-magnon interaction is turned-off ($J_\mathrm{sd} = 0$). Panels (c) and (d) are its counterpart when the electron-magnon interaction is turned-on (\mbox{$J_\mathrm{sd} = 1$ eV}). Panels (e)--(h) are same as (a)--(d), but for bosonic lesser GF $D^<_{ij}(t, t)$ [Eq.~\eqref{eq:mag_GF_less}] in the site basis. Note that some of the off-diagonal elements are not shown because they are either zero or identical to the ones shown in panels (b),(d) and panels (f),(h).}
	\label{fig:fig10}
\end{figure*}
Conversely, in Figs.~\ref{fig:fig10}(c) and ~\ref{fig:fig10}(d) we use \mbox{$J_\mathrm{sd} = 1$ eV} which leads to nontrivial time dependence of  diagonal and off-diagonal elements of $-i\hbar\bold{G}^<(t, t)$. Interestingly, the diagonal elements $-i\hbar G^<_{i\sigma, i\sigma}(t, t)$ in Fig.~\ref{fig:fig10}(c)  satisfy $-i\hbar[G^<_{i\uparrow, i\uparrow}(t, t) + G^<_{i\downarrow, i\downarrow}(t, t)] = Q_i$, where $Q_i$ the total charge density on site $i$ which is a constant in time. This means that no charge currents flows between sites $i$ and $j$. Instead, population of electrons with spin \mbox{$\sigma=\uparrow,\downarrow$} at site $i$ exchanges solely between spin \mbox{$\sigma=\uparrow, \downarrow$} states at that site. This is also accompanied by time-evolution of off-diagonal elements $-i\hbar G^<_{i\sigma, j\sigma'}(t, t)$ in Fig.~\ref{fig:fig10}(d).

The off-diagonal elements of lesser GF are also crucial for calculation of the so-called many-body lesser self-energy  $\bm{\Sigma}^<(t_1, t_2)$~\cite{Stefanucci2013,Schlunzen2020} satisfying the Keldysh equation
\begin{equation}\label{eq:gless_sc}
    \bold{G}^<(t, t') =
    \int\displaylimits_{-\infty}^{+\infty}
    \int\displaylimits_{-\infty}^{+\infty}
    dt_1 dt_2
    \bold{G}^r(t, t_1)\bold{\Sigma}^<(t_1, t_2)\bold{G}^a(t_2, t'),
\end{equation}
which encapsulates time evolution of quantum many-body systems in terms of solely one-particle quantities. A self-consistent solution to Eq.~\eqref{eq:gless_sc} can yield exact many-body lesser self-energy. Alternatively, one can systematically approximate it~\cite{Mahfouzi2014}
using conserving approximations of MBPT. One such conserving approximation for the lesser self-energy of electron-boson interacting systems
is the so-called self-consistent Born approximation (SCBA)~\cite{Mahfouzi2014,Frederiksen2007,Lee2009,Mera2016}. The SCBA ensures charge
conservation in nonequilibrium~\cite{Frederiksen2004}, and in steady-state nonequilibrium one can Fourier transform $\bm{\Sigma}^<(t_1, t_2)$
to work with energy-dependent~\cite{Luiser2009,Rhyner2014}
\begin{eqnarray}
    \Sigma_{ij}^<(E) &=&
    \sum_{\alpha\beta}\sum_{\lambda q}
     M^\lambda_{i \alpha}(\bold{q}) \Big\{ n_\lambda(\bold{q}) G^<_{\alpha \beta}[E \pm \hbar\omega_\lambda(\bold{q})] \nonumber \\
    &&+ [n_\lambda(\bold{q}) + 1] G^<_{\alpha \beta}[E \mp \hbar\omega_\lambda(\bold{q})] \Big\} M^{\lambda *}_{\beta j}(\bold{q}).
\end{eqnarray}
Here, $n_\lambda(\bold{q}) = [e^{\hbar\omega_\lambda(\bold{q})/k_bT} - 1]^{-1}$ is the Bose-Einstein distribution function for boson-frequency $\omega_\lambda(\bold{q})$ of a boson-mode $\lambda$; $\bold{q}$ is the wave-vector for bosons; $\Sigma_{ij}^<(E)$ and $G^<_{ij}(E)$ represent $2\times 2$ matrix blocks in electron spin space of $\sigma=\uparrow, \downarrow$ between site $i$ and $j$; and $M_\mathrm{i\alpha}^\lambda(\bold{q})$ is defined in Ref.~\cite{Luiser2009}.

To reduce computational complexity~\cite{Frederiksen2007} of calculations of $\Sigma^<_{ij}(E)$ and enable simulations of devices containing large number of atoms, the ``local self-energy approximation'' is often employed~\cite{Luiser2009,Rhyner2014,Cavassilas2016,Bescond2018}. In this approximation, one assumes $|\Sigma^<_{ii}(E)| \gg |\Sigma^<_{ij}(E)|$, i.e., that off-diagonal elements of self-energy are minuscule when compared to the diagonal ones and can, therefore, be set to zero. This is done in conjunction with discarding the off-diagonal elements of the electronic lesser GF, $G^<_{ij}(E) \approx G^<_{ii}(E)\delta_{ij}$. Thus,
which furnishes the local self-energy approximation ``local self-energy approximation''~\cite{Luiser2009,Rhyner2014,Cavassilas2016,Bescond2018} leads to
\begin{eqnarray}
    \Sigma_{ii}^<(E) & \approx &
    \sum_{\alpha}\sum_{\lambda q}
     M^\lambda_{i \alpha}(\bold{q}) \Big\{ n_\lambda(\bold{q}) G^<_{\alpha \alpha}[E \pm \hbar\omega_\lambda(\bold{q})] \nonumber \\
    &&+ [n_\lambda(\bold{q}) + 1] G^<_{\alpha \alpha}[E \mp \hbar\omega_\lambda(\bold{q})] \Big\} M^{\lambda *}_{\alpha j}(\bold{q}).
\end{eqnarray}

Consequently, one may be tempted to follow SCBA along with the local self-energy approximation even for electron-magnon interacting systems, however, our numerically exact calculation of electronic lesser GF shows that the average value of off-diagonal elements in Fig.~\ref{fig:fig10}(d) is not minuscule, and instead is roughly one-fifth of the diagonal elements in Fig.~\ref{fig:fig10}(c) (see $y$-axis scale). Therefore, ``local approximation'' for the self-energy and nonequilibrium GFs cannot be justified in the case of many-body electron--localized-spin interacting systems.  Lastly, in Figs.~\ref{fig:fig10}(e)--(h), we show the counterpart of Figs.~\ref{fig:fig10}(a)--(d) but for the bosonic GFs. Here we note that the off-diagonal elements of the lesser GF are always comparable to the diagonal ones whether or not electron-magnon interaction is turned on.

\section{Conclusions}\label{sec:conclusions}
Using numerically exact diagonalization, we study 1D finite size quantum many-body systems composed of $N \leq 7$ sites hosting spin-$S$ localized spins ($S=1$ or $S=5/2$) that interact with each other via NN exchange interaction of strength $J_\mathrm{H}$ along with conduction electrons spin via $sd$ exchange interaction of strength $J_\mathrm{sd}$; and investigate the range of validity of the truncated HP transformation [Eq.~\eqref{eq:hp}] by comparing to calculations from \emph{numerically exact} localized spin operators. Our main results can be summarized as the following:
\begin{enumerate}
\item  For a quantum many-body system composed of localized spins alone ($N_\mathrm{e} = 0$ electrons), Fig.~\ref{fig:fig4} shows that as more interacting HP bosons are introduced in the system, progressively larger value of truncation number $N_T$ is required to incorporate multi-magnon interaction terms and accurately predict the nonequilibrium dynamics of localized spins.

\item Figure~\ref{fig:fig5} shows that the breakdown-time $t_\mathrm{break}$ for truncated HP transformation power-series with a truncation number $N_T$ follows $t_\mathrm{break}\propto e^{N_T/C}$. Although, the exponential dependence of $t_\mathrm{break}$ on the truncation number $N_T$ is favorable, its typical value for low-order values of $N_T = 1$~to~$5$ is roughly \mbox{$t_\mathrm{break} \approx 1.5 \hbar/J_H \approx 0.1$ fs} ($J_\mathrm{H} = 1$ eV) which shows its poor performance on femtosecond time-scale.

\item Figures~\ref{fig:fig6}(a)--(c) show that when we introduce electron-magnon interaction of localized spins with conduction electron spins, the truncation number $N_T$ required to accurately predict nonequilibrium dynamics of localized spins is reduced due to transfer of spin angular momentum between the two subsystems, which effectively reduces the total number of interacting magnons within the localized spin subsystem. On the other hand, Figs.~\ref{fig:fig6}(d) and (e), show that the resummed HP transformation in Eq.~\eqref{eq:hp_resum} is numerically exact when $N_\mathrm{max}  = 2S$. The corresponding interacting electron-magnon Hamiltonian furnished by it in Eq.~\eqref{eq:em_ham} is more elaborate than tradition low-order truncated HP transformation of localized spin operators.

\item Figure~\ref{fig:fig7} shows that the truncated HP transformation with a truncation number $N_T = 6$ [Fig.~\ref{fig:fig7}(b)] produces an incorrect GS of the interacting electron-magnon system and only when $N_T = 20$ [Fig.~\ref{fig:fig7}(c)], it is able to match the correct interacting GS evaluated using localized spin operators [Fig.~\ref{fig:fig7}(a)]. As a consequence, the interacting electronic spectral function in Figs.~\ref{fig:fig8}(b) and (c) [(e) and (f)] for spin-1 localized spins [spin-$\frac{5}{2}$ localized spins] does not match the exact one in Fig.~\ref{fig:fig8}(a) [Fig.~\ref{fig:fig8}(d)].

\item In Fig.~\ref{fig:fig9}, we show that the interacting magnonic spectral function acquires substantial modification from its noninteracting counterpart in the form of additional peaks that are due to excited states $\ket{\Psi_k}$ whose structure is also shown.

\item In Fig.~\ref{fig:fig10}, we showcase the time-evolution of the diagonal and off-diagonal elements of the lesser GF (electronic and bosonic) which shows that the magnitude of off-diagonal elements of the electronic lesser GF is comparable to the diagonal ones, hence, they cannot be ignored as done in the so-called ``local self-energy approximation" which is routinely employed in electron-phonon calculations within SCBA.

\end{enumerate}
\begin{acknowledgments}
This research was supported by the US National Science Foundation (NSF) under Grant No. ECCS-1922689. This paper has originated from ``Research Projects Based Learning'' implemented within a graduate level course {\em PHYS814: Advanced Quantum Mechanics}~\cite{phys814}  at the University of Delaware.
\end{acknowledgments}

redpajama_set_name: RedPajamaArXiv