\section{Introduction}
In the research of condensed matter~\cite{mcmillan2004condensed}, pressure is a fundamental thermodynamic variable that determines the state of matter, and plays an important role in  the field. The discovery of new materials and thence application to industrial use makes up an important part of modern innovation. As a basic thermodynamic parameter, pressure demonstrates the ability to activate semi-core electrons, empty orbitals, and non-atom-centered quantum orbitals on interstitial sites, changing a given element's chemistry, and thus resulting in a plethora of novel and previously unexpected occurrences, for example the generation of new types of functional materials that deviate from those under atmospheric pressure (101.325 kPa)~\cite{miao2020chemistry}.

Currently, hydrogen-rich materials have garnered much attention in regards to obtaining superconductivity at high temperatures, and has led to much theoretical and experimental work on the search for high-temperature superconductivity of hydrides under high pressure~\cite{wang2018hydrogen}. The chemical pre-compression method, proposed by Neil Ashcroft in 2004, has led to the key realization that hydrogen-rich compounds are a new potential class of high-temperature superconductors~\cite{ashcroft2004hydrogen}.  As investigating the high-temperature superconductivity of metallic hydrogen is highly challenging~\cite{ashcroft1968metallic}, most scholars have instead redirected focus to the synthesis and properties of hydride-rich compounds instead. Recently, the realization of very high temperature superconductivity, at near room temperature, was discovered in hydrogen disulfide~\cite{nakao2019superconductivity,harshman2017compressed} and lanthanum hydrogen~\cite{sun2021high,kostrzewa2020lah,yi2021stability} an important set of milestones.

Hydrogen disulfide was previously believed to undergo dissociation under high pressures, and was not heavily considered as a potential superconductor. However, recent theoretical work has suggested that the dissociation would not occur, but that on the contrary the material becomes superconducting~\cite{li2014metallization}. Inspired by this discovery, Drozdov et al. compressed sulfur hydride in a diamond anvil cell, and has shown that hydrogen sulphide would indeed form a superconductor when compressed with a $T_c$ of up to 200 K~\cite{drozdov2015conventional}.

Armed with this new knowledge, researchers have expanded their search to include lanthanide hydrides under pressure, with notable examples such as La-H and Y-H. The Fermi density and thus superconductivity in La-H systems are determined not only by the $d$ orbital La electrons and the $s$ orbital H electrons, but also by the $f$ orbital electrons of lanthanum~\cite{sun2020second}.

Using DFT, studies have found that the face-centered cubic (FCC) form of LaH$_{10}$ is a good metal, with various bands crossing the Fermi level that forms a high electronic density at the Fermi level~\cite{liu2019microscopic}.
In contrast, for the YH$_{10}$ system, only the $d$ orbital Y electrons and the $s$ orbital H electrons become major contributors to the Fermi density at high pressures~\cite{heil2019superconductivity}. As external pressure destabilizes the localised La-$4f$ more than the other orbitals (La-$6s$,La-$5d$), the latter populates when pressure is applied, leading to possible novel emergent quantum states associated with the strong electronic correlations.
In particular, it was purported that the La-H system has a unique high-$T_c$ with $f$ electrons at the Fermi level.~\cite{song2020high}.

The effects of electronic localization and hybridization apply for all high-density compounds, and is a hallmark for a wide family of similar materials signified by transition metals, lanthanides, and rare earth elements. Many-body effects give rise to a number of unique and interesting phenomena in similar systems, such as high-temperature superconductivity in cuprates and the metal-insulator transition in vanadates at ambient temperature~\cite{sarma2013many}. The Zaanen-Sawatsky-Allen theory offers a categorization of transition metal periodic solids in terms of their correlation' strengths, and has provided a good knowledge base, though deep understanding of the hydride's characteristics distant from ambient conditions remains hitherto difficult and incomplete~\cite{olalde2011direct}.

Therefore, in order to identify new hydrogen-rich high-$T_c$ superconductors at the lowest possible pressure, the use of quantitative theoretical computations are required. Superconductivity in these potential compounds are mediated by the interaction between highly intense lattice vibrations of hydrogen atoms and the localised electrons. An accurate description of this interaction necessitates fine descriptions of the electronic characteristics, which are infamous for being difficult for such correlated $f$ systems, treating both itinerant and localised electrons on the same footing.

Several theoretical aspects, including the electron-phonon coupling strength $\lambda$, phonon dispersion relations, electron spectral weight, and cross terms between electron-electron, and electron-phonon interactions are corrected with electronic correlations.

Here, we propose a pragmatic first-principles calculation consistent platform that uses many-body corrections for the electronic spectral weight, which then feeds into adjusted estimates of $T_c$. The many-body corrections to phonon dispersion are typically less drastic than the corrections to the electronic spectra, as correlations effects can shift the f- spectral weight over several eV. Furthermore, since their full treatment is beyond reach, the scope of this work is restricted to only correcting the electronic spectra.

In this manuscript, we show that the many-body corrections in $f$ orbital systems could cause significant changes, with spectral weight shifts in the order of one electron volt. We evaluate the precise physical properties of La-H and Ce-H systems, giving special care to the effect of correlations on spectral properties. We find that LaH$_{16}$ and CeH$_{16}$ have stable P6/mmm space group crystal structures at pressures up to $250$ GPa.
Additionally, for the recently discovered Cerium hydride, we predict a comparatively high $T_c$, evaluated using our pre-established hierarchical approach~\cite{plekhanov2021computational}.

\section{Discussion}

The idea that hydrogen-rich compounds could be potential high $T_c$ superconductors dates back to the turn of the millennium, when chemical pre-compression was proposed as a viable way to reduce the metallization pressure of hydrogen in the presence of other elements, leading to observed $T_c$ exceeding 150 K in the LaH$_{16}$ system. This indicates compressed hydrogen-rich compounds as potential room-temperature superconductors.

The genesis of superconductivity in these hydrides is known to originate from electron-phonon interactions. The characteristic phonon frequency, the electron-phonon coupling strength, the density of states at Fermi level, and the Coulomb pseudopotential $\mu^\star$ are the four factors that define $T_c$ according to BCS theory.
Density functional theory, using conventional pseudo-potentials, such as PBE, is widely acknowledged to provide an accurate explanation of lattice dynamics. However, for compounds with weakly hybridized and localized $f$ electrons, such as La and Ce, where many-body corrections are required, DFT is known to have difficulties in dealing with strong electronic correlations.
We apply a density functional theory technique combined with dynamical mean-field theory (DMFT) in this study. The charge and spin local fluctuations, which are important for the local paramagnetic moment of lanthanide elements, are readily corrected by DMFT.
DMFT is used to account for changes in orbital character at the Fermi surface caused by spectral weight transfer related with Hubbard $f$ band splitting. This influences low energy electron-electron scattering processes via phonon momentum transfer, as stated in the Methods section, according to the Allen-Dynes formalism.

P6/mmm-LaH$_{16}$ and P6/mmm-CeH$_{16}$ shows very little dependence on pressure at a hundred GPa. The phonon DOS of P6/mmm-LaH$_{16}$ and P6/mmm-CeH$_{16}$ computed at $250$ GPa is shown in Fig.~\ref{Ce_LaH phonon}, panels (a) and (b). Then, we investigate the impact of various DMFT electronic charge self-consistency schemes (see Fig.~\ref{Ce_LaH phonon}(c) and Fig.~\ref{Ce_LaH phonon}(d)).We find that correlation effects have a significant impact on $\alpha^2F(\omega)$. In particular, we compare: i) PBE density functional theory, ii)  DFT+DMFT with the full charge self-consistent formalism (DFT+DMFT+CSC). We used the Koster-Slater interaction vertex for the La and Ce correlated manyfold, with typical values for $U=6$ eV, $J=0.6$ eV. Interestingly, the full charge self-consistent DMFT provides a large increase of the superconducting temperature
{(see panels (e) and (f) in Fig.~\ref{Ce_LaH phonon})}. This confirms that many-body effects have a sizeable contribution to the prediction of the superconducting temperatures in lanthanide hydrides. Note that we use $U=6$ eV, $J=0.6$ eV throughout the rest of the paper.

\begin{figure*
\includegraphics[width=\textwidth]{graph/LaCeH_phonon_new.pdf}
\caption{\textbf{Many-body corrections to the superconducting temperature}.
We report (a) phonon density of states, (b)Eliashberg function $\alpha^2F(\omega)$, the spectral weight at the Fermi level is obtained at different levels of approximation: i) DFT PBE (black line),  and ii) with the full charge self-consistent formalism (DFT+DMFT+CSC,  red line).
(c) the superconducting temperature $T_{c}$ obtained by the Allen and Dynes formalism. We obtain a theoretical estimate for LaH$_{16}$ is $T_c=166.2$K by DFT and $T_c=219.6$K by DMFT, also a theoretical estimate for CeH$_{16}$ is $T_c=69.7$K by DFT and $T_c=165.1$K by DMFT
 The full charge self-consistent DMFT provides a large increase of the superconducting temperature, and the physical value of the Hund's coupling for La and Ce is$U=6$eV, $J=0.6$eV. All calculations were performed in the P6/mmm phase of LaH$_{16}$ and CeH$_{16}$ at 250 GPa.}
\label{Ce_LaH phonon}
\end{figure*}

We focus on the LaH$_{16}$ and CeH$_{16}$ P6/mmm systems at $250$ GPa. Although DMFT readily provides important corrections to the electronic structure, it is worth investigating how DMFT affects structural qualities. Calculating forces with limited atomic displacement is not tractable due to the high computing overhead of doing many-body adjustments.
Recently, we have developed a method for calculation of forces  within DMFT, allowing for ultra-soft and norm-conserving pseudopotentials in the underlying DFT~\cite{plekhanov2021calculating}. This opens up new possibilities for systems with heavy components that are not well-suited to all-electron computations. The structure relaxation at $250$ GPa is seen in Fig.~\ref{DMFT_interation}. Typically, we obtain corrections for the bond lengths of the order of 5$\%$.

The many-body effects tend to slightly shorten the La-H and Ce-H bonds, somehow increasing the La-H and Ce-H covalency. At the same time, this La-H and Ce-H bond length reduction is accompanied by a moderate increase in the H-H distance. This behaviour is opposite to what we have previously found in CeH$_9$, where both Ce-H and H-H bonds were longer within DMFT treatment~\cite{plekhanov2021computational}.

\begin{figure*
\includegraphics[width=1\textwidth]{graph/DMFT_interation.pdf}
\caption{\textbf{Structural relaxation of clathrate lanthanides with many-body corrections}.
We report the structural relaxation of the LaH$_{16}$ (panel a and c) and prototype CeH$_{16}$ (panel b and d) compounds. All calculations are performed at $250$ GPa. The volume density is obtained by the equation of state in DFT+DMFT+CSC that provides very similar results to PBE (not shown). Internal coordinates are relaxed with DFT+DMFT+CSC, building upon the recent implementation of DFT forces for ultra-soft pseudo-potentials.
We report the forces and total energies obtained during the structural optimization, respectively, in panel (a) (b) for LaH$_{16}$ (CeH$_{16}$). Convergence is obtained within $25$ iterations. The shortest H-H and Ce-H bond lengths increase throughout the structural optimization
(see panel c and d for La and Ce hydrides, respectively).}
\label{DMFT_interation}
\end{figure*}

\begin{figure
\includegraphics[width=1\columnwidth]{graph/LaH16_dmft_new.pdf}
\caption{\textbf{Spectral weight transfer induced by many-body corrections}.(a) Electronic band structure and (b) density of states obtained by DFT calculations. $t DOS$ and $f DOS$ denote the spectral weight obtained by the imaginary part of respectively the lattice and $f$ impurity Green function, corresponding to the spectral weight traced over all orbitals and traced over the $f$ orbitals, respectively. In (c) and (d) we show the energy-resolved spectral weight, obtained respectively by the one-shot DFT+DMFT and the full charge self-consistent DFT+DMFT+CSC. All calculations were performed in the P6/mmm phase of LaH$_{16}$ at 250 GPa.}
\label{LaH16_dmft}
\end{figure}

\begin{figure
\includegraphics[width=1\columnwidth]{graph/CeH16_dmft_new.pdf}
\caption{\textbf{Spectral weight transfer induced by many-body corrections}.(a) Electronic band structure and (b) density of states obtained by DFT calculations. $t DOS$ and $f DOS$ denote the spectral weight obtained by the imaginary part of respectively the lattice and $f$ impurity Green function, corresponding to the spectral weight traced over all orbitals and traced over the $f$ orbitals, respectively. In (c) and (d) we show the energy-resolved spectral weight, obtained respectively by the one-shot DFT+DMFT and the full charge self-consistent DFT+DMFT+CSC. All calculations were performed in the P6/mmm phase of CeH$_{16}$ at 250 GPa.}
\label{CeH16_dmft}
\end{figure}

The changes highlighted above stem directly from a spectral weight transfer induced by many-body corrections (see Fig. \ref{LaH16_dmft}(a) and \ref{LaH16_dmft}(b)). In DFT, the La system is described by a two band system in absence of long-range magnetic order. We note that DFT is a single Slater determinant approach, and hence cannot capture the role of paramagnetism, with an associated magnetic multiplet (fluctuating magnetic moment). Such effects typically induce a splitting of spectral features into satellites, as observed in Figs. \ref{LaH16_dmft}(c) and \ref{LaH16_dmft}(d), with a resulting large increase of $f$ character at the Fermi level. As sharp La features occur near the Fermi level, we emphasize that a higher level of theory is required to capture correctly the superconducting properties. For instance, in our calculations the one-shot (DFT+DMFT) and full charge self-consistent approach (DFT+DMFT+CSC) induces a small shift of the sharp La feature at the Fermi level, which in turns mitigates the $f$ character increase at the Fermi level.

The role of $f$ orbitals appears to be very important for the superconducting properties in rare-earth hydrates~\cite{plekhanov2021computational}. That is why in the present paper, in addition to LaH$_{16}$ with formally empty La $f$ shell, we study also Cerium hydrate (CeH$_{16}$) with one $f$ electron in the atomic Ce configuration. We have checked that both systems remain stable at pressure up to at least $250$ GPa.

We report in Fig. \ref{CeH16_dmft} the DFT+DMFT+CSC framework applied to CeH$_{16}$ in the P6/mmm phase at 250 GPa. We attribute the decrease in $T_c$ to a higher $f$ occupation, which shifts the chemical potential away from the $f$ spectral features present near the Fermi level (see Fig.~\ref{CeH16_dmft}(c) and Fig.~\ref{CeH16_dmft}(d) ).

Our findings point to a possible path for increasing $T_c$ in lanthanide hydrides: an increase in $f$ character at the Fermi level in DFT+DMFT, which is associated with a lower degree of La-H and Ce-H covalency and a lower degree of hybridization, which, in turns, is a marker for a higher superconducting temperature in these systems.

\section{Methods}
\begin{figure
\includegraphics[width=0.9\columnwidth]{graph/workflow1.pdf}
\caption{\textbf{DFT inter-operability for a consistent many-body platform}. Schematic overview of the main modules of the theoretical platform and its interrelations. Firstly, structures are predicted by Crystal structure AnaLYsis by Particle Swarm Optimization (CALYPSO)), via Gibbs enthalpies for the equation of state and convex hull. The underlying core engines are the CASTEP and Quantum Espresso DFT software. Inter-operability between QE and CASTEP is achieved via format conversion of input files, pseudo-potentials and $\mathbf{k}$-point grids.
Core libraries are used to provide the many-body corrections, via quantum embedding, which in turns provides corrected forces and energies. In post-processing, the Eliashberg function and superconducting $T_c$ are obtained with the DMFT+a2F approach. Finally, data are archived for future usage in hdf5.}
\label{workflow1pdf}
\end{figure}

Figure \ref{workflow1pdf} depicts our theoretical approach. We present a schematic overview of the theoretical platform's principal elements and their interrelationships. Our method establishes a modular framework for high-pressure material screening. Firstly, Crystal structure AnaLYsis by Particle Swarm Optimization (CALYPSO) provides stoichiometric compositions via Gibbs enthalpies for the equation of state and convex hull~\cite{wang2012calypso}. Interoperability between Quantum Espresso(QE)~\cite{giannozzi2009quantum,clark2005first} is accomplished through input file format conversion, pseudopotentials, and $\mathbf{k}$-point grids. The many-body corrections are provided by core libraries via the DMFT quantum embedding, which results in corrected forces\cite{plekhanov2021calculating} and total free energies\cite{plekhanov2018many,lee2019mott}.

The underlying structure relaxations were carried out using the QE and CASTEP packages in the framework of DFT and using PBE-GGA (Perdew-Burke-Ernzerhof generalized gradient approximation)~\cite{perdew1996generalized,perdew1992atoms}.
Norm conserving pseudopotential were used to describe the core electrons and their effects on valence orbitals\cite{rappe1991erratum}. Valence electron configuration of $5s^2 5p^6 5d^1 6s^2$, $5s^2 5p^6 4f^1 5d^1 6s^2$ (i.e., with explicitly included $f$ electrons) and $1s^1$ was used for the La, Ce and H atoms, respectively. A plane-wave kinetic-energy cut-off of $1000$ eV  and dense Monkhorst-Pack $\mathbf{k}$-points grids with reciprocal space resolution of $12\times 12\times 12$ were employed in the calculation.

Phonon frequencies and superconducting critical temperature were calculated using density-functional perturbation theory as implemented in QE~\cite{baroni2001phonons}. The $\mathbf{k}$-space integration (electrons) was approximated by a summation over a $20\times 20\times 12$ uniform grid in reciprocal space, with the Methfessel-Paxton smearing scheme, using a temperature of $k_B$T = 0.05 eV for self-consistent cycles and relaxations; the same grid ($20\times 20 \times 12$) was used for evaluating DOS and coupling strength. Dynamical matrices and $\lambda$ were calculated on a uniform $5\times 5\times 3$ grid in $\mathbf{q}$-space for P6/mmm-LaH$_{16}$ and P6/mmm-CeH$_{16}$.

In post-processing, the superconducting transition temperature $T_c$ was estimated using the Allen-Dynes modified McMillan equation~\cite{dynes1972mcmillan}:

\begin{equation}
T_c = \frac{\omega_{\log}}{1.2} \exp \left [
         \frac{-1.04(1+\lambda)}{\lambda - \mu^\star(1+0.62\lambda)}\right ]
\end{equation}
where $\mu^{*}$ is the Coulomb pseudopotential. The electron-phonon coupling strength $\lambda$ and $\omega_{\log}$ were calculated as:
\begin{equation}
\omega_{\log} = \mbox{exp} \left [ \frac{2}{\lambda} \int \frac{d\omega}{\omega}
                                  \alpha^2F(\omega) \log(\omega) \right ],
\end{equation}

\begin{equation}
\lambda = \sum_{\mathbf{q}\nu} \lambda_{\mathbf{q}\nu} =
2 \int \frac{\alpha^2F(\omega)}{\omega} d\omega.
\end{equation}

In the Allen-Dynes formalism, the Eliashberg function $\alpha^2F(\omega)$ is obtained by summing over all scattering processes at the Fermi level mediated by phonon momentum transfer and reads \cite{allen1975transition}:

\begin{equation}
\alpha^2F(\omega) = N(\epsilon_F) \frac{\sum_{\mathbf{k_1},\mathbf{k_2}}{
                    \left|
                    M_{\mathbf{k_1},\mathbf{k_2}}
                    \right|^2
                    \delta(\omega-\omega_{\mathbf{q}\nu}) \delta(\epsilon_{\mathbf{k_1}}) \delta(\epsilon_{\mathbf{k_2}})}}
                    {\sum_{\mathbf{k_1},\mathbf{k_2}}{
                     \delta(\epsilon_{\mathbf{k_1}}) \delta(\epsilon_{\mathbf{k_2}})}}.
\end{equation}

Here, $N(\epsilon_F)$ is the DOS at Fermi level, $\omega_{\mathbf{q}\nu}$ is the phonon spectrum of a branch $\nu$ at momentum $\mathbf{q}=\mathbf{k_2}-\mathbf{k_1}$, $\epsilon_{\mathbf{k_1}}$ and $\epsilon_{\mathbf{k_2}}$ are electronic band energies referred to the Fermi level, while $M_{\mathbf{k_1},\mathbf{k_2}}$ are the electron-phonon coupling matrix elements. Many-body effects introduce a change of spectral character at the Fermi level, where electronic correlations induce a mass enhancement and introduce a finite lifetime, due to incoherence. In this spirit of the DMFT scissors,
we correct the DFT bands with the renormalised DMFT band picture:
\begin{equation}
\alpha^2F(\omega) = \mathcal{A}_{tot} \frac{\sum_{\mathbf{k_1},\mathbf{k_2}}{
                    \left|
                    M_{\mathbf{k_1},\mathbf{k_2}}
                    \right|^2
                    \delta(\omega-\omega_{\mathbf{q}\nu}) \mathcal{A}(\mathbf{k_1}) \mathcal{A}(\mathbf{k_2})
                    }}
                    {\sum_{\mathbf{k_1},\mathbf{k_2}}{
                     \mathcal{A}(\mathbf{k_1}) \mathcal{A}(\mathbf{k_2})}},
\end{equation}
where $\mathcal{A}_{tot}$ and $\mathcal{A}(\mathbf{k})$ are respectively the total and $\mathbf{k}$-momentum resolved spectral weights at Fermi level. This approach is
denoted as \textit{DMFT+a2F} in the workflow.

Within the DFT+DMFT quantum embedding approach, the DFT Kohn-Sham eigenstates are used in the calculation of the DMFT Green function\cite{weber2021unifying,plekhanov2021computational}. We use atomic projectors to define the Anderson Impurity Model (AIM), which is successively solved within the Hubbard-I approximation.
A breadth of quantum solvers is readily available in the TRIQS open-source platform\cite{parcollet2015triqs, AICHHORN2016200}. In the full charge self-consistent approach (DFT+DMFT+CSC), the Kohn-Sham potentials are calculated from the DMFT electronic density, obtained using the DMFT occupancies. Upon DMFT convergence, total energies and forces are calculated using the Green function and self-energy. All DFT calculations in this work are carried out using the pseudopotential formalism.

\section{Conclusion}

We created a methodology for estimating the superconducting temperature in lanthanide hydrides with many-body corrections, and investigated a novel type of stable high-temperature superconducting material, CeH$_{16}$. The DMFT charge self-consistency, which involves many-body adjustments to the local charge density in first-principles calculations, is used to restore a consistent theoretical framework. A change in the spectral weight of the $f$ states causes a rise in predicted superconducting temperature, which influences the spectral character at the Fermi level. We discussed the capabilities for relaxing lanthanide hydrides within the DMFT formalism, built on our recent developments providing DMFT forces for underlying ultra-soft and norm-conserving pseudopotentials, despite the fact that many-body corrections have so far been limited to the electronic contributions to the Eliashberg function.
The latter gives structural insights, and we find that many-body effects on the lattice have negligible impacts at high pressures, since DFT structures and pressures are similar to their DMFT counterparts. Although a comprehensive treatment of phonons at the DMFT level is out of reach for such complicated materials, the latter shows that many-body adjustments to the electronic component are responsible for the strong corrections of $T_c$.
We investigated the aliovalent effect and discovered that compared to iso-structural La hydride and Ce hydride, an increase in $f$ character at the Fermi level leads to an increase in superconducting temperature, which is a compelling observation for future explorations of $f$ systems as high $T_c$ superconductors. Our method is general and provides a modular framework for interoperating common first-principles software, such as the freely accessible CASTEP+DMFT and Quantum Espresso codes, with a tiny numerical footprint and ease of implementation.

\section*{Data Availability}
The codes are available at url \textbf{dmft.ai} under the GPL 3.0 license.
\quad

\section*{Acknowledgements}
YW are supported from the China Scholarship Council, CW, NB and EP are supported by the grant [EP/R02992X/1] from the UK Engineering and Physical Sciences Research Council (EPSRC).
This work was performed using resources provided by the ARCHER UK National Supercomputing Service and the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk).

\quad

\section*{Additional Information}
Correspondence should be addressed to Evgeny Plekhanov (evgeny.plekhanov@kcl.ac.uk), Nicola Bonini (nicola.bonini@kcl.ac.uk) and Cedric Weber (cedric.weber@kcl.ac.uk).

\bibliographystyle{apsrev}


\section{Introduction} \label{sec:intro}
Super-massive black holes (SMBH) co-evolve with their host galaxies. This is evidenced by the strong correlations between the black hole mass and global properties (bulge mass, velocity dispersion and luminosity) of the host galaxy \citep{2009ApJ...698..198G,2011MNRAS.413.1479S,2013ARA&A..51..511K}.
Radiation, jets and outflows from an active galactic nucleus (AGN) can impose an enormous amount of energy into the environment of the SMBH on various physical scales, suppressing cooling flows and quenching star formation. Both theoretical and observational studies have suggested that AGN feedback is the primary driver for the co-evolution between SMBHs and their host galaxies (see reviews by \citealp{2012ARA&A..50..455F} and \citealp{2013ARA&A..51..511K}).

There exist two modes of black hole accretion, known as the cold mode and the hot mode \citep{2014ARA&A..52..529Y}, which correspond to two respective modes of AGN feedback: radiative mode and kinetic mode.
The radiative mode (or quasar mode) operates mainly in luminous AGNs whose  accretion rates are close to the Eddington limit. Radiatively driven ultra-fast outflows (UFOs) have been detected in the X-ray spectrum of up to 35\% of local luminous, radio-quiet AGNs, which are characterized by Fe XXV/XXVI absorption lines with blueshifted line-of-sight velocities of $\gtrsim10^4\rm~km~s^{-1}$ \citep{2010A&A...521A..57T}.
However, in the local universe luminous AGNs are found in $\lesssim10\%$ of normal galaxies \citep{2013ApJ...768....1M},
while the majority of present-day SMBHs grow their mass at low rates ($\lesssim$1\% Eddington limit) via a radiatively inefficient, hot accretion flow \citep{2014ARA&A..52..529Y}, manifesting themselves as low-luminosity active galactic nuclei (LLAGNs; \citealp{2008ARA&A..46..475H}).
The kinetic mode (or jet mode) operates in these LLAGNs, which are typically radio-loud sources due to relativistic jets symbiotic with the hot accretion flow.

Theories \citep{1999MNRAS.303L...1B} and numerical simulations \citep{2012MNRAS.426.3241N,2012ApJ...761..130Y,2012ApJ...761..129Y,2015ApJ...804..101Y,2021ApJ...914..131Y} matured over the past two decades predict that an energetic outflow in the form of an uncollimated hot wind must also be generated from the hot accretion flow. The wind has a mass outflow rate much larger than the accretion rate at the black hole horizon, thus playing a crucial role in the accretion process. Moreover, by depositing its energy and momentum to the environment, the hot wind can provide an efficient means of AGN kinetic feedback, in addition to the conventional jet-driven feedback. In fact, numerical simulations involving both the jet and the hot wind have shown that the momentum flux of the wind is always larger than that of jet \citep{2021ApJ...914..131Y}.
The role of the hot wind in star formation and black hole growth has been investigated in detail by \citet{2019ApJ...885...16Y}.
In the influential cosmological simulations of galaxy formation and evolution, IllustrisTNG \citep{2017MNRAS.465.3291W, 2018MNRAS.473.4077P}, winds from weakly accreting SMBHs are invoked to quench star formation in intermediate- to high-mass galaxies.

However, direct observational evidence of the putative hot wind remains scarce.
Using {\it Chandra} non-dispersed spectrum, \cite{2013Sci...341..981W} showed that the accretion flow onto the Galactic center SMBH, Sgr A*, is consistent with an inward decreasing accretion rate, inferring the existence of an outflow, but no direct information about the outflow was yielded.
\citet{2016Natur.533..504C} inferred the existence of nuclear winds in a sample of quiescent galaxies hosting LLAGNs, based on optical spectroscopic observations that reveal kpc-scale bipolar emission features co-aligned with
strong velocity gradients of ionized gas.
Most recently, \citet{2021NatAs...5..928S} reported evidence for a high-velocity ($\sim$3000 km~s$^{-1}$), high-temperature ($\sim$ 10 keV) outflow from M81*, a nearby prototypical LLAGN, based on a pair of symmetrically blueshifted and redshifted Fe XXVI Ly$\alpha$ emission lines in the high-resolution {\it Chandra} grating spectrum of this bright X-ray source.
This outflow was interpreted as the long-sought hot wind produced in the hot accretion flow onto M81$^*$, which propagates out to $\gtrsim10^{6}$ times the gravitational radius of the SMBH.

The finding of \citet{2021NatAs...5..928S} has begun to bridge a gap between observations of LLAGNs and the successful theory of hot accretion flows, as well as modern cosmological simulations of galaxy formation that require a new kinetic feedback mode from weakly accreting SMBHs, i.e., in the form of a wind.
However, M81* remains the only LLAGN in which compelling evidence of a hot wind is found.
In this work, we are thus motivated to search for signatures of the putative hot wind in more LLAGNs.
X-ray spectra of high spectral and angular resolutions afforded by {\it Chandra} have proved to be the key of revealing the hot wind. We have examined existing {\it Chandra} grating observations of local LLAGNs (see details in Appendix \ref{sec:sample}). Only one LLAGN, which resides in NGC\,7213, stands out as a promising candidate for a detailed analysis presented below.

Lying at a distance of 22.75 Mpc ($z$=0.00584, from the NASA/IPAC Extragalactic Database\footnote{http://ned.ipac.caltech.edu}),
NGC 7213 is classified as a Sa galaxy with both Syfert and LINER (low-ionization nuclear emission line region) signatures, presumably powered by a central SMBH with a mass of $\sim10^8\rm~M_{\odot}$ estimated from stellar velocity dispersion \citep{2002ApJ...579..530W}.
Using optical integral-field spectroscopic observations of the inner galaxy,
\cite{2014MNRAS.438.3322S} estimated a black hole mass of $8_{-6}^{+16}\times10^7\rm~M_{\odot}$ from the $M_{\rm BH}-\sigma$ relation of \cite{2009ApJ...698..198G}. Hereafter we adopt $M_{\rm BH}=10^8\rm~M_{\odot}$ as the fiducial black hole mass.
The bolometric luminosity of the LLAGN is estimated to be $9\times10^{42}\rm~erg~s^{-1}$ by \cite{2005MNRAS.356..727S} and $1.7\times10^{43}\rm~erg~s^{-1}$ by \cite{2012MNRAS.424.1327E} based on broadband spectral energy distribution fitting. These yield an average Eddington ratio of  $\lambda_{\rm Edd}=L_{\rm bol}/L_{\rm Edd}\sim10^{-3}$ ($L_{\rm Edd} = 1.3\times10^{38}[M_{\rm BH}/{\rm M_\odot}]\rm~erg~s^{-1}$ is the Eddington limit), which ensures the LLAGN classification.

The LLAGN in NGC\,7213 has been extensively studied in the X-ray band.
Spectra obtained by XMM-Newton and {\it Chandra} observations exhibit a number of prominent emission lines, in particular the K$\alpha$ line from neutral or weakly-ionized Fe and  highly-ionized iron lines including Fe XXV K$\alpha$ and Fe XXVI Ly$\alpha$ bwteen 6--7 keV \citep{2003A&A...407L..21B,2005MNRAS.356..727S,2008MNRAS.389L..52B,2013MNRAS.429.3439E}.
In the {\it Chandra} high-resolution grating spectrum, the Fe XXVI Ly$\alpha$ and XXV K$\alpha$ lines appear blueshifted at a substantial velocity of $\sim900\rm~km~s^{-1}$ and were suggested to originate from a collisionally ionized gas \citep{2008MNRAS.389L..52B}.
We shall revisit these lines in detail below to provide a physical understanding of their origin.
The underlying continuum shows no significant Compton reflection component \citep{2003A&A...407L..21B,2010MNRAS.408..551L,2015MNRAS.452.3266U}, as constrained by {\it BeppoSAX}, {\it Suzaku}, {\it Swift} and {\it NuSTAR} observations that cover the X-ray spectrum to above tens of keV.
Together with the narrow neutral Fe K$\alpha$ line, this suggests that
the standard optically-thick accretion disk may be truncated at a radius of $10^3-10^4~r_{\rm g}$ \citep{2010MNRAS.408..551L}, where $r_{\rm g}=GM_{\rm BH}/c^2$ is the gravitational radius of the black hole.
A hot accretion flow is likely responsible for the bulk of the X-ray emission from this LLAGN \citep{2010MNRAS.408..551L,2016MNRAS.463.2287X}, although contribution from a relativistic jet is also possible and indirectly suggested by the presence of a compact radio nuclear source \citep{2005MNRAS.356..727S}, which exhibits weakly correlated long-term variability with the X-ray nucleus \citep{2011MNRAS.411..402B}.
The four-decade-long X-ray light curve of the nucleus has also been associated with a fast-rise-exponential-decay behavior \citep{2018MNRAS.475.1190Y}, which was argued to be due to a delayed tidal disruption event.

The remainder of this paper is organized as follows. In Section \ref{sec:data}, we describe the relevant X-ray observations and our procedures of data reduction. Our X-ray spectral analysis, including a blind search of emission lines and characterization of the underlying thermal plasma, are presented in Section \ref{sec:result}. In Section \ref{sec:wind}, the observed X-ray spectrum exhibiting blueshifted lines from highly-ionized iron, is confronted with magnetohydrodynamical simulations to test the hot wind scenario. In Section \ref{sec:discuss},
we show that alternative origins of the observed Fe lines, such as stellar activity, AGN photoionization and hot accretion flow, are quantitatively disfavored, leaving the hot wind as the most plausible scenario.
We also compare the cases of NGC\,7213 and M81* and address general implications for LLAGN feedback.
In Section \ref{sec:sum} we summarize our results.
Quoted errors throughout this work are at 90\% confidence level, unless otherwise stated.

\section{Observations and data preparation} \label{sec:data}
In this work, we primarily utilized two sets of X-ray data of NGC\,7213: i) {\it Chandra} spectra obtained by the High Energy Transmission Grating (HETG) in combination with the Advanced CCD Imaging Spectrometer (ACIS), which are able to resolve emission lines expected from the hot plasma around the SMBH, at a velocity resolution of 300--2000 km~s$^{-1}$ over an energy range of 0.5--8 keV, and ii) {\it NuSTAR} spectra to constrain the broadband continuum over 3--79 keV.

NGC\,7213 was observed by {\it Chandra}/HETG on August 6 and 9, 2007 (ObsID 7742, 8590), with an exposure of 115 and 35 ks, respectively.
The data was downloaded from the public archive and reprocessed following the standard procedure with CIAO v4.13 and calibration files CALDB v4.9.5.
The $\rm\pm1$-order grating spectra of both the high energy grating (HEG) and medium energy grating (MEG) were extracted using \emph{tgextract}, for each of the two observations.
For the source spectra,
we adopted the default cross-dispersion half-width of 2.4\arcsec (corresponding to 265 pc at the assumed distance of NGC\,7213), which ensures the same enclosed energy fraction ($\sim$97\%) for different wavelengths\footnote{https://cxc.cfa.harvard.edu/proposer/POG/}.
The background spectra were extracted from two adjacent regions above and below the source region each with the default full-width of 19.1\arcsec.
The spectral extraction regions are illustrated in Figure~\ref{fig:merge} displaying a 0.5--8 keV counts image combining the two ACIS/HETG observations.
Our visual examination of this image found no significant point sources that may contaminate the source and background spectra.
The source shows no significant intra-observation and inter-observation variability.
Hence we co-added the $\rm\pm1$-order spectra from the two observations, respectively for the HEG and MEG, which result in a total clean exposure time of 148 ks.
The background spectra were also co-added, which contribute to $\lesssim2\%$ of the source counts.
We grouped the spectra to have at least one count per bin, to preserve the maximally possible spectral resolution for the subsequent spectral analysis.

NGC 7213 was targeted by {\it NuSTAR} on October 5, 2014, for an effective exposure of 101.6 ks (ObsID 6001031002). The data was downloaded and reprocessed following the standard \emph{nupipeline} in the software package NuSTARDAS v1.9.2.
The spectra of NGC\,7213 were extracted for both focal plane modules A and B (FPMA and FPMB) with \emph{nuproducts}. The source region was a circle with a radius of 100\arcsec\ centered on the galaxy, while the background spectra were extracted from a concentric annulus with an inner radius of 105\arcsec\ and an outer radius of 175\arcsec\ (Figure~\ref{fig:merge}).
The spectra were similarly grouped to have at least one count per bin.
Six off-nucleus point sources are present within the {\it NuSTAR} spectral extraction region. We find that their collective 0.5--8 keV flux accounts for less than 1\% of the flux of the nuclear source. Thus these off-nucleus sources should have a negligible contribution to the {\it NuSTAR} spectra.
We note that the source light curve shows no significant intra-observation variability.

\begin{figure}[htb!]
\centering\includegraphics[width=0.8\textwidth]{NGC7213_zoom_211206.pdf}
\caption{A 0.5--8 keV ACIS/HETG image combining the two {\it Chandra} observations of NGC\,7213.
The roll angles of the two observations differ by 180$^\circ$, hence the HEG/MEG dispersed arms of the two observations are indistinguishable here.
The blue solid rectangles mark the spectral extraction region, with the adjacent dashed rectangles defining the background.
The red solid and dashed circles mark the source and background spectral extraction regions for the {\it NuSTAR} observation.
The insert is a zoomed-in view of the circumnuclear region. The zeroth-order spectrum of the off-nucleus emission is extracted from the solid annulus, while the corresponding background is defined by the adjacent dashed annulus.
\label{fig:merge}}
\end{figure}

\section{Spectral analysis} \label{sec:result}
\subsection{Broadband continuum} \label{subsec:cont}
Spectral analysis is carried out with Xspec v.12.9.1, adopting the Cash statistics (\emph{C}-stat) to determine the best-fit models.
To characterize the broadband X-ray continuum, we first consider a simple absorbed power-law model \emph{tbabs*powerlaw}, as we find no hint of excess foreground absorption in the spectra.
Hereafter the absorption column density is fixed at the Galactic foreground value, $N_{\rm H}=1.08\times10^{20}\rm~cm^{-2}$.
The co-added MEG spectrum over 0.5--5 keV, the co-added HEG spectrum over 1--8 keV,
and the FPMA/FPMB spectra over 3-79 keV are jointly fitted, excluding energies between 6--7.5 keV energy to avoid influence by the known Fe lines.
The photon-index $\Gamma$ is tied between HEG and MEG and between FPMA and FPMB, but allowed to be different between the {\it Chandra} and {\it NuSTAR} spectra.
The normalization is allowed to vary among the four spectra, to account for possible flux variation and difference in the enclosed-energy fraction.
The best-fit photon-index is $\Gamma=1.723\pm0.013$ for HEG/MEG and $\Gamma=1.838\pm0.012$ for FPMA/FPMB, which are consistent with previous studies \citep{2008MNRAS.389L..52B,2015MNRAS.452.3266U}.
The point-spread function (PSF) corrected, unabsorbed 2-10 keV flux from the {\it Chandra} spectra is $\sim40\%$ higher than that from the {\it NuSTAR} spectra, in agreement with the flux decreasing trend found by \citet{2018MNRAS.475.1190Y}.
The best-fit absorbed power-law models are plotted against the observed spectra in Figure \ref{fig:baseline}, and the fitted parameters are listed in Table \ref{tab:mod_fit}.

\begin{figure}[htb!]
\centering\includegraphics[width=0.7\textwidth]{baseline_4data_20211130.pdf}
\caption{Observed broadband X-ray spectra of NGC\,7213.
Blue: co-added {\it Chandra} MEG spectrum;
Green: co-added {\it Chandra} HEG spectrum;
Red: {\it NuSTAR} FPMA spectrum;
Black: {\it NuSTAR} FPMB spectrum.
Here the spectra are grouped to have a signal-to-noise ratio greater than 3 per bin for better visualization.
Error bars are of 1$\sigma$.
The spectra are fitted with an absorbed power-law.
The lower panel shows the residual-to-error ratio. The energy range of 6--7.5 keV, where significant Fe lines are present, is excluded from the spectral fit.
\label{fig:baseline}}
\end{figure}

We then try to replace the power-law model with a cut-off power-law (\emph{cutoffpl} in Xspec), in order to test whether there is significant spectral softening at higher energies. We fit the {\it NuSTAR} spectra alone to constrain the high-energy cutoff, which results in a 3\,$\sigma$ lower limit, $E_{\rm cut}>93\rm~keV$, beyond the {\it NuSTAR} spectral coverage.
This suggests no significant spectral softening and that a single power-law is sufficient to account for the broadband continuum, especially at energies below 8 keV.
To further examine the possible presence of a reflection component, we employ the \emph{pexmon} model in Xspec, which can self-consistently calculate the reflected continuum from neutral gas and the fluorescent lines of neutral iron.
In this model, the intrinsic continuum assumes a power-law with an exponential cutoff.
We find a best-fit $\Gamma = 1.786^{+0.013}_{-0.005}$, $E_{\rm cut} = 184^{+32}_{-60}\rm~keV$,
and a reflection fraction $R=\frac{\Omega}{2\pi}<0.06~(3\,\sigma)$. During the fit the inclination angle has been fixed at $45^{\circ}$ to reduce model degeneracy, and
we have verified that changing this parameter has little effect on the derived reflection fraction.
The lack of a reflection component, also found by \citet{2010MNRAS.408..551L} based on {\it Suzaku} observations, indicates a truncated accretion disk.

We conclude that the broadband X-ray spectra of the LLAGN in NGC\,7213 can be well characterized by a power-law model.

\begin{deluxetable*}{lcccccccccc}
\footnotesize
\tablecaption{Broadband spectral fit results\label{tab:mod_fit}}
\tablewidth{0pt}
\tablehead{
\colhead{Data and model} & \colhead{$C$/d.o.f} & \colhead{$\Gamma$} & \colhead{$N_{\rm pl}$} & \colhead{$kT_{\rm h}$} &
\colhead{$kT_{\rm l}$} &
\colhead{$v_{\rm h}$} &
\colhead{$v_{\rm l}$} &
\colhead{$N_{\rm h}$} &
\colhead{$N_{\rm l}$} & \colhead{$N_{\rm wind}$}
}
\colnumbers
\startdata
PL & & & & & & & & & & \\
\hline
HEG (1--8 keV) & 4367/ 4212 & \multirow{2}*{$1.723^{+0.013}_{-0.013}$} & $5.69^{+0.09}_{-0.09}$ & - & - & - & - & - & - & -\\
MEG (0.5--5 keV) & 3750/ 3881 & {                   } & $5.96^{+0.06}_{-0.06}$ & - & - & - & - & - & - & -\\
FPMA (3--79 keV) & 1665/ 1592 &  \multirow{2}*{$1.838^{+0.012}_{-0.012}$} & $4.76^{+0.12}_{-0.12}$ & - & - & - & - & - & - & -\\
FPMB (3--79 keV) & 1581/ 1570 & {            } & $4.80^{+0.12}_{-0.12}$ & - & - & - & - & - & - & -\\
\hline
PL+apec$_{\rm h}$ & & & & & & & & & & \\
\hline
HEG (5--8 keV) & 354/356 & 1.723 & $4.1^{+0.9}_{-1.1}$ & $16^{+6}_{-6}$ & - & $-1.1^{+0.3}_{-0.4}$ & - & $4.4^{+1.9}_{-1.8}$ & - & -  \\
\hline
PL+wind & & & & & & & & & &\\
\hline
HEG (5--8 keV) & 356/358 & 1.723 & $2.4_{-0.2}^{+0.2}$ & - & - & - & - & - & - & $0.2^{+0.1}_{-0.1}$ \\
\hline
PL+apec$_{\rm h}$+apec$_{\rm l}$ & & & & & & & & & & \\
\hline
HEG (1--8 keV) & 4287/4206 & \multirow{2}*{1.723} & \multirow{2}*{$5.1^{+0.2}_{-0.2}$}  & \multirow{2}*{$13^{+5}_{-2}$} & \multirow{2}*{$0.8^{+0.2}_{-0.1}$} &  \multirow{2}*{$-1.2^{+0.1}_{-0.2}$} & \multirow{2}*{$0.04^{+0.12}_{-0.03}$} & \multirow{2}*{$1.9^{+0.6}_{-0.7}$} & \multirow{2}*{$0.05^{+0.03}_{-0.01}$} & - \\
MEG (0.5--5 keV) & 3740/3875 & { } & { } & { } & { } & { } & { } & { } & { } & -  \\
\enddata
\tablecomments{(1) The fitted data set and the adopted model; `h' and `l' denote high- and low-temperature, respectively. `PL' stands for power-law and `wind' denotes a synthetic wind model; (2) Cash statistics over degree of freedom; (3) Power-law photon-index; (4) Normalization of the power-law, in units of $10^{-3}\rm~ph~s^{-1}~cm^{-2}~keV^{-1}$ at 1 keV; (5)-(6) Plasma temperature of high- and low-temperature \emph{apec} models; (7)-(8) Line-of-sight velocity of high- and low-temperature \emph{apec} model, in units of $10^3\rm~km~s^{-1}$; (9)-(10) Normalization of high- and low-temperature \emph{apec} model, in units of $10^{-3}\rm~cm^{-5}$; (11) Normalization of the synthetic wind model, in units of $10^{-3}\rm~cm^{-5}$.
Errors are quoted at 90\% confidence level.
}
\end{deluxetable*}

\subsection{Blind search of lines} \label{subsec:blsearch}
Now we perform a blind search of emission/absorption lines in the co-added HEG and MEG spectra, following the method described in \citet{2010A&A...521A..57T} and \citet{2021NatAs...5..928S}.
Specifically, the putative line is modelled by a Gaussian profile, whose centroid energy $E_0$ and normalization $N$ are free to vary, and the dispersion is fixed at half of the HETG resolution (0.0056 $\rm \AA$ for HEG and 0.0111 $\rm \AA$ for MEG).
The Gaussian is added to the baseline power-law model derived in Section \ref{subsec:cont} (i.e., with a fixed photon-index of 1.723), progressively increasing $E_0$ channel by channel over the energy range of interest and seeking for the best-fit $N$.
The improvement of \emph{C}-stat, $\Delta C$, with respect to the baseline model is evaluated at each step.
For each pair of ($E_0$, $N$), confidence level of the line can be calculated according to ${\Delta}C$, given the fact that $C$-stat is asymptotically distributed as $\chi^{2}$.
The resultant contour plots are shown in Figure \ref{fig:bl_search} for the HEG spectrum between 5--8 keV.

\begin{figure}[hb!]
\centering\includegraphics[width=0.75\textwidth]{line_search_20211207.pdf}
\caption{Significant lines between 5--8 keV in the co-added HEG spectrum from a blind search.
The photon energy has been corrected for the systemic redshift of NGC\,7213.
Contours denote the improvement of $C$-stat ($\Delta C$) for a Gaussian line when added to the baseline continuum model. The blue, green, red contours represent values of $\Delta C$ equivalent to a confidence level of 99\%, 95\% and 68\%, respectively.
Three emission lines are identified at a confidence level $>95\%$.
\label{fig:bl_search}}
\end{figure}

Three emission lines, all with a confidence level $>95\%$, are identified from the HEG spectrum between 5--8 keV. We note that these lines were also identified by \citet{2008MNRAS.389L..52B} based on the same dataset.
The centroid energy of the three lines, after correcting for the systemic redshift ($z$=0.00584) of NGC\,7213, is 6.391 keV, 6.709 keV and 7.001 keV, respectively.
We have checked the robustness of line detection by varying the photon-index of the baseline power-law model between 1.70--1.74.
The resultant confidence levels of the three lines all remain above 95\%.
To refine the characterization of each line, we again apply the Gaussian model, this time setting the line width as a free parameter. The Gaussian fit results, including line centroid, line width, flux and equivalent width, are listed in Table \ref{tab:blsearch}.

\begin{deluxetable}{CLCLLL}
\tablecaption{Spectral fit results of the Fe lines\label{tab:blsearch}}
\tablewidth{0pt}
\tablehead{
\colhead{$E_{\rm 0}$} & \colhead{Significance} & \colhead{$E_{c\rm }$} & \colhead{$\sigma_{\rm E}$} & \colhead{$F$} & \colhead{EW}
}
\colnumbers
\startdata
6.391 & $>99.99\%$ & $6.395^{+0.007}_{-0.008}$ & $20^{+10}_{-8}$ & $3.1^{+0.8}_{-0.7}$ & $126.9^{+0.5}_{-0.8}$ \\
6.709 & $96.7\%$ & $6.717^{+0.016}_{-0.013}$ & $<55$ & $0.7^{+0.5}_{-0.5}$ & $26.0^{+0.4}_{-0.2}$ \\
7.001 & $>99.99\%$ & $6.996^{+0.012}_{-0.011}$ & $<99$ & $1.4^{+0.7}_{-0.6}$ & $61.2^{+0.8}_{-0.3}$ \\
\enddata
\tablecomments{(1) Most probable line centroid energy, in units of keV, determined by the blind search and corrected for the systemic redshift of NGC\,7213; (2) Significance level of the line; (3) Line centroid in units of keV, derived from the Gaussian fit; (4) Line width in units of eV; (5) Line flux in units of $10^{-13}\rm~erg~s^{-1}~cm^{-2}$; (6) Equivalent width in units of eV. Error bars are given at 90\% confidence level and upper limits of the line width are given at 3\,$\sigma$. }
\end{deluxetable}

The 6.391 keV line can be identified as the K$\alpha$ line of neutral or weakly ionized Fe.
This line is marginally resolved, with a best-fit line width of $20^{+10}_{-8}$ eV, corresponding to a velocity dispersion of $1.0^{+0.4}_{-0.4}\times10^3\rm~km~s^{-1}$.
The relatively low velocity dispersion and lack of significant Doppler shift suggest that this line arises from a truncated accretion disk or remote cold gas illuminated by the central LLAGN. A similar conclusion was drawn by \citet{2010MNRAS.408..551L}.
The 6.709 keV line is the weakest among the three lines and is consistent with the Fe XXV (He-like Fe) K$\alpha$ triplet, which is not fully resolved by HEG.
The 7.001 keV line may be
a blueshifted Fe XXVI (H-like Fe) Ly$\alpha$ line, as previously suggested by \citet{2008MNRAS.389L..52B}.
For the rest-frame energy of Fe XXVI Ly$\alpha$\footnote{Rest-frame line centroid energies throughout this work are taken from ATOMDB (http://www.atomdb.org).}, $6.966\rm~keV$, this implies for a line-of-sight velocity of $\sim -1.3\times10^3\rm~km~s^{-1}$.
We further consider the possibility that this line is Fe I K$\beta$, which has a rest-frame energy of 7.058 keV. However, this would imply for a line-of-sight velocity of $\sim2400\rm~km~s^{-1}$, which is inconsistent with the observed centroid of Fe I K$\alpha$. Moreover, the flux ratio between Fe I K$\beta$ and K$\alpha$ should have a canonical value of $\sim0.13$ \citep{2003A&A...410..359P}, while the observed ratio between the 7.001 keV and 6.391 keV lines is substantially larger ($0.4\pm0.2$). Hence we conclude that the 7.001 keV line cannot be associated with Fe I K$\beta$.
Guided by the implied blueshift of the Fe XXVI Ly$\alpha$ line, we refit the blended Fe XXV K$\alpha$ using a double Gaussian profile (Figure~\ref{fig:gauss_fit}a), obtaining line centroids at $6.660_{-0.019}^{+0.015}\rm~keV$ and $6.722_{-0.016}^{+0.017}\rm~keV$.
If the two Gaussians could be respectively associated with the forbidden transition (rest-frame energy at 6.637 keV) and resonant transition (rest-frame energy at 6.700 keV), they would both be blueshifted with a line-of-sight velocity of $\sim -1.0\times10^3\rm~km~s^{-1}$, nicely compatible with the putative blueshifted Fe XXVI Ly$\alpha$. This lends further support to the identification of the latter.

\begin{figure*}[htb!]
\centering\includegraphics[width=0.49\textwidth]{fit_4gauss_bin1_20211207.pdf}
\centering\includegraphics[width=0.49\textwidth]{fit_apec_20211207.pdf}
\caption{(a) The co-added HEG spectrum (black data points) fitted with a baseline power-law (magenta dotted line) plus four Gaussian lines (black dotted curves) representing, from low to high energies, Fe I K$\alpha$, Fe XXV K$\alpha$ (forbidden transition), Fe XXV K$\alpha$ (resonant transition) and Fe XXVI Ly$\alpha$. The total model is represented by the magenta solid curve.  The spectrum is binned to have a signal-to-noise ratio greater than 3 for better illustration. The error bars are of 1\,$\sigma$. The lower panel gives the residual-to-error ratio. (b) Similar to panel (a), but the fitted model consists of the baseline power-law, a blueshifted \emph{apec} (blue dotted curve) and two Gaussian lines representing Fe I K$\alpha$ and K$\beta$.
\label{fig:gauss_fit}}
\end{figure*}

\begin{deluxetable}{lcccccc}[hb!]
\tablecaption{Additional significant lines\label{tab:softband}}
\tablewidth{0pt}
\tablehead{
\colhead{Line identification} & \colhead{$E_{\rm lab}$} & \colhead{Significance} & \colhead{$E_{\rm c}$} & \colhead{$\sigma_{\rm E}$} & \colhead{$F$} & \colhead{EW}\\
\colhead{(1)} & \colhead{(2)} &
\colhead{(3)} & \colhead{(4)} & \colhead{(5)} & \colhead{(6)} & \colhead{(7)}
}
\startdata
\multirow{2}*{Mg XII Ly$\alpha$} & \multirow{2}*{1.473} & 97.2\% & $1.471^{+0.006}_{-0.006}$ & $<71$ & $0.7^{+0.4}_{-0.4}$ & $1.1^{+0.1}_{-0.1}$ \\
{               } & {   } & - & 1.471 & - & $<1.2$ & $0.97^{+0.01}_{-0.01}$ \\
\hline
\multirow{2}*{Ar XVIII Ly$\alpha$} & \multirow{2}*{3.321} & 96.9\% & $3.325^{+0.019}_{-0.010}$ & $<40$ & $2.3^{+1.8}_{-1.6}$ & $5.83^{+0.04}_{-0.10}$\\
{                 } & {   } & - & 3.325 & - & $<5.5$ & $2.6^{+0.1}_{-0.1}$ \\
\hline
\multirow{2}*{Ca XIX K$\alpha$ (r)} & \multirow{2}*{3.902} & - & 3.911 & - & $<5.5$ & $4.4^{+0.2}_{-0.2}$ \\
{              } & {   } & 96.1\% & $3.911^{+0.004}_{-0.003}$ & $<11$ & $2.5^{+1.8}_{-1.6}$ & $7.3^{+0.2}_{-0.1}$\\
\enddata
\tablecomments{(1)-(2) Identified line transition and rest-frame centroid energy in units of keV. Each line has two rows, the upper row for the MEG spectrum and the lower row for the HEG; (3) Significance of the line in the blind search; (4) Line centroid in units of keV, derived from the Gaussian fit; (5) Line width in units of eV; (6) Line flux in units of $10^{-14}\rm~erg~s^{-1}~cm^{-2}$; (7) Equivalent width in units of eV. Errors are quoted at 90\% confidence level, while all the upper limits are of 3\,$\sigma$.
}
\end{deluxetable}
\begin{figure}[ht!]
\centering\includegraphics[width=1\textwidth]{blind_search_1-5_20211209.pdf}
\caption{Blind search of lines in the co-added MEG (upper panel) and HEG (lower panel) spectra between 1--4 keV. Contours denote the improvement of $C$-stat (${\Delta}C$) for a Gaussian line when
added to the baseline continuum model. The blue, green, red contours represent values of ${\Delta}C$ equivalent to a confidence level of 99\%, 95\% and 68\%, respectively. An emission line is considered real if its significance were greater than 95\% in at least one of the two spectra and their flux or 3\,$\sigma$ upper limit were consistent between the two spectra. Three significant lines are thus detected and labeled with the tentatively identified atomic transition.
\label{fig:blsearch_soft}}
\end{figure}

Our blind search also finds several significant emission lines between 1--4 keV in the MEG and HEG spectra (Figure~\ref{fig:blsearch_soft}), where the two gratings have a comparable line sensitivity.
For these lines to be considered {\it real}, we have required that they are detected at a significance above 95\% in at least one of the two spectra and that their flux or 3\,$\sigma$ upper limit be consistent between the two spectra.
The basic properties of these lines are again derived from a Gaussian fit and listed in Table \ref{tab:softband}.
From the line centroid, we identify Ly$\alpha$ lines of Mg XII and Ar XVIII. Both lines are unresolved and have a Doppler shift  $\lesssim450\rm~km~s^{-1}$.
We also tentatively associate an unresolved 3.911 keV line with the resonant transition of Ca XIX K$\alpha$ (rest-frame energy 3.902 keV), which would imply for a Doppler shift $\lesssim700\rm~km~s^{-1}$.
We note that no significant emission lines are found between 4--6 keV. Nor do we find any significant absorption line in the MEG and HEG spectra.

\subsection{Modelling the hot plasma} \label{subsec:spec}
The blueshifted Fe XXVI Ly$\alpha$ and Fe XXV K$\alpha$ lines indicate a hot plasma with a substantial bulk motion. We use a phenomenological model to characterize such a plasma. Specifically, an optically-thin plasma model, {\it apec} in Xspec, is applied to fit the 5--8 keV HEG spectrum, with the plasma temperature and redshift as free parameters and the abundance fixed at solar.
The portion of the spectrum at lower energies is temporarily neglected, which may require an additional component as suggested by the emission lines of the low-$Z$ elements and a small excess to the power-law continuum seen at energies $\lesssim$1 keV in the MEG spectrum (Figure~\ref{fig:baseline}).
Foreground absorption is again fixed at the Galactic value.
Two Gaussian lines are added to account for the Fe I K$\beta$ and K$\alpha$ lines, with a fixed line ratio of 0.13 between the two.
The baseline power-law continuum is also included, with the photon-index fixed at 1.723.
This composite model results in a reasonable fit (Figure~\ref{fig:gauss_fit}b), with a plasma temperature $kT = 16\pm6\rm~keV$, which is mainly driven by the flux ratio between the Fe XXVI Ly$\alpha$ and Fe XXV K$\alpha$ lines.
After correcting for the systemic redshift of NGC\,7213, the fitted model requires a line-of-sight velocity of $-1.1^{+0.3}_{-0.4}\times10^3\rm~km~s^{-1}$, consistent with the expectation from the line centroids (Section~\ref{subsec:blsearch}).
The unabsorbed 2--10 keV of this 16-keV plasma is found to be $\sim4.0\times10^{41}\rm~erg~s^{-1}$.

While the HEG spectrum shows no sign of a redshifted counterpart to the blueshifted Fe lines, it is tempting to provide constraint on a potential redshifted component.
To do so, we add to the above composite model a second \emph{apec} component, requiring it to have the same temperature and a line-of-sight velocity exactly opposite to the blueshifted component.
The upper limit of this redshifted \emph{apec} is then determined by deviating the $C$-stat until the 3$\sigma$ level is arrived at. The 2--10 keV flux ratio between the redshifted and blueshifted components is thus constrained to be  $<$0.45.

The 16-keV plasma, however, cannot be responsible for the detected K$\alpha$ and Ly$\alpha$ emission lines of the low-$Z$ elements (Table~\ref{tab:softband}), since at such a high temperature these low-$Z$ elements would be fully ionized.
To account for these lines,
we again add a second {\it apec} component to the above composite model, this time allowing both the plasma temperature and redshift to vary freely.
It turns out that this new composite model provides a reasonable simultaneous fit to the 1--8 keV HEG spectrum and 0.5--5 keV MEG spectrum.
Compared to the single power-law model, the composite model improves the $C$-stat by 80 for the HEG spectrum and 10 for the MEG spectrum for 6 less degree-of-freedom.
The best-fit temperature of the second {\it apec} is $0.8_{-0.1}^{+0.2}\rm~keV$,
and the best-fit redshift is nearly zero, suggesting little bulk motion.
The unabsorbed 0.5-10 keV luminosity of this component is $8.7\times10^{39}\rm~erg~s^{-1}$.
This signifies a low-$T$ plasma that might be associated with a circumnuclear hot gas.
It is noteworthy that with the addition of the low-$T$ component, the best-fit temperature of the blueshifted high-$T$ component becomes $13^{+5}_{-2}$ keV, still consistent with the value (16 keV) determined from fitting the 5--8 keV spectrum alone.

The spectral fit results in this section are summarized in Table \ref{tab:mod_fit}.

\subsection{Zeroth-order spectrum of the off-nucleus region} \label{subsec:circum_reg}
We have also examined the X-ray spectrum of an off-nucleus region (Figure~\ref{fig:offcenter}), which is extracted from the zeroth-order ACIS/HETG image and provides meaningful constraint on the spatial extent of the thermal components found in the first-order dispersed spectra.
The off-nucleus region is defined as an annulus with an inner-to-outer radius of 2.5\arcsec--5\arcsec\ (275--550 pc; see the insert in Figure~\ref{fig:merge}), which is immediately outside the first-order spectral extraction region
and also avoids potential effect of pileup caused by the bright nucleus.
The corresponding background spectrum is extracted from a concentric annulus with an inner-to-outer radius of 5\arcsec--8\arcsec.
Spectra extracted from the two {\it Chandra} observations are co-added.
We find that a power-law model subjected to the Galactic foreground absorption can well describe this spectrum.
The best-fit photon index of $1.7\pm0.1$ and an unabsorbed 0.5--10 keV luminosity of $2.5\times10^{40}\rm~erg~s^{-1}$
can be understood as the collective emission from unresolved X-ray binaries plus PSF-scattered photons from the bright nucleus.

While a thermal component is not formally required by the apparently featureless spectrum, we add an \emph{apec} model to the fit, fixing the abundance at solar and the plasma temperature at $0.8$ keV as found for the low-$T$ component in the first-order spectra (Section~\ref{subsec:spec}).
The maximum normalization of this \emph{apec} is found such that the deviation of $C$-stat arrives at the 3\,$\sigma$ level, which is indicated by the blue dotted curve in Figure \ref{fig:offcenter}.
The corresponding unabsorbed 0.5--10 keV luminosity is $\sim4\times10^{38}\rm~erg~s^{-1}$, which is $\sim$20 times lower than that of the 0.8 keV component detected in the first-order spectrum.
In terms of surface brightness, the 0.8 keV plasma is a factor of 70 times lower in this off-nucleus region, considering the underlying projected area of zeroth- and first-order spectra. This indicates that the 0.8 keV plasma is concentrated within the central $\sim$300 pc and strengthens our above suggestion that it is tracing a circumnuclear hot gas.

\begin{figure}[ht!]
\centering\includegraphics[width=0.6\textwidth]{fit_offcenter_20211207.pdf}
\caption{Co-added zeroth-order spectrum of the off-nuclear region in NGC\,7213.
The spectrum is binned to have a signal-to-noise ratio greater than 3 for better illustration.
The spectrum is fitted by an absorbed power-law (black dashed line) plus a maximally-allowed \emph{apec} model (blue dotted curve). The black solid curve marks the total model. The lower panel shows the residual-to-error ratio. Error bars are of 1\,$\sigma$.
\label{fig:offcenter}}
\end{figure}

\section{A hot wind from the LLAGN as the origin of the Fe lines}
\label{sec:wind}
The X-ray spectra presented in the previous section reveal the presence of a high-velocity, hot plasma in the LLAGN of NGC\,7213.
This is highly similar to the case of M81$^*$, in which a high-velocity, hot plasma is traced by highly-ionized Fe lines and interpreted as an energetic wind driven by the hot accretion flow onto the central SMBH \citep{2021NatAs...5..928S}.
Here we consider this same scenario for NGC\,7213, which is illustrated in Figure \ref{fig:AD_absorb}.
In this schematic diagram of an inflow-outflow system around a weakly accreting SMBH, a cold thin disk dictates the accretion inflow at large radii, but becomes truncated at a radius $r_{\rm tr} \sim 10^3~r_{\rm g}$. Inside $r_{\rm tr}$, the disk is replaced by a geometrically-thick hot inflow.
Close to the event horizon, a highly-collimated, relativistic jet is expected to launch.
Also symbiotic with the hot accretion flow is a bipolar, non-relativistic wind, which ultimately carries away the bulk of the inflowing mass.
Compared to the jet, the wind originates from a wide range of radius in the accretion flow, having a much wider opening angle and much smaller radial velocity.
Due to adiabatic expansion, the wind has its temperature decreasing from that of the hot accretion flow and spanning a range well suited for producing the highly-ionized Fe lines.
We are thus motivated to perform magnetohydrodynamical (MHD) simulations of the putative wind, tailored to the physical conditions of the SMBH in NGC\,7213, to test this favorable origin of the observed Fe lines.

\begin{figure}[ht!]
\centering\includegraphics[width=0.4\textwidth]{sketch_211209.pdf}
\caption{Schematic diagram of the inflow-outflow system around the central SMBH in NGC\,7213 (represented by the black circle). The pair of grey triangles represent the geometrically-thick hot accretion flow. The black stripes along the equatorial plane mark the truncated thin disk, while the vertical yellow stripe denotes a relativistic jet. The blue and red cones with large opening angles represent the hot wind propagating towards and away from the observer, where emission lines of highly-ionized Fe are generated. In the case of NGC 7213, a high viewing angle (defined as away from the jet axis) is inferred, such that a large fraction of the redshifted hot wind is obscured by the truncated disk.
\label{fig:AD_absorb}}
\end{figure}

\subsection{Wind simulation and synthetic wind spectrum}
\label{subsec:simulation}
We perform a 2.5-dimension MHD simulation with the ZEUS-MP/2 code version 1.5.19 \citep{2006ApJS..165..188H} in spherical coordinates.
The black hole mass is set to be $10^8\rm~M_{\odot}$ based on observations \citep{2002ApJ...579..530W}.
The simulation setup is similar to  \citet{2021NatAs...5..928S}, which is briefly described here.
The wind starts from the inner boundary, which is set at 2000 $r_{\rm g}$, i.e., the assumed disk truncated radius, appropriate for the moderate Eddington ratio of the LLAGN in NGC\,7213 and also consistent with the narrowness of the Fe I K$\alpha$ line (Section~\ref{subsec:blsearch}).
We ignore a possible wind from the cold disk outside the truncated radius, which is expected to have too low temperatures to produce the highly ionized Fe.
The rotational energy, kinetic energy, thermal energy and magnetic energy of the wind are taken into account.
Following \citet{2015ApJ...804..101Y} and \citet{2021ApJ...914..131Y}, we set the rotation velocity as $0.8v_{\rm k}$, wind velocity as $0.6v_{\rm k}$, temperature as $0.5kT_{\rm vir}$, and the plasma $\beta=p_{\rm gas}/p_{\rm mag}=0.2$ at the inner boundary, where $v_{\rm k}=(GM_{\rm BH}/r)^{0.5}$, $kT_{\rm \rm vir}=GM_{\rm BH}/(3r)$, $p_{\rm gas}$, and $p_{\rm mag}$ are the Keplerian velocity, virial temperature, gas pressure and magnetic pressure, respectively.
Only the toroidal magnetic field is taken into account, since the poloidal field is typically much smaller than the toroidal component at large radii as found in wind simulations \citep{2021ApJ...914..131Y}. The grid of the simulation is extended from $2\times10^3~r_{\rm g}$ to $10^6~r_{\rm g}$ in radial direction. A radially injecting boundary condition is set for the inner boundary, while an outflow boundary condition is adopted at the outer boundary.
The grid in latitudinal direction is extended from $\theta = 10^{\circ}$ to $50^{\circ}$ to avoid the singularity of spherical coordinate and to neglect the jet. The existence of a relativistic jet does not affect the dynamics of the wind, nor does it produce emission lines.
A reflective boundary condition is set for the latitudinal direction, since we expect that the pressure from the truncated disk will confine the hot wind.
\begin{figure*}[ht!]
\centering\includegraphics[width=0.49\textwidth]{temp_distribution_20211209.pdf}
\centering\includegraphics[width=0.485\textwidth]{density_distribution_20211209.pdf}
\caption{Plasma temperature (a) and density (b) distribution of the hot wind from the MHD simulation. The radial velocity field is overlaid with arrows. The y-axis follows the jet axis and the x-axis defines the equatorial plane of the truncated disk. Regions occupied by the jet or the disk have been excluded.
\label{fig:distri}}
\end{figure*}

A direct outcome of the simulation is the two-dimensional (2D) distributions of density, temperature and velocity, as shown in Figure~\ref{fig:distri}.
We note that the absolute scale of the density is determined by matching the synthetic wind spectrum to the observed spectrum (see Section~\ref{subsec:synthetic}).
The synthetic wind spectrum is produced based on the simulated  quantities.
Specifically, we first expand the 2D grid into a 3D grid, given the axisymmetry about the jet axis and mirror symmetry with respect to the equatorial plane.
We then calculate the density-weighted X-ray spectrum of each grid in the black hole rest frame, utilizing ATOMDB for a given metal abundance to include atomic transitions and bremsstrahlung.
The synthetic spectrum is then calculated by integrating along a given viewing angle, i.e., the angle between the jet axis and line-of-sight.
In this step we take into account the effect of gravitational lensing, which will boost the redshifted component, and the effect of relativistic Doppler shift, which will boost (reduce) the blueshifted (redshifted) component.
Self-absorption in the hot wind is neglected, since the hot wind is expected to be optically-thin, which we verify in Appendix~\ref{sec:depth}. However, we have considered potential absorption by the truncated thin disk, which could be optically-thick to X-ray photons from behind the disk, for certain light-of-sight (see illustration in Figure~\ref{fig:AD_absorb}).

\subsection{Comparison with the observed HEG spectrum}
\label{subsec:synthetic}
The synthetic wind spectrum is compared with the observed HEG spectrum over 5--8 keV.
Added to the wind model are again the baseline power-law and two Gaussian lines representing Fe I K$\alpha$ and K$\beta$.
We emphasize that this comparison is not a formal fit to the observed spectrum. Rather, our aim is to examine whether the synthetic spectrum of the hot wind, naturally arising from the standard paradigm of hot accretion flow onto a weakly accreting SMBH, can provide a satisfactory description to the observed spectrum, in particular the blueshifted Fe emission lines.

By construction of our wind simulation and synthetic spectrum, the primary free parameter is the viewing angle. We find that a viewing angle of 80$^{\circ}$ results in a synthetic spectrum that best matches the observed spectrum (Figure \ref{fig:wind}).
In this case, the bulk of the highly-ionized Fe lines are produced in the radial range of $10^4-10^5~r_{\rm g}$, where the plasma temperature ranges between $3\times10^7-3\times10^8$ K and the radial velocity ranges between 1500--8000$\rm~km~s^{-1}$ (Figure~\ref{fig:distri}).

\begin{figure}[hb!]
\centering\includegraphics[width=0.6\textwidth]{fit_wind_bin1_20211206.pdf}
\caption{The synthetic wind spectrum (green dotted curve) for a viewing angle of 80$^\circ$ and a Fe abundance of twice solar is confronted  with the observed spectrum. Other symbols are the same as in Figure~\ref{fig:gauss_fit}.
\label{fig:wind}}
\end{figure}

The large viewing angle has two main effects. First, the projected line-of-sight velocity of the Fe line prevalent region takes values of $\sim1100\rm~km~s^{-1}$, compatible with the observed Doppler shift.
Second, the redshifted component, which arises from behind the equatorial plane, is largely obscured by the truncated disk.
This breaks the intrinsic symmetry between the blueshifted and redshifted components as a result of the intrinsic symmetry of the bipolar wind about the equatorial plane by design, and provides a natural explanation to an unseen redshifted counterpart to the blueshifted Fe XXVI Ly$\alpha$ line in the observed spectrum (Section~\ref{subsec:spec}).
Quantitatively, the thin disk has a large hydrogen column density of $4.3\times10^{26}\left(\frac{\alpha}{0.1}\right)^{-4/5}\left(\frac{M_{\rm BH}}{10^8\rm~M_{\odot}}\right)^{1/5}\left(\frac{\dot{m}}{0.001}\right)^{7/10}\left(\frac{r}{2000~r_{\rm g}}\right)^{-3/4}\rm~cm^{-2}$ along the direction perpendicular to the disk plane, for the black hole mass $M_{\rm BH}$, the accretion rate normalized by the Eddington limit $\dot{m}$, and viscosity parameter $\alpha$ suitable for NGC\,7213 \citep{2008bhad.book.....K}.
The corresponding effective optical depth of electron scattering,  $\tau_{*}\sim4.7\times10^3(\frac{\alpha}{0.1})^{-4/5}(\frac{M_{\rm BH}}{10^8\rm~M_{\odot}})^{1/5}(\frac{\dot{m}}{0.001})^{1/5}$, is sufficiently high to block the X-ray photons from behind the disk.

Alternatively, the weak or absent redshifted component may be due to an intrinsic asymmetry of the hot wind.
As discussed in \citet{2021NatAs...5..928S}, the gas density of the wind above and below the equatorial plane can deviate from each other at a level of $\sim25$\% on the dynamical timescale ($\sim$10 yr) of the wind, which can lead to a $\sim50$\% difference between the observed redshifted and blueshifted components.
While such an explanation was offered for the moderate inequality between the observed redshifted and blueshifted Fe XXVI Ly$\alpha$ lines in M81*, the substantially weaker redshifted component in the present case of NGC\,7213 disfavors an intrinsic variability as the primary cause.

A second tunable parameter in the synthetic spectrum is the Fe abundance. Initially, we have assumed a solar abundance as adopted in the \emph{apec} models (Section~\ref{subsec:spec}).
However, the equivalent width of the Fe XXVI Ly$\alpha$ line derived from the synthetic spectrum is 32 eV, which is only about half of the observed value (61.2 eV).
Hence we adopt a synthetic spectrum with the Fe abundance being twice solar, which leads to an equivalent width of 45 eV for Fe XXVI Ly$\alpha$ with the $C$-stat improved by ${\Delta}C =7$.
We emphasize that the abundance cannot be strongly constrained, due to a degeneracy in the continuum between the wind  bremsstrahlung and contribution from the hot accretion flow (and/or a jet, if present) as represented by the power-law.

Given the best-matching synthetic wind spectrum, the outflow rate of the hot wind can be evaluated according to the simulated  density and radial velocity distributions.
We find an outflow rate of $\sim0.08\rm~M_{\odot}~yr^{-1}$.
Interestingly, this value is comparable to the inflow rate of $0.07-0.2\rm~M_{\odot}~yr^{-1}$ estimated for warm ionized gas within the central $\sim100$ pc traced by H$\alpha$ emission \citep{2014MNRAS.438.3322S}.
Theoretically, the hot wind is expected to carry away the bulk of the inflowing mass in the hot accretion flow \citep{2015ApJ...804..101Y}. The good agreement between the inferred outflow and inflow rates lends strong support to the wind scenario.

\section{Discussion} \label{sec:discuss}

\subsection{Alternative origins of the high-velocity hot plasma} \label{subsec:alter}
\subsubsection{Stellar activities} \label{subsec:stellar}
In principle, the high-velocity Fe lines may be produced by stellar activities such as young supernova remnants (SNRs) or massive star binaries with strong colliding winds.
However, the X-ray luminosity of the 16-keV plasma, $4\times10^{41}\rm~erg~s^{-1}$, is orders of magnitude higher than possible for a single young SNR \citep{2012MNRAS.419.1515D}.
Theoretical investigation \citep{2016MNRAS.456.2537R} also shows that the cumulative X-ray luminosity from a collection of SNRs evolving within the sphere of influence of a $\sim$10$^8\rm~M_{\odot}$ SMBH ($\sim$100 pc in radius) does not exceed $\sim$10$^{38}\rm~erg~s^{-1}$.
Similarly, colliding wind binaries cannot reach the observed X-ray luminosity and plasma temperature \citep{2012ASPC..465..301G}.

To further evaluate the collective X-ray emission from the circumnuclear stellar populations, we first estimate the stellar mass and star formation rate (SFR) in the HETG spectral extraction region.
The stellar mass within a projected radius of 265 pc is found to be $\sim1\times10^{10}\rm~M_{\odot}$, based on the $V$-band surface brightness distribution of NGC\,7213 \citep{2005AJ....129.2138L} and assuming a $V$-band mass-to-light ratio of 5.4 typical of Sa galaxy \citep{1979ARA&A..17..135F}.
The SFR is estimated to be $\lesssim0.02\rm~M_{\odot}~yr^{-1}$, based on the measured luminosity of narrow H$\alpha$ line from the nucleus of NGC\,7213 \citep{1984ApJ...285..458F} and the empirical relation between SFR and H$\alpha$ luminosity \citep{1998ARA&A..36..189K}.
These estimates then provide an upper limit in the X-ray luminosity contributed by both old and young stellar populations, $L_{\rm 2-10~keV}^{\rm star} \lesssim 9.4\times10^{38}\rm~erg~s^{-1}$, according to the empirical relation given by \citet{2010ApJ...724..559L}: $L_{\rm 2-10~keV}^{\rm star} = \alpha M_* + \beta {\rm SFR}$, where $\alpha = (9.05\pm0.37)\times10^{28}{\rm~erg~s^{-1}~M_\odot^{-1}}$ and $\beta = (1.62\pm0.22)\times10^{39}{\rm~erg~s^{-1}~(M_\odot~yr^{-1})}^{-1}$.
This is orders of magnitude lower than the observed X-ray luminosity of the 16-keV plasma.
Moreover, the circumnuclear star formation is expected to produce a 0.5--2 keV luminosity of  $L_{0.5-2}^{\rm SFR} = 9\times10^{37}({\rm SFR}/0.02{\rm~M_{\odot}~yr^{-1}})\rm~erg~s^{-1}$, based on the empirical relation of \citet{2003A&A...399...39R}. This is about two orders of magnitude lower than the observed X-ray luminosity of the low-$T$ component (Section \ref{subsec:spec}).
Hence we conclude that stellar activities cannot account for the thermal components in the HETG spectrum.

\subsubsection{AGN Photoionization} \label{subsec:phoion}
Alternatively, the highly ionized Fe lines may be produced by AGN photoionization.
We evaluate this possibility using the photoionization code {\it Cloudy} (ver c17.02; \citealp{2017RMxAA..53..385F}). The code calculates the ionization state of an isotropic and uniform gas cloud photoionized by a central source, which has an intrinsic X-ray spectrum, $\phi(\nu)$, same as the baseline power-law continuum (Section~\ref{subsec:cont}), i.e., with a photon-index of 1.723 and 0.5--10 keV luminosity of $2.4\times10^{42}\rm~erg~s^{-1}$. We note that the X-ray luminosity has varied by a factor of $\sim4$ over a period of $\sim$38 yr \citep{2018MNRAS.475.1190Y}, but our conclusion below is not affected by this moderate variability.
The inner radius of the cloud is set to be 0.001 pc ($\sim260~r_{\rm g}$), while the outer radius varies with different assumed values of column density $N_{\rm c} = n_{\rm c}(r_{\rm out}-r_{\rm in})$, where $n_{\rm c}$ is the cloud density assumed to be constant for simplicity.
We note that adopting a radially-decreasing density distribution does not significantly affect our following conclusion.
The ionization parameter at the illuminated face of the cloud is expressed by $U_{\rm X}=\int^{\rm 10keV}_{\rm 2keV}\frac{\phi(\nu)d\nu}{4\pi r_{\rm in}^2cn_{\rm c}}$.
The lack of reflection component in the broadband X-ray spectrum indicates a Compton-thin case.
Hence we explore the plausible parameter space of $N_{\rm c} =10^{21}-10^{24}\rm~cm^{-2}$ and ionization parameters $U_{\rm X}
=1-7$.
For a given set of $N_{\rm c}$ and $U_{\rm X}$,
the outward luminosity of a certain atomic line is then calculated using the theoretical line emissivity from ATOMDB and assuming an abundance of twice solar, to be consistent with the favored abundance in the synthetic wind spectrum.

The predicted line luminosities are shown in Figure~\ref{fig:phoion}.
We find that the predicted Fe XXVI Ly$\alpha$ line from a photoionized cloud falls short by a factor of a few to a few hundred to the observed luminosity, for any reasonable combination of $N_{\rm c}$ and $U_{\rm X}$.
The predicted Fe XXV K$\alpha$ line is also typically weak and only comes closer to the observed luminosity at small values of $U_{\rm X}$.
Therefore, the hot plasma traced by the highly-ionized Fe lines are highly unlikely due to photoionization by the LLAGN. We note that \citet{2008MNRAS.389L..52B} also disfavored an origin of photoionization for the Fe lines based on the argument that the resonant transition dominates the Fe XXV triplet.
Line luminosities are also predicted for the three low-$Z$ lines detected: Mg XII Ly$\alpha$, Ar XVIII Ly$\alpha$ and Ca XIX K$\alpha$.
Among them, the Mg XII Ly$\alpha$ line can have a predicted luminosity comparable to the observed value, whereas the other two predicted lines typically fall short to the observed value except with the lowest $U_{\rm X}$. This suggests that photoionization by the LLAGN can at best account for a fraction of the low-$Z$ lines.
\begin{figure*}[ht!]
\centering\includegraphics[width=0.49\textwidth]{phoion_N23A2_211229_Fe.pdf}
\centering\includegraphics[width=0.49\textwidth]{phoion_N24A2_211229_Fe.pdf}
\\
\centering\includegraphics[width=0.49\textwidth]{phoion_N23A2_211229_other.pdf}
\centering\includegraphics[width=0.49\textwidth]{phoion_N24A2_211229_other.pdf}
\caption{The predicted line luminosity induced by AGN photoionization, as a function of the ionization parameter. The upper panels are for the Fe XXV K$\alpha$ and Fe XXVI Ly$\alpha$ lines, and the lower panels are for the Mg XII Ly$\alpha$, Ar XVIII Ly$\alpha$ and Ca XIX K$\alpha$ lines. The left and right panels assume a gas column density of $10^{23}\rm~cm^{-2}$ and $10^{24}\rm~cm^{-2}$, respectively.
The horizontal dotted lines and the adjacent stripes mark the observed luminosity and 90\% uncertainty.
\label{fig:phoion}}
\end{figure*}

\subsubsection{The hot accretion flow} \label{subsec:inflow}
In principle, the hot accretion inflow inside the truncation radius can also produce highly-ionized Fe lines with a substantial Doppler shift.
To test such a possibility, we utilize a set of 3D general relativity (GR) MHD simulation of the hot accretion flow with Kerr metric performed with the Atheta++ code \citep{2016ApJS..225...22W}.
The simulation is scale-free in black hole mass, hence we can scale the simulation data to $M_{\rm BH}=10^8\rm~M_{\odot}$. Although the 3D GRMHD simulation only extends from the event horizon to 100 $r_{\rm g}$, all physical quantities of the hot accretion flow, such as density, temperature and velocity have a tight scaling relation in the radial direction, which allow us to extrapolate these quantities from 100 $r_{\rm g}$ to 2000 $r_{\rm g}$.
An inflow rate of $0.08\rm~M_{\odot}~yr^{-1}$ at $2000~r_{\rm g}$ is adopted to normalize the gas density, which is consistent with the derived outflow rate in the hot wind (Section~\ref{subsec:compare}).

We then follow the same procedure for the wind simulation to produce synthetic spectrum of the hot accretion flow for a given line-of-sight.
The Fe abundance is set to be twice solar.
We find that within $1000~r_{\rm g}$ the plasma temperature reaches  $>8\times10^{9}\rm~K$, which is so high that the gas becomes fully ionized; the bulk of Fe XXV and Fe XXVI ions exist in the radial range of $1000-2000~r_{\rm g}$, where the temperature drops to $\sim10^9$ K.
However, the predicted equivalent width of the Fe XXVI Ly$\alpha$ line is only $\sim15\rm~eV$, much low than the observed value of 61.2 eV.
Moreover,
due to the wide spread of the projected velocity in this radial region, the resultant line profile takes a broad and flat-topped shape, which is inconsistent with the observed narrow line profile.
We also find it difficult to produce a large flux ratio as observed between the blueshifted and redshifted components, since the inflow itself is optically-thin and is barely obscured by the truncated disk.
These discrepancies between the synthetic and observed spectra persist for a number of viewing angles tested.
Hence we conclude that the blueshifted highly-ionized Fe lines are not originated from the hot accretion flow.

\subsection{Comparison with M81* and implications for LLAGN feedback} \label{subsec:compare}
In the previous sections, we have shown that the blueshifted lines from highly ionized Fe are best understood as arising from a hot wind lunched from the hot accretion flow onto the SMBH in NGC\,7213.
A mass outflow rate of $\sim0.08\rm~M_{\odot}~yr^{-1}$ is inferred for this hot wind (Section~\ref{subsec:compare}).
We further estimate an associated kinetic power of $\sim3\times10^{42}\rm~erg~s^{-1}$.
This is equivalent to $\sim15$\% of the SMBH's bolometric luminosity \citep{2018MNRAS.475.1190Y}, a fraction very similar to that found for the hot wind in M81* \citep{2021NatAs...5..928S}.
The momentum flux of the hot wind is estimated to be $\sim4\times10^{33}\rm~g~cm~s^{-1}$, which is $\sim$6 times the photon momentum flux from the LLAGN.
This is again similar to the case of M81*.
In view that M81* and NGC\,7213 encompass a substantial range in the Eddington ratio ($3\times10^{-5} -10^{-3}$), this similarity may suggest a quasi-universal mechanical output of the hot wind with respect to the SMBH's radiative output, which deserves further examination when detection of the hot wind in more LLAGNs becomes available.
In any case, the kinetic power and momentum flux of the hot wind can provide an efficient and long-lasting means of feedback from weakly accreting SMBHs that are prevalent in normal galaxies.

Among the crucial open questions about LLAGN feedback is whether the hot wind can deposit its kinetic energy and momentum into a substantial volume of the host galaxy, as required to effectively quench star formation in cosmological simulations  \citep{2017MNRAS.465.3291W,2018MNRAS.473.4077P}.
Numerical simulations predict that the hot wind will interact with the circumnuclear medium, releasing its kinetic energy and momentum and potentially creating strong shocks \citep{2019ApJ...885...16Y}.
The 0.8-keV plasma found in the HETG spectrum might be evidence for such an interaction, which is likely the manifestation of a circumnuclear hot gas shock-heated by the outward propagating hot wind.
We note that circumnuclear hot gas with a similar temperature is also evident in the case of M81* \citep{2021NatAs...5..928S}.
The internal energy of this circumnuclear hot gas in NGC\,7213 is $\lesssim1\times10^{54}\rm~erg$, assuming that it fills a spherical volume with a radius of 265 pc (i.e., the HETG spectral extraction region). This translates to a reasonably short time of $\sim1\times10^{4}\rm~yr$ if the entire kinetic power of the hot wind were converted into this internal energy.

The hot wind could have propagated to larger distances from the LLAGN.
The zeroth-order spectrum of the off-nucleus region examined in Section~\ref{subsec:circum_reg} allows for the presence of a 0.8-keV plasma out to $\sim550$ pc, but with a substantially lower surface brightness.
This may be understood as either the hot wind having been impeded and had the bulk of its kinetic power exhausted by the circumnuclear medium, or a paucity of gas in the off-nucleus region.
Interestingly, recent ALMA CO observations reveal the presence of non-rotating molecular gas features at distances of few hundred pc from the nucleus \citep{2020A&A...641A.151S}, which were interpreted as a molecular gas outflow with a velocity of $\lesssim100$ km~s$^{-1}$ and a mass outflow rate of $0.03\pm0.02\rm~M_\odot~yr^{-1}$.
While \citet{2020A&A...641A.151S} favored stellar activity as the mechanism responsible for the outflow,
the possibility that this molecular gas is accelerated and entrained by the hot wind deserves further investigation.

At present, direct evidence for the hot wind from a weakly accreting SMBH is only found in M81* and NGC\,7213. This may be attributed to two limiting factors. On the one hand, high-resolution X-ray spectra are currently only available for a handful of nearby LLAGNs (Appendix~\ref{sec:sample}). Among these sources, M81* and NGC\,7213 stand out with better data quality to allow for the detection of highly-ionized Fe lines.
On the other hand, even with a deep X-ray spectrum, the presence of a strong X-ray continuum from  a relativistic jet (e.g., in the case of Cen A and M87) may easily mask the emission lines.
Future high-resolution X-ray spectroscopic observations by {\it Chandra}, as well as planned missions such as XRISM and Athena, hold promise to confirm the universality of a hot wind from LLAGNs.

\section{Summary}
\label{sec:sum}
We have searched the {\it Chandra} archive for ACIS/HETG observations of LLAGNs that may reveal potential evidence for a putative hot wind driven by a hot accretion flow and characterized by highly-ionized iron lines, as recently established in the case of M81*.
Only one source, the LLAGN in NGC\,7213, manifests itself with X-ray spectral features similar to those found in M81*.
We have performed an in-depth X-ray spectral analysis for NGC\,7213, as well as MHD simulations of the hot wind tailored to the empirical conditions of the weakly accreting SMBH in this galaxy.
Our main findings are as follows.

\begin{itemize}
    \item[(i)] A blind search results in the detection of a number of emission lines from the {\it Chandra} 1st-order HEG/MEG spectra. In particular, Fe XXVI Ly$\alpha$ and Fe XXV Ka$\alpha$ lines with a blue-shifted velocity of $\sim1100\rm~km~s^{-1}$ are identified. The flux ratio between the two lines suggest that they are tracing a hot plasma with a characteristic temperature of $16\pm6$ keV. The spectra show no significant redshifted counterpart to this plasma. In addition, several emission lines from low-$Z$ elements indicate the presence of a $\sim$0.8 keV plasma, which may be tracing a circumnuclear hot gas.
    \item[(ii)] The synthetic X-ray spectrum of the hot wind from the MHD simulation with a viewing angle of $\sim80^{\circ}$ and an Fe abundance twice solar can well match the blueshifted, highly-ionized Fe lines. In this case, the unseen redshifted counterpart is most likely due to obscuration by the truncated disk along the line-of-sight. All plausible alternative origins considered for the Fe lines, such as stellar activity, AGN photoionization and the hot accretion flow itself, are quantitatively disfavored. This makes a strong case that a hot wind exists in the  nucleus of NGC\,7213, presumably driven by the hot accretion flow.
    \item[(iii)] The mass outflow rate of the hot wind is estimated to be $\sim0.08\rm~M_{\odot}~yr^{-1}$, quite comparable to the inflow rate of ionized gas estimated for warm ionized gas within the central $\sim$100 pc traced by H$\alpha$ emission. This agreement between the outflow and inflow rates lends further support to the hot wind scenario for NGC\,7213.
     \item[(iv)] The wind carries a kinetic energy of $\sim3\times10^{42}\rm~erg~s^{-1}$ and a momentum flux of $\sim4\times10^{33}\rm~g~cm~s^{-1}$, thus providing an important means of kinetic feedback to the environment.
\end{itemize}

Future X-ray and optical spectroscopic observations in a spatially-resolved fashion will be crucial to determining to what extent and at what efficiency the hot wind can affect the interstellar medium of the host galaxy.

\begin{acknowledgments}
This research made use of observations taken from the {\it Chandra} X-ray Observatory, software package CIAO provided by the {\it Chandra} X-ray Center and spectral fitting package XSPEC,
and the High Performance Computing Resource in the Core Facility for Advanced Research Computing at Shanghai Astronomical Observatory.
Z.L. and F.S. acknowledges support by the National Key Research and Development Program of China (2017YFA0402703) and Natural Science Foundation of China (grant 11873028).
F.Y. and B.Z. are supported in part by the Natural Science Foundation of China (grants 11633006, 12133008, 12192220, and 12192223).
\end{acknowledgments}

\software{CIAO \citep{2006SPIE.6270E..1VF}, XSPEC \citep{1996ASPC..101...17A}, NuSTARDAS (v1.9.2), ZEUS-MP/2 (v1.5.19; \citealp{2006ApJS..165..188H}), Cloudy (vc17.02; \citealp{2017RMxAA..53..385F}), Atheta++ code \citep{2016ApJS..225...22W}}