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

Primordial black holes (PBHs) can form in the early Universe, prior to galaxies and stars, through a variety of mechanisms and can constitute a significant fraction of the dark matter (DM)~(e.g.~\cite{Zeldovich:1967,Hawking:1971ei,Carr:1974nx,Khlopov:1980mg,Khlopov:1985jw,GarciaBellido:1996qt,Khlopov:2008qy,Frampton:2010sw,Kawasaki:2016pql,Carr:2016drx,Inomata:2016rbd,Pi:2017gih,Inomata:2017okj,Garcia-Bellido:2017aan,Georg:2017mqk,Inomata:2017vxo,Kocsis:2017yty,Ando:2017veq,Cotner:2016cvr,Cotner:2017tir,Cotner:2018vug,Sasaki:2018dmp,Carr:2018rid,Banik:2018tyb,1939PCPS...35..405H,Cotner:2019ykd,Kusenko:2020pcg,Flores:2020drq,Domenech:2021uyx}). Depending on the formation mechanism, the PBH mass can span many orders of magnitude. Those PBHs with mass above $M \gtrsim 10^{15}$~g survive Hawking evaporation until present. In a sizable mass range, $M \sim 10^{-16} - 10^{-10} M_{\odot}$, PBHs can constitute all of the DM~\cite{Katz:2018zrn,Smyth:2019whb,Montero-Camacho:2019jte,Sugiyama:2019dgt}. PBHs have been associated with a variety of astrophysical phenomena and puzzles, e.g. PBHs could be the seeds for the formation of super-massive black holes~(e.g.~\cite{Bean:2002kx,Kawasaki:2012kn,Clesse:2015wea}) and could lead to new signatures~from compact star disruptions from PBH capture~(e.g.~\cite{Fuller:2017uyd,Takhistov:2017nmt,Takhistov:2017bpt,Takhistov:2020vxs}).

PBHs in the mass window of $\sim 10-10^3 M_\odot$ have been suggested as progenitors \cite{Nakamura:1997sm,Clesse:2015wea,Bird:2016dcv,Raidal:2017mfl,Eroshenko:2016hmn,Sasaki:2016jop,Clesse:2016ajp} for the gravitational wave (GW) signals detected by LIGO~(e.g. \cite{Abbott:2016blz,Abbott:2016nmj,Abbott:2017vtc}). The first definitive observation of intermediate-mass $150 M_\odot$ BH has been claimed by LIGO \citep{Abbott:2020tfl}. In conventional stellar evolution astrophysics such BHs are difficult to form from the collapse of a single star, akin to the formation of stellar-mass BHs. While this PBH DM parameter space region is already under pressure from a variety of constraints
(see Ref.~\citep{Carr:2020gox} for review), these rely on and are sensitive to uncertainties of underlying assumptions, which often include the cosmological history.

In Ref.~\cite{Lu:2020bmd} we have put forth a novel cosmology-independent observable, based on PBHs heating of the surrounding interstellar medium (ISM) gas, as a powerful probe for stellar and intermediate-mass PBHs contributing to the DM. In this work we further explore this observable and provide a detailed analysis of PBH gas heating as well as various underlying processes, such as PBH emission mechanisms. Following Ref.~\cite{Lu:2020bmd}, we focus on the Leo T irregular dwarf galaxy as a particularly favorable astrophysical system target. Further, we extend the analysis to include systems with bulk velocity relative to the DM.

The ``heating'' of near-by stars by passing massive BHs has been previously considered~\cite{Carr:1999dyn,Lacey:1985,Totani:2010}. Theoretically expected X-ray emission from traversing PBHs interacting with the surrounding medium has been also employed to set limits on the PBH DM fraction~\cite{Gaggero:2016dpq,Inoue:2017csr}. Constraints from gas heating and system cooling have been considered in the context of particle DM candidates~\cite{Bhoonah:2018wmw,Farrar:2019qrv,Wadekar:2019xnf}, which however have heating mechanisms different from PBHs. Further, while the focus for particle DM is on target astrophysical systems with high relative particle velocities resulting in high collision rates, for PBHs the astrophysical systems with slower relative velocities are favored due to increased PBH accretion and resulting emission.

While our discussion focuses on present-day effects, emission from gas-accreting PBHs is also relevant for cosmology, e.g. for the reionization and recombination epochs (e.g. \cite{Ricotti:2007au}). X-ray emission from gas accreting onto PBHs would result in ionization and increased quantities of molecular hydrogen, which would in turn provide cooling and enhance earlier star formation~\cite{Ricotti:2007au}. Further, emission from accreting PBHs could result in spectral $\mu-$ and $y-$distortions of the cosmic microwave background (CMB).
Heating and ionization of intergalactic medium gas from PBH emission could modify 21cm power spectrum in an observable way~\cite{Mena:2019nhm}.
In this work we do not discuss cosmological effects, which provide complementary ways of exploring PBH DM.

\section{PBHs in ISM}
\label{sec:accretion}

\subsection{Bondi-Hoyle-Lyttleton accretion}
\label{ssec:bondihoyle}

X-ray emission from gas-accreting isolated BHs traversing the ISM of the Milky Way was explored in~Ref.~\cite{Agol:2001hb}. Ref.~\cite{Inoue:2017csr} adapted the analysis for PBHs constituting part of the DM. In this section we broadly follow the same method. PBHs traveling at low Mach number $\mathcal{M} = v/c_s$ accrete via the Bondi-Hoyle-Lyttleton process, with an accretion rate~\cite{1939PCPS...35..405H,1944MNRAS.104..273B,1952MNRAS.112..195B}
\begin{align}
\label{eq:bondihoyle}
    \dot{M} = 4 \pi r_B^2 \tilde{v} \rho = \frac{4\pi G^2 M^2 n \mu m_p}{\Tilde{v}^3}~.
\end{align}
Here $M$ is the PBH mass, $r_B = G M/\tilde{v}^2$ is the Bondi radius\footnote{Bondi radius sets
the scale at which gravity becomes important and approximates how far from object material starts to be drawn in and accreted.}, $\mu$ is the mean molecular weight, $n$ is the ISM gas number density, $m_p$ is the proton mass and
 $\tilde{v} \equiv (v^2+c_s^2)^{1/2}$, where $v$ is the PBH velocity relative to the ISM gas and $c_s$ is the  temperature-dependent sound speed in gas, which we take to be approximately $c_s \sim 10$~km/s~\cite{Inoue:2017csr}. At high Mach number, the BH accretion rate is expected to be limited to $\sim10$--$20$\% of the canonical Bondi-Hoyle-Lyttleton accretion rate of Eq.~\eqref{eq:bondihoyle}~\cite{Guo:2020qkf}, as shown by recent 3D hydrodynamical simulations.
 However, the resulting system behavior is sensitive to assumed initial conditions. We find that it is the low-velocity regime that predominantly contributes to our effects of interest. Hence, we adopt the standard accretion rate of Eq.~\eqref{eq:bondihoyle}.

The bolometric luminosity of the accretion disk of a BH can be related to its accretion rate as $L=\epsilon(\dot{M})\dot{M}$. Here $\epsilon(\dot{M})$ is an effective phenomenological parameter describing radiative efficiency, which varies according to different accretion regimes.
Assuming a radiative efficiency of $\epsilon_0 = 0.1$ characteristic of thin disks\footnote{The radiative efficiency $\epsilon_0$ can vary from 0.057 for a non-rotating Schwarzschild BH to
0.42 for an extremal Kerr BH~(see e.g.~Ref.~\cite{Kato:book2008}).}, the Eddington accretion rate in terms of Eddington luminosity is
\begin{equation}
\label{eq:eddingtonrate}
    \dot{M}_\textrm{Edd} = \frac{L_\textrm{Edd}}{\epsilon_0 c^2} = 6.7\times10^{-16}\left(\frac{M}{M_\odot}\right)~M_\odot\text{/s}~.
\end{equation}
To determine the accretion regimes of interest, it is useful to introduce a dimensionless form of the accretion rate
\begin{equation}
\label{mdot}
    \dot{m} = \frac{\dot{M}}{\dot{M}_\textrm{Edd}} = 2.64\times10^{-7} \left(\frac{M}{M_\odot}\right)\left(\frac{n}{1\textrm{ cm}^{-3}}\right) \left(\dfrac{\mu}{1}\right) \left(\frac{\Tilde{v}}{10\textrm{ km}/\textrm{s}}\right)^{-3}~.
\end{equation}

\subsection{Accretion disk formation}

An accretion disk can form around a BH from the gas infalling with sufficient angular momentum.
For accretion rates $\dot{m}\gtrsim 0.07\alpha$, a standard thin $\alpha$-disk~\cite{Shakura:1972te} likely forms. Here the disk's viscosity is described by a phenomenological parameter $\alpha \sim \mathcal{O}(0.1)$.
Perturbations in the density or velocity distributions of the surrounding gas can supply angular momentum to the disk. Assuming a Kolmogorov spectrum of density fluctuations in the ISM, the resulting density differential on opposing sides of the accretion ``cylinder'' can be calculated~\cite{1976ApJ...204..555S}. We equate the angular momentum from this differential with the angular momentum of the outer edge of the thin disk~\cite{Agol:2001hb},
\begin{equation}
\label{eq:router}
    r_{\rm out} \simeq 5.4\times10^{9} r_s \left(\frac{M}{100 M_\odot}\right)^{\frac{2}{3}}\left(\frac{\Tilde{v}}{10 \textrm{ km/s}}\right)^{-\frac{10}{3}}~,
\end{equation}
where $r_s = 2 G M/c^2$ is the BH Schwarzschild radius. On the other hand, for inefficient accretion described by the Advection-Dominated Accretion Flow (ADAF) regime~\cite{Narayan:1994xi,Yuan:2014gma}, the inflowing matter does not follow a Keplerian description. In this case, the outer accretion disk radius is uncertain and we consider a range of parameter values as described in Sec.~\ref{sec:heating}.

The inner edge of the accretion disk is assumed to follow
the innermost stable circular orbit (ISCO) of a test particle, which for a Schwarzschild BH is $r_{\rm ISCO} = 3 r_s$. If the outer edge of the BH accretion disk, at $r_{\rm out}$, is found to be located within the ISCO radius, i.e. $r_{\rm out}/r_{\rm ISCO} < 1$, we do not expect efficient disk formation.
Thus, using Eq.~\eqref{eq:router} we can infer an approximate condition on the mass of a BH for accretion disk formation as
\begin{equation}
M > 2.5 \times 10^{-13} M_{\odot} \left(\frac{\Tilde{v}}{10 \textrm{ km/s}}\right)^{5}~.
\end{equation}
This condition is well-satisfied for all BHs in our parameter range of interest. On the other hand, it shows that small PBHs are generally not expected to have an accretion disk.

\subsection{Gas and PBH distributions}
\label{ssec:pbhdist}

The total number of PBHs of mass $M$ within a gas system of volume $V$ and DM density $\rho_{\rm DM}$, assuming PBHs make up a fraction $f_{\rm PBH}$ of the DM, is given by
\begin{equation}
N_{\rm PBH}(M) = f_{\rm PBH} \frac{\rho_{\rm DM}V}{M}~.
\end{equation}
Our considerations apply in case there is at least one PBH in the gas system. This results in a lower ``incredulity" limit
\begin{equation}
\label{eq:incredlimit}
    f_\textrm{bound} > \frac{3M}{4\pi r_\textrm{sys}^3 \rho_\textrm{DM}}~,
\end{equation}
where $r_\textrm{sys}$ is the size of the system\footnote{For simplicity, we assume spherical DM halo distribution as convention, however other distributions are also possible~(e.g.~\cite{Hayashi:2012si})}.

Here we assume a monochromatic PBH mass function, however our results can be readily applied to an extended PBH mass distribution. The velocity of PBHs contributing to the DM can be described by a Maxwell-Boltzmann distribution
\begin{equation}
\label{eq:maxwellian}
    f(v) = \sqrt{\frac{2}{\pi}}\frac{v^2}{\sigma_v^3}\exp\left(-\frac{v^2}{2\sigma_v^2}\right)~,
\end{equation}
where $\sigma_v$ is the velocity dispersion within a given system.

We are interested in the gas heating due to PBHs. For a gas system in thermal equilibrium, we estimate the total amount of heating $H_{\rm tot}$ from PBHs of mass $M$ as
\begin{equation}
\label{eq:totalheateq}
    H_{\rm tot}(M) = N_{\rm PBH}(M) H(M) \simeq \int_{n_\textrm{min}}^{n_\textrm{max}}\int_{v_\textrm{min}}^{v_\textrm{max}} dn dv \frac{df_n}{dn}\frac{df_v}{dv}\mathcal{H}(M,n,v)~,
\end{equation}
where $df_n/dn$ is the gas number density distribution, $df_v/dv$ is the relative PBH velocity distribution and $\mathcal{H}(M,n,v)$ is the amount of heat being deposited into the system from a single PBH and also takes into account other potential contributing effects (e.g. photon absorption). $\mathcal{H}$ represents the cumulative contribution of all three heating mechanisms we consider below. For photon emission as well as outflows, there is an additional integration to treat absorption. For the gas systems of interest the density is approximately constant, for which $df_n/dn$ can be taken as a delta function in $n$.

\section{Gas Heating Mechanisms}
\label{sec:heating}

We consider three distinct gas heating mechanisms due to PBHs: absorption of photons emitted from the accretion disk, direct heating due to dynamical friction and energy deposit by protons of outflows that are inherent to the ADAF regime. Of these three, dynamical friction has the least uncertainty but typically constitutes a subdominant contribution. Since for our relevant parameters of interest the PBHs are in the stellar and intermediate-mass BH-range, the heating associated with Hawking evaporation is negligible. Following the original proposal of Ref.~\cite{Lu:2020bmd}, PBH gas heating has been recently considered in the context of small evaporating PBHs as well~\cite{Kim:2020ngi,Laha:2020vhg}.

\subsection{Accretion photon emission}
\label{ssec:xrays}

Photon emission from PBH accretion disks is one of the major contributors to the ISM gas heating. Emission in the X-ray band is most readily absorbed by the surrounding hydrogen ISM gas, becoming efficient when the mass accretion rates are significant.

For different accretion disk regimes we follow Ref.~\cite{Yuan:2014gma}.  We consider that the accretion flow results in a thin disk formation when $\dot{m}>\dot{M}=0.07\alpha$, and is described by ADAF for less efficient accretion. We do not consider the ``slim disk'' and other solutions relevant for near-Eddington or super-Eddington accretion when $\dot{m}\sim 1$, as the PBHs in our parameter space of interest are not expected to typically achieve such high rates. Furthermore, the slim disk emission spectrum is not significantly different from the spectrum of the thin disk for $\dot{m}\lesssim 10$~\cite{1999ApJ...522..839W}. In Fig.~\ref{fig:emissionregions}, we display
the primary photon emission band from accreting PBHs over a wide range of PBH masses and for different densities of the surrounding gas, assuming Bondi-Hoyle-Lyttleton accretion of Eq.~\ref{eq:bondihoyle} with a characteristic flow velocity of $\tilde{v} \simeq 10$~km/s, and emission results of Sec.~\ref{sssec:thindisk} and Sec.~\ref{sssec:ADAF}. Fig.~\ref{fig:emissionregions} suggests which PBH observational signatures could be generically expected for a particular environment.

\begin{figure}[tb]
\begin{center}
\includegraphics[trim={5mm 0mm 40 0},clip,width=0.65\textwidth]{plots/emissionregions.pdf}

\caption{  \label{fig:emissionregions}
Dominant photon emission from PBH accretion over a wide mass-range of PBHs $M$, for different number densities $n$ of the surrounding gas, assuming Bondi-Hoyle-Lyttleton accretion and a characteristic flow velocity of $\tilde{v} \simeq 10$~km/s in Eq.~\eqref{eq:bondihoyle}. Regions of different accretion regimes $\dot{m}$ (black diagonal lines) resulting in slim disk, thin disk and ADAF (including LHAF, ADAF, eADAF sub-regimes, see text) accretion flows are shown.}
\end{center}
\end{figure}

Hydrogen gas absorbs ionizing radiation efficiently, but below the threshold of $E_i = 13.6~\textrm{eV}$, hydrogen is optically thin to continuous emission. For non-relativistic gas velocity dispersions, Doppler broadening of the emission spectra is not significant so that we can ignore scattering with photons that have energies lower than $E_i$. The majority of photons are absorbed and a significant fraction of the energy is deposited as heat when the medium is optically thick. For $13.6\textrm{ eV}<E<30\textrm{ eV}$, we employ the photo-ionization cross-section of Refs.~\cite{Bethe1957,1990A&A...237..267B},
\begin{equation}
\label{eq:ionizationformula}
    \sigma(E) = \sigma_0 y^{-\frac{3}{2}}\left(1+y^{\frac{1}{2}}\right)^{-4}~,
\end{equation}
where $y=E/E_0$, $E_0 = 1/2 E_i$ and $\sigma_0 = 605.73 \textrm{ Mb} = 6.06\times10^{-16}\textrm{ cm}^2$. The optical depth for a gas system of size $l$ and number density $n$ is then given by
\begin{equation}
    \tau (n, E) = \sigma(E) n l~.
\end{equation}
Above $30\textrm{ eV}$, we use attenuation length data from Fig.~(32.16) of Ref.~\cite{Olive_2014}.

The resulting heating power due to accretion disk emission photons is
\begin{equation}
\label{eq:heatingphotons}
    \mathcal{H}(M,n,v) = \int_{E_i}^{E_\textrm{max}}L_\nu (M,n,v) f_h \left(1-e^{-\tau}\right)d\nu~,
\end{equation}
where $L_{\nu}(\nu)$ is the photon luminosity at frequency $\nu$ and $f_h$ is the fraction of energy deposited as heat that we estimate to be $\sim 1/3$.
Both the ADAF and thin disk emission spectra decrease exponentially at high energies, so we evaluate the integral up to a maximum energy $E_{\rm max}$, above which the contributions are negligible\footnote{For thin disk, we choose the maximum temperature to be $E_{\rm max} = 5 T_i$, where $T_i$ is the characteristic temperature of the inner disk. For ADAF, we take $E_{\rm max} = 3 T_e$, with $T_e$ being the temperature of the electron plasma.}.

\subsubsection{Thin disks}
\label{sssec:thindisk}

The standard (geometrically) thin  $\alpha$-disk, characterized by the viscosity parameter $\alpha$, is optically thick and efficiently emits blackbody radiation~\cite{Shakura:1972te}. Thin disk allows for a fully analytic description. The temperature of the disk varies with radius as~\cite{Pringle:1981ds}
\begin{equation}
\label{eq:alphatemp}
    T(r) = T_i \left(\frac{r_i}{r}\right)^{3/4} \left[1-\left(\frac{r_i}{r}\right)^{\frac{1}{2}}\right]^{\frac{1}{4}}~,
\end{equation}
where $r_i$ is the inner disk radius taken to be the ISCO radius $r_{\rm ISCO} = 3 r_s$ as before. Here, using Eq.~\eqref{eq:bondihoyle},
\begin{equation} \label{eq:tialpha}
    T_i = \left(\frac{3GM\dot{M}}{8\pi r_i^3 \sigma_B}\right)^{1/4} = 53.3\textrm{ eV}\left(\frac{n}{1\textrm{ cm}^{-3}}\right)^{\frac{1}{4}} \left(\frac{\Tilde{v}}{10\textrm{ km}/\textrm{s}}\right)^{-\frac{3}{4}}~,
\end{equation}
where $\sigma_B$ is the Stefan-Boltzmann constant.
For Bondi-Hoyle-Lyttleton accretion, the disk temperature is seen to be independent of the BH mass. The disk temperature $T(r)$ decreases at large radii as $r^{-3/4}$, but also diminishes near the inner radius. The disk attains a maximum temperature of $T_{\rm max} = 0.488 T_i$ at a radius of $r = 1.36 r_i$.

The thin disk emission spectrum is a combined contribution of the blackbody spectra from each radius. Using the scaling relations in Ref.~\cite{Pringle:1981ds} and requiring continuity, the
resulting spectrum can be approximately described as
\begin{align}
\begin{split}
\label{eq:thinspectrum}
    \nu < T_o :~~~~~& L_\nu = c_\alpha \left(\frac{T_{\rm max}}{T_o}\right)^{\frac{5}{3}}\left(\frac{\nu}{T_\textrm{max}}\right)^2 \\
    T_o < \nu < T_{\rm max} :~~~~~& L_\nu = c_\alpha \left(\frac{\nu}{T_\textrm{max}}\right)^{\frac{1}{3}} \\
    T_{\rm max} < \nu :~~~~~& L_\nu = c_\alpha \left(\frac{\nu}{T_\textrm{max}}\right)^2 e^{1-\frac{\nu}{T_{\rm max}}}
\end{split}
\end{align}
where
\begin{equation}
c_\alpha = 1.27\times10^{29}\textrm{ erg}\textrm{ eV}^{-1}\textrm{s}^{-1} \Big(\dfrac{M}{M_\odot}\Big)^2 \Big(\dfrac{n}{1\textrm{ cm}^{-3}}\Big)^{3/4}\Big(\dfrac{\Tilde{v}}{10\textrm{ km}/\textrm{s}}\Big)^{-9/4}~.
\end{equation}
Here $c_\alpha$ is normalized so that the total integrated luminosity (power) is $\int L_{\nu} d\nu = 0.057\dot{M}c^2$, i.e.~emission has the maximum possible efficiency for a Schwarzschild black hole. Using Eq.~\eqref{eq:bondihoyle}, \eqref{eq:router} and \eqref{eq:alphatemp}, the temperature at the outer disk radius $T_o = T(r_{\rm out})$ is given by
\begin{equation}
\label{eq:touter}
    T_o \simeq 6.13\times10^{-5}\textrm{ eV} \left(\frac{n}{1\textrm{ cm}^{-3}}\right)^{\frac{1}{4}}\left(\frac{M}{M_\odot}\right)^{-\frac{1}{2}}\left(\frac{\Tilde{v}}{10\textrm{ km}/\textrm{s}}\right)^{\frac{7}{4}}~.
\end{equation}
For all relevant parameter values of interest for $n,M,\Tilde{v}$, the temperature is always $T_o \ll 13.6~\textrm{eV}$. Hence, we only need to consider the temperature regimes in which $T_o < \nu$.

Eq.~\eqref{eq:tialpha} shows that thin disks efficiently emit at X-ray frequencies, a regime where the hydrogen ionization cross section is large. Hence, the emission from thin disks would be well absorbed and could contribute significantly to gas heating.

\subsubsection{ADAF - emission power}
\label{sssec:ADAF}

When the accretion is significantly sub-Eddington, ADAF disks form~\cite{Narayan:1994xi,Yuan:2014gma}.
Here, the heat generated by viscosity during accretion is not efficiently radiated out and much of the energy is advected via matter heat capture into the BH event horizon, along with the gas inflow. In contrast to the thin disk, the ADAF disk is geometrically thick and optically thin.

An ADAF disk has a complex multi-component emission spectrum, of which we will consider the three dominant ones associated with electron cooling\footnote{We neglect additional possible contribution of synchrotron radiation from non-thermal electrons, which is present in ADAF models of Sgr A$^\ast$~\cite{2003ApJ...598..301Y}.}: synchrotron radiation, inverse Compton (IC) scattering and bremsstrahlung.

To describe the ADAF spectrum, we will employ approximate analytic expressions obtained in Ref.~\cite{Mahadevan:1996jf} in combination with the updated values for the phenomenological input parameters. The accretion inflow forms a hot two-temperature plasma~\cite{1995ApJ...452..710N}, with the hotter ions at $T_\textrm{ion}\simeq 10^{11} K$ and the cooler radiating electrons with $T_e \simeq 10^{9}-10^{10} K$. We consider the temperature to be spatially uniform over the entire emission region. To characterize the flow, we use the following parameter values consistent with recent numerical simulations and observations~\cite{Yuan:2014gma}: the fraction of the viscously dissipated energy that directly heats
electrons\footnote{The remainder $(1-\delta)$ goes into the ions.} $\delta = 0.3$, the ratio of gas pressure to total pressure $\beta = 10/11$, the minimum flow radius is the ISCO radius with $r_\textrm{min}=3r_s$, and the viscosity parameter $\alpha = 0.1$. For our results below, we employ $x_M \simeq 4\times 10^3 \dot{m}^{1/4}$ and constants $c_1=0.5$, $c_3=0.3$ in the formalism of Ref.~\cite{Mahadevan:1996jf} to characterize the ADAF solution~\cite{Narayan:1994xi}.

For a more robust analysis, we will further divide ADAF into three distinct sub-regimes, depending on the accretion rate, following Ref.~\cite{Yuan:2014gma}.  Namely, we shall consider luminous hot accretion flow (LHAF) for accretion rates $\dot{M}_\textrm{ADAF}=0.1\alpha^2 < \dot{m} < \dot{M}_\textrm{LHAF} =0.07\alpha $, standard ADAF in the range $10^{-3} \alpha^2=\dot{M}_\textrm{eADAF} < \dot{m} < \dot{M}_\textrm{ADAF}$, and ``electron" ADAF (eADAF) for $\dot{m} < \dot{M}_\textrm{eADAF}$. As discussed below in Section~\ref{sssec:ADAFtemp}, the temperature determination is handled differently for each of these ADAF regimes.

The synchrotron emission is self-absorbed and produces a rising $L_\nu \propto \nu^{5/2}$ spectrum, similar to the $\nu^{1/3}$ spectrum in the optically thin regime and also to that of the thin disk, which peaks at
\begin{equation}
\label{eq:nup}
    \nu_p =1.83\times10^{-2}\textrm{ eV} \left(\frac{\alpha}{0.1}\right)^{-\frac{1}{2}}\left(\frac{1-\beta}{1/11}\right)^{\frac{1}{2}}\theta_e^2\left(\frac{r_\textrm{min}}{3 r_s}\right)^{-\frac{5}{4}}\left(\frac{M}{M_\odot}\right)^{-\frac{1}{2}}\left(\frac{\dot{m}}{10^{-8}}\right)^\frac{3}{4}~,
\end{equation}
where
\begin{equation}
\theta_e=\frac{kT_e}{m_ec^2}=\frac{T_e}{5.93\times10^9 K}
\end{equation}
is the temperature in units of the electron mass $m_e$. The peak luminosity is given by
\begin{equation}
\label{eq:lnup}
    L_{\nu_p} = 5.06\times10^{38}\frac{\textrm{erg}}{\textrm{s}\cdot\textrm{eV}}~ \alpha^{-1}(1-\beta)\left(\frac{M}{M_\odot}\right)\dot{m}^\frac{3}{2}\theta_e^5r_\textrm{min}^{-\frac{1}{2}}~.
\end{equation}
The total resulting synchrotron power is
\begin{equation}
P_\textrm{syn}= \int_0^{\nu_p} L_{\nu} d\nu \simeq 0.71\nu_p L_{\nu_p}~.
\end{equation}
For our typical parameters of interest, the synchrotron emission is poorly absorbed by the surrounding medium, since it peaks well below the hydrogen ionization threshold. Therefore, it does not significantly contribute to gas heating.

The synchrotron photons undergo inverse-Compton scattering with the surrounding electron plasma. The thermal IC spectrum can be well approximated by a power law in the frequency range $\nu_p\leq\nu\lesssim3 k T_e$ as
\begin{equation}
\label{eq:ICspectrum}
L_{\nu,\textrm{IC}} = L_{\nu_p}\left(\frac{\nu}{\nu_p}\right)^{-\alpha_c}~,
\end{equation}
where
 \begin{equation} \label{eq:alphac}
    \alpha_c = -\frac{\ln{\tau_{es}}}{\ln{A}}~.
\end{equation}
The exponent $\alpha_c$ depends on the electron scattering optical depth $\tau_{es}=12.4\dot{m}\alpha^{-1}r_\textrm{min}^{-1/2}$ and the amplification factor (mean amplification per scattering) $A=1+4\theta_e+16\theta_e^2$. The total IC power
\begin{equation}
\label{eq:ICpower}
    P_\textrm{IC}=\int_{\nu_p}^{3T_e}L_{\nu,\textrm{IC}}d\nu = \dfrac{\nu_p L_{\nu_p}}{1-\alpha_c}\left(\left[\frac{6.25\times 10^{7} \left(\dfrac{T_e}{10^{9}\textrm{ K}}\right)}{\left(\dfrac{\nu_p}{10^{12}}\right)}\right]^{1-\alpha_c}-1\right)
\end{equation}
is sensitive to the considered temperature due to the $T^5$ dependence of $L_{\nu_p}$ and the temperature dependence implicit in $\alpha_c$.

Depending on the accretion rate, different components of the ADAF spectrum can dominate. At low mass accretion rates (eADAF regime), the exponent in IC luminosity becomes $\alpha_c>2$, resulting in $P_\textrm{IC}\simeq P_\textrm{syn}/0.71(\alpha_c-1)$, and the IC power is a fraction of the synchrotron power. For moderate accretion rate (ADAF regime), $1<\alpha_c\lesssim 2$ and while IC power is the dominant emission component, the synchrotron emission is still comparable. At the highest accretion rates right below the thin disk threshold (LHAF regime), $\alpha_c < 1$ and the IC radiation becomes highly efficient. We give in Section~\ref{sssec:ADAFtemp} details of how the temperature $\theta_e$ is treated for each regime.

The third component of the ADAF emission spectrum comes from bremsstrahlung of the thermal electrons, which provides a flat spectrum with an exponential drop at $T_e$
\begin{equation}
\begin{split}
\label{eq:bremslum}
    L_{\nu,\textrm{brems}} = 1.83 \times 10^{17} \left(\frac{\alpha}{0.1}\right)^{-2} \left(\frac{c_1}{0.5}\right)^{-2} \ln{\left(\frac{r_\textrm{max}}{r_\textrm{min}}\right)} F(\theta_e) \\ \left(\frac{T_e}{5\times10^{9}~\textrm{K}}\right)^{-1}e^{-(h\nu/kT_e)} m \dot{m}^2\textrm{ ergs} \textrm{ s}^{-1}\textrm{Hz}^{-1}~,
\end{split}
\end{equation}
where $r_\textrm{max}$ is $\mathcal{O}(10^3)$ and where $F(\theta_e) \sim \mathcal{O}(1)$ is given by~\cite{Mahadevan:1996jf}
\begin{align}
\begin{split}
\label{eq:ftheta}
        F(\theta_e)=
    \begin{cases}
    4\left(\dfrac{2\theta_e}{\pi^3}\right)^\frac{1}{2} (1+1.78 \theta_e^{1.34}) + 1.73\theta_e^\frac{3}{2}(1+1.1\theta_e+\theta_e^2-1.25\theta_e^\frac{5}{2}),& \theta_e \leq 1~,\\
    \left(\dfrac{9\theta_e}{2\pi}\right)\left(\ln{[1.12\theta_e+0.48]}+1.5\right)+2.30\theta_e(\ln{[1.12\theta_e]} + 1.28),& \theta_e \geq 1~.
\end{cases}
\end{split}
\end{align}
While bremsstrahlung produces the majority of hard X-rays and gamma rays, IC typically generates more power in the relevant range of high hydrogen optical depth $10\textrm{ eV}$ to $100\textrm{ eV}$. Since the total bremsstrahlung power is also negligible relative to the combined synchrotron and IC power in the relevant parameter space, it will thus not contribute significantly to gas heating.

\subsubsection{ADAF - temperature considerations}
\label{sssec:ADAFtemp}

For ADAF, the electron temperature $\theta_e$ is determined by balancing the heating and radiation processes~\cite{Mahadevan:1996jf}. At low $\dot{m}$, direct viscous electron heating $\delta q_+$ is dominant, and is superseded by ion-electron collisional heating $q_{ie}$ at high accretion rates. This results in different temperature considerations for each of the ADAF sub-regimes.

The heating contributions are given by
\begin{align}
\begin{split}
\label{eq:adafheating}
    \delta q_+ ~=&~   9.39\times10^{38}~\dfrac{\textrm{erg}}{\textrm{s}}~ \delta\left(\frac{1-\beta}{f}\right)c_3 m\dot{m}r_\textrm{min}^{-1} \\
    q_{ie} ~=&~ 1.2\times10^{38}~\dfrac{\textrm{erg}}{\textrm{s}}~g(\theta_e)\alpha^{-2}c_1^{-2}c_3 \beta m \dot{m}^2 r_\textrm{min}^{-1}~,
\end{split}
\end{align}
where $f$ is the fraction of energy advected that is $\sim 1$ for ADAF flows, $\delta = 0.3$ is the fraction of the viscously dissipated energy that directly heats
electrons, $\beta = 10/11$ is the ratio of gas pressure to total pressure, $r_\textrm{min}=3r_s$ is the minimum flow radius set by the ISCO radius, $\alpha = 0.1$ is the viscosity parameter and constants $c_1=0.5$, $c_3=0.3$ in the formalism of Ref.~\cite{Mahadevan:1996jf} characterize the ADAF solution~\cite{Narayan:1994xi} as before. Here, $g(\theta_e)$ is an $\mathcal{O}(1)$
function
\begin{equation}
\label{eq:gtheta}
    g(\theta_e) = \frac{1}{K_2(1/\theta_e)}\left(2+2\theta_e + \frac{1}{\theta_e}\right) e^{-(1/\theta_e)}
\end{equation}
 that takes on the values of $4.5 - 0.63$ in the temperature range $2\times10^{9} \leq T_e \leq 10^{10}$.
For computational simplicity, we use the approximation $g(\theta_e) \simeq 1.19 ~\theta_e^{-1.15}$, which has a maximum error of $\sim 20\%$~\cite{Mahadevan:1996jf} over our temperature range of interest. When the two electron heating processes are equal, $\delta q_+ = q_{ie}$, the accretion rate is
\begin{equation}
\label{eq:heatingequality}
    \dot{m}_\textrm{eq} = 2.0\times10^{-3} \left(\frac{\alpha}{0.1}\right)^2\left(\frac{\delta}{0.1}\right)\left(\frac{1-\beta}{\beta}\right)\left(\frac{c_1}{0.5}\right)^2 \left(\frac{f}{1.0}\right)^{-1} g(\theta_e)^{-1}~.
\end{equation}
This is similar to our condition $\dot{M}_{\textrm{ADAF}}=10^{-3}(\alpha/0.1)^2$ of Sec.~\ref{sssec:ADAF}. Thus, in the ``eADAF" and ``ADAF" regimes, we use $\delta q_+$ heating and $q_{ie}$ heating for ``LHAF".

For $\dot{m}<\dot{M}_\textrm{ADAF}$, direct viscous heating balances the synchrotron and IC emission, which are of comparable magnitude. The temperature can be estimated as~\cite{Mahadevan:1996jf}
\begin{equation}
\label{eq:lowmdottemp}
    \theta_e \simeq 0.2{A_c^{-\frac{1}{7}}}\delta^\frac{1}{7}\alpha^\frac{3}{14}(1-\beta)^{-\frac{1}{14}}r_\textrm{min}^{\frac{3}{28}}\left(\frac{M}{M_\odot}\right)^\frac{1}{14}\dot{m}^{-\frac{5}{28}}~.
\end{equation}
Here $A_c$ is prefactor determined by the relative contributions of synchrotron and IC cooling in the total cooling rate, ranging from $0.95$ to $1.4$. As $\dot{m}$ increases, the accretion flow radiates more efficiently through IC scattering, cooling the electron plasma. We approximate $A_c^\frac{1}{7} = 1.1$ for $\dot{m}<\dot{M}_\textrm{eADAF}$ , and $A_c^\frac{1}{7} = 1.3$ for $\dot{M}_\textrm{eADAF}<\dot{m}<\dot{M}_\textrm{ADAF}$.

In the LHAF regime, the dominant IC emission is highly sensitive to temperature variations and balanced by ion electron $q_{ie}$ heating. We use an iterative scheme modified from that of Ref.~\cite{Mahadevan:1996jf} to self-consistently determine the resulting temperature\footnote{Eq. (45) in Ref.~\cite{Mahadevan:1996jf} is used in lieu of Eq. (46) (which contains a sign error), the approximations for $x_M \simeq 4\times 10^3 \dot{m}^{1/4}$ and $g(\theta_e)$ are used in Eq. (47), and the new and old values of the temperature are averaged to improve convergence.}. Starting with an initial guess of $\alpha_c = 0.75$, Eq.~\eqref{eq:ICspectrum} can be solved for the temperature,
\begin{equation}
\label{eq:tempiter}
    \theta_e = \frac{\Big(\sqrt{4\tau_{es}^{-1/\alpha_c}-3}\Big)-1}{8}~.
\end{equation}
This result is then used to solve iteratively for a new $\alpha_c$ in Eq.~\eqref{eq:alphac} as
\begin{equation}
\label{eq:alphaiter}
    \alpha_{c,\textrm{new}} = 1-\frac{\ln\left[(1-\alpha_{c,\textrm{old}})\dfrac{q_{ie}}{\nu_p L_{\nu_p}}+1\right]}{\ln C_f}~,
\end{equation}
where
\begin{align}
\begin{split}
\label{eq:alphaiter2}
    \frac{q_{ie}}{\nu_p L_{\nu_p}} ~=&~ 1.78\times10^{-5}\theta_e^{-8.15}\alpha^{-\frac{1}{2}}\beta(1-\beta)^{-\frac{3}{2}}r_\textrm{min}^{\frac{3}{4}}\left(\frac{M}{M_\odot}\right)^\frac{1}{2}\dot{m}^{-\frac{1}{4}}~,\\
    C_f ~=&~ 20.2\alpha^{\frac{1}{2}}(1-\beta)^{-\frac{1}{2}}r_\textrm{min}^{\frac{5}{4}}\theta_e^{-1}\left(\frac{M}{M_\odot}\right)^{\frac{1}{2}}\dot{m}^{-\frac{3}{4}}~.
\end{split}
\end{align}
After each iteration, the new and old values of the temperature are averaged to speed up convergence.

\subsection{Dynamical friction}
\label{ssec:dynfric}

Dynamical friction due to gravitational interactions of traversing PBHs with the surrounding medium can both heat the gas and change the PBH velocity distribution. In contrast to photon emission and outflows, dynamical friction has the advantage of directly depositing kinetic energy, thus heating the gas without any additional considerations of attenuation length for emitted particles and stopping power.

While some of the other heating mechanisms we have considered have a high degree of uncertainty, dynamical friction follows a straightforward description via the ``gravitational drag'' force given by~\cite{2008gady.book.....B}
\begin{align}
\begin{split}
\label{eq:dynfricgen}
    F_\textrm{dyn} =& -\frac{4\pi G^2 M^2 \rho}{v^2} I \\
    =& -3.95\times10^{12}\textrm{kg~m~s}^{-2} \left(\frac{M}{M_\odot}\right)^2 \left(\frac{\rho}{1\textrm{ GeV}\textrm{ cm}^{-3}}\right)\left(\frac{v}{10\textrm{ km}\textrm{ s}^{-1}}\right)^{-2} I~,
\end{split}
\end{align}
where $I$ is a velocity-dependent geometrical factor. For a collisionless medium, such as DM, this factor takes the form
\begin{equation}
\label{eq:colllessI}
    I_\textrm{cl} = \ln\left(\frac{r_\textrm{max}}{r_\textrm{min}}\right)\left[\textrm{erf}\left(\frac{v}{\sqrt{2}\sigma_{\rm DM}}\right)-\frac{2v}{\sqrt{2\pi}\sigma_{\rm DM}}e^{-\frac{v^2}{2\sigma_{\rm DM}^2}}\right]~,
\end{equation}
where $\sigma_{\rm DM}$ is the DM velocity dispersion,  $r_\textrm{max}$ denotes the size of the surrounding affected system, which we identify with a gas cloud size and $r_\textrm{min}$ is the size of the perturbing body, which we take to be the Schwarzschild radius of the BH. For a gaseous medium, $I$ depends on the Mach number $\mathcal{M} = v/c_s$, where $c_s$ is the medium sound speed, and whether the motion of the perturbing PBH is subsonic ($\mathcal{M}<1$)  or supersonic ($\mathcal{M}>1$)~\cite{1999ApJ...513..252O}
\begin{align}
\begin{split}
\label{eq:gasI}
        I=
    \begin{cases}
    \frac{1}{2}\ln\left(\frac{1+\mathcal{M}}{1-\mathcal{M}}\right)-\mathcal{M},& \mathcal{M}<1~,\\
    \frac{1}{2}\ln\left(1-\frac{1}{\mathcal{M}^2}\right)+\ln\left(\frac{r_\textrm{max}}{r_\textrm{min}}\right),& \mathcal{M}>1~.
    \end{cases}
\end{split}
\end{align}
The resulting output power due to dynamical friction for PBH traversing gas is also in this case the gas heating power $\mathcal{H}$ as the gas is directly heated, is given by
\begin{equation}
\label{eq:dynpowergen}
    \mathcal{H}_{\rm dyn} = P_\textrm{dyn} = F_{\rm dyn} v =  3.95\times10^{23} \textrm{ erg}\textrm{ s}^{-1} \left(\frac{M}{M_\odot}\right)^2\left(\frac{\rho}{1 \textrm{ GeV}\textrm{ cm}^{-3}}\right)\left(\frac{v}{10\textrm{ km}\textrm{ s}^{-1}}\right)^{-1} I~.
\end{equation}

We can estimate the change in PBH velocity due to dynamical friction using Eq.~\eqref{eq:colllessI}. With dwarf galaxies in mind, let us consider a gaseous system with DM density $\rho_{\rm DM}$ in the central region of the system that is significantly higher than the gas density $\rho = n m_p$. The effect of dynamical friction on the PBH velocity can be seen from comparing the total work done on the PBH traversing the system to its kinetic energy
\begin{equation}
\label{eq:dynfricratio}
\frac{F_\textrm{dyn}r_\textrm{max}}{\frac{1}{2}Mv^2} =-6.13\times10^{-7}\left(\frac{M}{M_\odot}\right)\left(\frac{\rho}{1\textrm{ GeV}\textrm{ cm}^{-3}}\right)\left(\frac{v}{10\textrm{ km}\textrm{ s}^{-1}}\right)^{-4}\left(\frac{r_\textrm{max}}{500\textrm{ pc}}\right) I~.
\end{equation}
The effect on the PBH velocity is insignificant throughout most of the parameter space of interest, unless we consider gaseous environments of significantly higher density (in either gas or DM) or super-massive PBHs with mass $M\gtrsim10^5 M_\odot$. We therefore ignore dynamical friction effects on PBH velocity. We note that lower relative velocities will generally increase the PBH accretion rate, resulting in an increase of the photon and outflow power delivered to the gaseous system. Hence, ignoring the velocity modification yields more conservative bounds on PBH gas heating.

\subsection{Accretion mass outflows/winds}
\label{ssec:outflows}

Mass outflows (winds) composed primarily of protons are expected to be significant for hot ADAF flows, contributing to the heating of the surrounding medium~\cite{Yuan:2014gma}. Unlike the highly beamed jets\footnote{As jets are typically associated with Kerr black holes, they would require a separate treatment and we do not consider them here.}, outflows are not highly relativistic and cover wider angular distributions. However, there is significant uncertainty in the description of outflows.

The outflows reduce the in-falling matter at smaller radii and the accretion rate in the presence of outflows can be approximately modelled by a self-similar power-law~\cite{Blandford:1998qn}
\begin{equation}
\label{eq:inflowgen}
    \dot{M}_\textrm{in}(r) = \dot{M}_\textrm{in}(r_\textrm{out})\left(\frac{r}{r_\textrm{out}}\right)^s~,
\end{equation}
where the index $0 \leq s < 1$ is limited by energy and mass conservation. Here, $\dot{M}_\textrm{in}(r_\textrm{out})$ is the Bondi-Hoyle-Lyttleton accretion rate of Eq.~\eqref{eq:bondihoyle}. The corresponding outflow (wind) rate is~\cite{2006PASJ...58..965T}
\begin{equation}
\label{eq:outflowgen}
    \dot{M}_w(r) = r\frac{d\dot{M}_\textrm{in}}{dr} = s\dot{M}_\textrm{in}(r_\textrm{out})\left(\frac{r}{r_\textrm{out}}\right)^s~.
\end{equation}
From simulations, the exponent $s$ is estimated to be $0.4-0.8$~\cite{2012ApJ...761..130Y,Yuan:2014gma}.

A wide range of values has been considered in the literature for the outer radius $r_\textrm{out}$, including $100 r_s$~\cite{Xie:2012rs}, $500 r_s$~\cite{2012ApJ...761..130Y},  $r_B$~\cite{2011MNRAS.415.1228P}, to $10 r_B$~\cite{2013ApJ...767..105L}. The resulting outgoing wind has a velocity that is a fraction $f_k\simeq0.1-0.2$~\cite{2013ApJ...767..105L,2013IAUS..290...86Y,2012ApJ...761..130Y} of the Keplerian velocity $v_k$ at the radius at which it is ejected (due to $Be/v_k^2$, where $Be$ is the Bernoulli parameter, being roughly constant~\cite{2012ApJ...761..130Y}), i.e. $v(r) \simeq f_k \sqrt{G M/r}$. Most of the kinetic energy of the wind comes from the inner region, with the kinetic energy per particle determined at the ISCO radius with $v_{\rm in} = v(r_{\rm in} = 3 r_s)$ being
\begin{equation}
E_\textrm{kin} \simeq \frac{1}{2} m_p v_{\rm in}^2 = 782\textrm{ keV} \left(\dfrac{r_\textrm{in}}{3r_s}\right)^{-1}\left(\dfrac{f_v}{0.1}\right)^2~,
\end{equation}
which is independent of PBH mass.

We derive the energy distribution of the wind, $f(E)$, and integrate it with the energy deposited per proton $\Delta(E)$. The heat deposited is limited by the value of the energy loss per unit length $x$ stopping power $dE/dx = nS(E)$, which we take from Ref.~\cite{2019PhRvA..99d2701B} (see their Fig.~9),
\begin{equation}
\label{eq:energydepo}
    \Delta E = \int \frac{dE}{dx} dx \simeq \textrm{min}(E,n S(E) r_\textrm{max})~.
\end{equation}
The total heat deposited in the cloud is then
\begin{equation}
\label{eq:outflowint}
    \mathcal{H}_\textrm{outflow} = \int_{r_\textrm{in}}^{r_\textrm{out}} f_h \Delta E \frac{\mathcal{D}}{\mu m_p}\frac{d\dot{M}_\textrm{in}}{dr} dr = \int_{E(r_\textrm{min})}^{E(r_\textrm{out})} f_h \Delta E f(E) dE~,
\end{equation}
where we estimate the fraction of energy deposited as heat to be $f_h\sim1/3$ and
where the distribution function is given by
\begin{align}
\begin{split}
    \label{eq:fe}
    f(E) ~=&~ s\dot{M}_\textrm{in}(r_\textrm{out})\left(\frac{r}{r_\textrm{out}}\right)^{s}\frac{1}{E} \\
    ~=&~ 3.54\times10^{22}~\frac{\textrm{ erg}}{\textrm{ eV}^{2}\textrm{ sec}}~ s(2.35\times10^{8})^s \left(\frac{M}{M_\odot}\right)^2 \left(\frac{n}{1\textrm{ cm}^{-3}}\right)\\
    &~\times\left(\frac{\Tilde{v}}{10\textrm{ km}/\textrm{s}}\right)^{-3} \left(\frac{\mathcal{D}}{1}\right)\left(\frac{r_\textrm{out}}{r_s}\right)^{-s}f_v^{2s}E^{-s-1}~.
\end{split}
\end{align}
We have chosen a duty cycle factor of $\mathcal{D}= 1$, as discussed below.

Magnetic fields in the ISM can trap the protons from the accretion disk wind, effectively increasing the stopping power and the heat deposition. Since the strength, orientation, and distribution of magnetic fields in Leo T are not known\footnote{See related discussion for positron propagation, which even in the case of Milky Way is highly uncertain~\cite{Panther:2018xvc}.}, we conservatively omit any such factors.

\subsubsection{Duty cycle}
\label{sec:duty}

The wind outflow feedback can affect the PBH accretion and hence the emission efficiency. If the outflow feedback is very strong, heated ISM at the Bondi radius will be blown away, resulting in reduced accretion. On the other hand, if accretion is terminated so will be the wind and feedback.
Wind duty cycle can be estimated following~Ref.~\cite{Ioka:2016bil} by comparing the dynamical and accretion timescales. If the accretion timescale is significantly longer than the dynamical timescale, we do not expect accretion to be efficient.

The dynamical time of the system (i.e., the crossing time) is defined as
\begin{equation}
   t_{\rm dyn} = \frac{r_{\rm sys}}{v},
\end{equation}
while the accretion time scale (i.e., free-fall time scale) is defined as
\begin{equation}
	t_{\rm ff} = \frac{r_B}{v_{\rm ff}}=\frac{GM}{\sqrt{2}v^3}~,
\end{equation}
where $v_{\rm ff} = \sqrt{2GM/r_B}$ is the free-fall velocity, and $r_B=GM/v^2$ is the Bondi radius as before. Then we can consider the duty cycle $\mathcal{D}$ to be
\begin{equation}
\mathcal{D} = \max\left[1, \frac{t_{\rm dyn}}{t_{\rm ff}}\right]	= \max\left[1, \frac{\sqrt{2}r_{\rm sys}v^2}{GM}\right]~.
\end{equation}
We plot the duty cycle $\mathcal{D}$ on Fig.~\ref{fig:dutyc}, assuming characteristic $v=10$~km/s (here we take $\tilde{v} \sim v$) and $r_{\rm sys}=1$~pc values. For $M\lesssim 10^{4} M_\odot$, we estimate the factor to be 1, i.e. the duty cycle has no significant effect on Eq.~\eqref{eq:fe}. Above $M \gtrsim 10^4 M_{\odot}$, we find PBHs to be constrained by the incredulity limit (see Eq.~\eqref{eq:incredlimit}), which justifies us in taking $\mathcal{D}=1$ throughout.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=.5\textwidth]{plots/duty_cycle.pdf}
\caption{PBH accretion wind outflow duty cycle $\mathcal{D}$, for characteristic $v = 10$~km/s and system radius $r_{\rm sys} = 1$~pc.}\label{fig:dutyc}
\end{center}
\end{figure}

\section{Astrophysical Systems}

We now apply our analysis to two examples of astrophysical systems, namely Milky Way gas clouds and dwarf galaxies. Although we limit our analysis to these systems, the methods outlined here can be readily extended to other gas clouds.  Requiring the thermal balance of the gas cloud between standard cooling processes and heating from PBH, we constrain the PBH mass fraction $f_\textrm{PBH}$. A similar approach has been used to set limits on particle DM~\cite{Bhoonah:2018wmw,Farrar:2019qrv,Wadekar:2019xnf}. We find that for PBHs the different heating processes are particularly favorable to explore within the dwarf galaxy Leo T.

Natural heating sources (e.g. stellar radiation) are for simplicity neglected in this study, which results in a more conservative upper bound on PBH heating. Since our argument relies on maintaining thermal equilibrium, only gas systems that are expected to be stable on sufficiently long timescales $\tau_\textrm{sys}$ are used (see discussion in Ref.~\cite{Farrar:2019qrv}).  We thus require the timescale for which the gas system remains stable to be longer than the cooling timescale of the gas $\tau_\textrm{therm}$
\begin{equation}
\label{eq:timescales}
    \tau_\textrm{sys} \gg \tau_\textrm{therm} = \frac{3nkT}{2\dot{C}}~,
\end{equation}
where $k$ is the Boltzmann constant and $\dot{C}$ is the gas cooling rate per volume.

Furthermore, the disk formation timescale needs to be smaller than the system crossing timescale. This is directly analogous to the duty cycle considerations of Section~\ref{sec:duty} and is satisfied for our systems of interest. In the case of Leo T, this is not strictly necessary, since PBHs may orbit about the center of the dwarf galaxy. On the other hand, this could affect Milky Way cloud considerations, as PBHs do not spend all of their time inside them.

In principle a multitude of processes can participate in gas temperature exchange, and there exist detailed numerical analyses involving a full chemistry network~\cite{Smith:2016hsc}. We employ approximate results obtained in Ref.~\cite{Wadekar:2019xnf} for our gas systems of interest. The volumetric cooling rate (erg cm$^{-3}$ s$^{-1}$) of atomic hydrogen gas is given by
\begin{equation}
\label{eq:coolinggen}
    \dot{C} = n^2 10^{[\textrm{Fe/H}]}\Lambda(T)~,
\end{equation}
where [Fe/H]$ \equiv \log_{10}(n_{\rm Fe}/n_{\rm H})_{\rm gas} - \log_{10}(n_{\rm Fe}/n_{\rm H})_{\rm Sun}$ is the metallicity, and $\Lambda(T) \propto 10^{[{\rm Fe/H}]}$ is the cooling function which can be fitted numerically. Using the results from Ref.~\cite{Smith:2016hsc} library, Ref.~\cite{Wadekar:2019xnf} obtained $\Lambda(T) = 2.51\times10^{-28}T^{0.6}$~(erg cm$^3$ s$^{-1}$), valid in the range $300~\text{K} < T < 8000~\text{K}$.

The total heat deposited in the cloud by the PBH
\begin{equation}
  H_{\rm tot} = N_\textrm{PBH} H(M)= f_\textrm{PBH}\dfrac{\rho_\textrm{DM}V H(M)}{M}
\end{equation}
is required to be less than the total cooling rate $\dot{C} V$, yielding the condition on the PBH mass fraction
\begin{equation}
\label{eq:genbound}
    f_\textrm{PBH} < f_\textrm{bound} =\frac{M\dot{C}}{\rho_\textrm{DM}H(M)}~.
\end{equation}
Here $H(M)$ is the average heat deposited by one PBH of mass $M$.

\subsection{Gas systems with bulk relative velocity with respect to the DM}
\label{app:relvel}

A gas system with a bulk relative velocity $\vec{v}_b$ with respect to the Galactic frame, and thus the PBH in the DM, can be treated as follows.
Without loss of generality, consider the frame in which the bulk velocity is rotated to the z direction. Then the one dimensional Maxwellian velocity distributions are
\begin{align}
\begin{split}
\label{eq:1dmaxwellian}
    \begin{cases}
    f_x(v_x)=f_y(v_y)=\dfrac{1}{\sqrt{2\pi}\sigma_v}e^{-\frac{v_x^2}{2\sigma_v^2}}\\ \\
    f_z(v_z) = \dfrac{1}{\sqrt{2\pi}\sigma_v}e^{-\frac{(v_z-v_b)^2}{2\sigma_v^2}}
    \end{cases}
\end{split}
\end{align}
Imposing the constraint $v^2 = \sqrt{v_x^2+v_y^2+v_z^2}$, the relative velocity distribution is given by the integral
\begin{align}
\label{eq:maxwellintegral}
    f(v) ~=&~ \frac{1}{(2\pi)^{\frac{3}{2}}\sigma_v^3}\int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y\int_{-\infty}^{\infty}dv_z \delta(v^2-v_x^2-v_y^2-v_z^2)e^{-\frac{v_x^2+v_y^2+(v_z-v_b)^2}{2\sigma_v^2}} \\
    =&~ \frac{e^{-\frac{v^2+v_b^2}{2\sigma^2}}}{(2\pi)^{\frac{3}{2}}\sigma_v^3}\int_{-v}^{v}dv_x\int_{-\sqrt{v^2-v_x^2}}^{\sqrt{v^2-v_x^2}}dv_y \frac{v}{\sqrt{v^2-v_x^2-v_y^2}}  \left[e^{-\frac{v_b\sqrt{v^2-v_x^2-v_y^2}}{\sigma_v^2}} + e^{-\frac{v_b\sqrt{v^2-v_x^2-v_y^2}}{\sigma_v^2}}\right]~. \notag
\end{align}
Switching to polar velocity coordinates and simplifying,
\begin{align}
\begin{split}
\label{eq:modmaxwellian}
    f(v) ~=&~ \frac{2\pi v}{(2\pi)^\frac{3}{2}\sigma_v^3}e^{-\frac{v^2+v_b^2}{2\sigma_v^2}}\int_0^{v}dv_r \frac{v_r}{\sqrt{v^2-v_r^2}}\left[e^{-\frac{v_b\sqrt{v^2-v_r^2}}{\sigma_v^2}}+e^{\frac{v_b\sqrt{v^2-v_r^2}}{\sigma_v^2}}\right] \\
    ~=&~ \frac{v}{\sqrt{2\pi}\sigma v_b}\left[e^{-\frac{(v-v_b)^2}{2\sigma^2}}-e^{-\frac{(v+v_b)^2}{2\sigma^2}}\right]~.
\end{split}
\end{align}
In the limit $v_b\rightarrow 0$, Eq.~\eqref{eq:modmaxwellian} reproduces the standard Maxwellian distribution Eq.~\eqref{eq:maxwellian}.

\subsection{Milky-Way gas clouds}

\begin{figure}[tb]
\begin{center}
\includegraphics[trim={5mm 0mm 40 0},clip,width=.475\textwidth]{plots/LuminosityG33480.pdf}
\includegraphics[trim={5mm 0mm 40 0},clip,width=.475\textwidth]{plots/fallowed2G33480.pdf}

\caption{  \label{fig:G33480graphs}  \textbf{Left:}
Amount of heat absorbed by the Milky-Way gas cloud G33.4-8.0 from a single PBH of mass $M$. \textbf{Right:} Constraints from G33.4-8.0 on the DM fraction in PBHs. The reach of constraints is bounded by the positive slope solid black line indicating the incredulity limit. In both panels three heating components are shown: photon emission (red), dynamical friction (green), mass outflows (blue), as well as the total heat (black dashed). Both Model 1 and Model 2 of emission are shown, and the variation between them is shaded in blue. The observed line features are due to transitions between the eADAF, ADAF, LHAF, and thin disk regimes.}
\end{center}
\end{figure}

A variety of interstellar neutral hydrogen (HI) gas clouds exist in the inner Galaxy within a few hundred pc of the Galactic plane~\cite{2015ApJS..219...16P}.
In general, such Milky-Way clouds are not ideal candidates for studying PBHs due to the high relative rotational velocity, high velocity dispersion, and lower total DM mass (relative e.g. to dwarf galaxies) due to their smaller size. The last of these shortcomings limits the maximum mass considered for PBHs due to the incredulity limit of Eq.~\eqref{eq:incredlimit}. Thus, even if a cloud had favorable properties otherwise, it would not be able to provide significant constraints in the intermediate PBH mass range.

To demonstrate our analysis, we consider as an example a particular Milky-May cloud, G33.4-8.0. To describe G33.4-8.0, we use parameters in Refs.~\cite{2015ApJS..219...16P,2015ApJS..216...29B,2004ApJ...602..162H,2002ApJ...580L..47L,2010ApJ...722..367F,2001ApJ...559..318R,Wadekar:2019xnf}. Namely, the atomic hydrogen gas density is taken to be $n=0.4~\pm 0.1\textrm{ cm}^{-3}$, the cloud radius $r\simeq 30\textrm{ pc}$, the bulk relative velocity $v_{\rm bulk} = 220\textrm{ km/s}$, the temperature $T=400~\pm 90$~K, and the corresponding adiabatic sound speed $c_s = 2.4 \textrm{ km/s}$, and the cooling rate of $\simeq 2.1\times10^{-27} \textrm{ erg}\textrm{ cm}^{-3}\textrm{ s}^{-1}$. The DM density of the cloud is $\rho = 0.64\textrm{ GeV}/ \textrm{cm}^{3}$, assuming Navarro-Frenk-White (NFW) profile~\cite{Navarro:1995iw} for the Galactic DM halo, with the DM velocity dispersion of $124.4\textrm{ km/s}$.

In Fig.~\ref{fig:G33480graphs}, we display the heating from a single PBH for G33.4-8.0, including contributions due to photon emission, dynamical friction and mass outflows. To obtain the limits in Fig.~\ref{fig:G33480graphs}, we numerically integrate Eq.~\eqref{eq:totalheateq}. To treat the PBH emission and its uncertainty, we consider two distinct input parameter sets\footnote{We do not consider a variation in PBH emission parameter $\delta$ between $0.1-0.5$ for our analysis, since it only corresponds to a $\mathcal{O}(\%)$ variation in the photon luminosity.} given in Table~\ref{tab:param}. The observed irregularities in the figure are due to transitions between the eADAF, ADAF, LHAF, and thin disk regimes.

We also display in Fig.~\ref{fig:G33480graphs} the combined resulting constraints from G33.4-8.0 gas heating, bound by the incredulity limit
\begin{equation}
    f_\textrm{bound} > 5.25\times10^{-4}\left(\frac{M}{M_\odot}\right)\left(\frac{r_\textrm{sys}}{30\textrm{ pc}}\right)^{-3}\left(\frac{\rho_\textrm{DM}}{0.64\textrm{ GeV}\textrm{ cm}^{-3}}\right)^{-1}~.
\end{equation}
The largest PBH mass that can be bounded by G33.4-8.0 with this method is $1.9\times10^{3} M_\odot$. The weakness of Milky Way gas cloud bounds is due in large part to the high velocity dispersion as well as the co-rotation of the gas clouds in the Milky Way, which adversely affects accretion disk formation. In addition, the low metallicity of Leo T contributes to reduced cooling rates relative to G33.4-8.0 enhancing the effect of PBH heating.

\begin{table}[tb]
  \setlength{\extrarowheight}{2pt}
  \setlength{\tabcolsep}{5pt}
  \begin{center}
	\begin{tabular}{c|c|c|c|c|c|c|c}  \hline\hline
	   & $\delta$ & $\beta$ & $\alpha$ & $s$ & $r_{\rm outer}$ & $f_v$ & $\mathcal{D}$\\
	\hline
	Model~1 & 0.3 & 10/11 & 0.1 & 0.5 & $r_B$ & 0.1 & 1\\  \hline
	Model~2 & 0.3 & 10/11 & 0.1 & 0.7 & $100r_s$ & 0.2 & 1\\   \hline\hline
	\end{tabular}
  \end{center}
\caption{\label{tab:param} Input parameters for PBH emission.}
\end{table}

\subsection{Dwarf galaxies}

Dwarf DM-rich galaxies provide a good setting to study DM-gas heating. In contrast to Milky-Way gas clouds, their lack of relative bulk velocity with respect to the DM as well as significant size that allows to weaken the incredulity limit makes them favorable systems for PBH heating.

We will focus on Leo T, a transition type of galaxy between a dwarf irregular and a spheroidal, which is
a few hundred kpc away from the Milky Way. Leo T has been well studied and modeled theoretically and it has the desired properties, including low baryon velocity dispersion that is also favorable for PBH heating.

In the inner $r\lesssim 350\textrm{ pc}$ of Leo T, the gas is mostly composed of atomic hydrogen, but consists of highly ionized hydrogen outside this central region~\cite{2013ApJ...777..119F}. To minimize the cooling rate, we avoid the efficiently cooled ionized region (see e.g.~Ref.~\cite{Smith:2016hsc}) and focus our analysis on the central region. The atomic hydrogen number density ranges from a maximum of $\sim 0.2\textrm{ cm}^{-3}$ at the center decreasing to $0.03\textrm{ cm}^{-3}$ at the boundary, at $r=350\textrm{ pc}$~\cite{2013ApJ...777..119F}. We estimate the density to be a constant $0.07\textrm{ cm}^{-3}$ (roughly the average density in the central region), and justify this by noting that both the heating processes and cooling rate scale approximately as $n^2$. This average density of gas produces the column density of $n r_\textrm{sys} = 7.56\times10^{19}\textrm{ cm}^{-2}$, in the right range for Leo T.  Similarly, we assume a constant DM density of $1.75\textrm{ GeV}/\textrm{cm}^3$ which accounts for the total DM mass in Leo T of $10^{7} M_\odot$.

The hydrogen gas in the central region of Leo T is comprised both of a non-rotating warm component with velocity dispersion of $\sigma_g=6.9~ \textrm{km}/\textrm{s}$ and temperature $T\simeq 6000$~K~\cite{2008MNRAS.384..535R,2013ApJ...777..119F} as well as a subdominant cold component. For simplicity, we treat the entire hydrogen density as warm gas which results in a more conservative bound. From the adiabatic formula, the sound speed is estimated to be $c_s\simeq 9\textrm{ km/s}$ for a $6000 K$ gas. The DM velocity dispersion is expected to be the same as the gas dispersion, so that $\sigma_{v}=\sigma_g$. Using the above input parameters as well as metallicity of $\textrm{[Fe/H]} \simeq -2$~\cite{2008ApJ...685L..43K} and assuming that gas metallicity is related to stellar metallicity in Eq.~\eqref{eq:coolinggen}, the resulting Leo T's cooling rate is taken to be $\dot{C} = 2.28\times10^{-30}\textrm{ erg}\textrm{ cm}^{-3}\textrm{ s}^{-1}$.

\begin{figure}[tb]
\begin{center}
\includegraphics[trim={5mm 0mm 40 0},clip,width=.475\textwidth]{plots/Luminosityleot.pdf}
\includegraphics[trim={5mm 0mm 40 0},clip,width=.475\textwidth]{plots/fallowed2leot.pdf}
\caption{  \label{fig:leotgraphs}  \textbf{Left:}
The amount of heat absorbed by Leo T dwarf galaxy from a single PBH of mass $M$. \textbf{Right:} Constraints from Leo T on the DM fraction in PBHs. The reach of constraints is bounded by the positive slope solid black line indicating the incredulity limit. The variation in the constraint from Models 1 and 2 are shaded in blue. The black line shows the combined constraint from Model 2. For both panels three heating components are shown: photon emission (red), dynamical friction (green), mass outflows (blue), as well as the total heat (black). The observed line features are due to transitions between the eADAF, ADAF, LHAF, and thin disk regimes. Both Model 1 and Model 2 of emission are shown, and the difference is shaded in blue.}
\end{center}
\end{figure}

In Fig.~\ref{fig:leotgraphs}, we display the heating from a single PBH for Leo T, including contributions due to photon emission, dynamical friction as well as mass outflows. As before, we numerically integrate Eq.~\eqref{eq:totalheateq} and consider PBH emission parameter variation as described in Tab.~\ref{tab:param}. We note that for the Leo T parameters, the condition for thin disks, $\dot{m}>0.07\alpha$, is satisfied only for very massive PBHs $M\gtrsim 10^{5}M_\odot$.

We also display in Fig.~\ref{fig:leotgraphs} the combined resulting Leo T gas heating constraints, bound by the incredulity limit
\begin{equation}
    f_\textrm{bound} > 1.21\times10^{-7}\left(\frac{M}{M_\odot}\right)\left(\frac{r_\textrm{sys}}{350\textrm{ pc}}\right)^{-3}\left(\frac{\rho_\textrm{DM}}{1.75\textrm{ GeV}\textrm{ cm}^{-3}}\right)^{-1}~.
\end{equation}
The largest PBH mass that can be bound by Leo T by this method is $8.26\times10^{6}M_\odot$.

\begin{figure}[tb]
\begin{center}
\includegraphics[trim={5mm 0mm 40 0},clip,width=0.65\textwidth]{plots/fallowedleot.pdf}

\caption{  \label{leotfinalgraph} Constraints from Leo T dwarf galaxy on PBH gas heating are shown in blue. The light blue shaded band denotes variation due to PBH emission parameters. These constraints are bounded by the PBH incredulity limit (solid blue positive slope line). Other existing constraints (see Ref.~\citep{Carr:2020gox})  are shown by dashed lines, including Icarus~\cite{2018PhRvD..97b3518O}
(I) caustic crossing in purple, Planck~\cite{Ali-Haimoud:2016mbv,Serpico:2020ehh}
(P) in yellow, X-ray binaries~\cite{Inoue:2017csr}
(XRB) in green, dynamical friction of halo objects
(DF) in red, Lyman-$\alpha$~\cite{2019PhRvL.123g1102M}
(Ly-$\alpha$) in maroon, combined bounds from the survival of astrophysical systems in Eridanus II~\cite{Brandt:2016aco},
 Segue 1~\cite{2017PhRvL.119d1102K},
 and disruption of wide binaries~\cite{2014ApJ...790..159M}
(S) shown in magenta, large scale structure~\cite{Carr:2018rid}
(LSS) in cyan, and X-ray/radio~\cite{Ricotti:2007au}
(X/R) in brown.}
\end{center}
\end{figure}

\section{Conclusions}
\label{sec:conclusions}

PBHs in the stellar as well as intermediate $\sim 10-10^3 ~M_{\odot}$ mass-range have garnished significant interest and could potentially be associated with LIGO gravity wave events. As PBHs traverse the ISM, they would interact with and heat the surrounding gas. We have considered several generic heating mechanisms, including accretion photon emission, dynamical friction as well as mass outflows/winds. The balance of system heating and cooling allows for a general test of PBHs in a cosmology-independent way. As a demonstration, we have applied our analysis to Milky-Way gas clouds and dwarf galaxies. Considering constraints on gas heating for a given gas system, we have derived new bounds from Leo T dwarf galaxy on PBHs contributing to DM across a broad mass-range of $M_{\rm PBH} \sim \mathcal{O}(1) M_{\odot}-10^7 M_{\odot}$. These limits do not depend on cosmology and establish a novel independent test of PBHs constituting part of the DM in the stellar and intermediate PBH mass range. Our analysis can be readily applied to other astrophysical systems of interest.

\acknowledgments
\addcontentsline{toc}{section}{Acknowledgments}

The work of G.B.G., A.K., V.T., and P.L. was supported by the U.S. Department of Energy (DOE) grant No. DE-SC0009937. This work was also supported by the MEXT
Grant-in-Aid for Scientific Research on Innovative Areas (No.~20H01895, No.~21K13909 for K.H.). A.K. was also supported by Japan Society for the Promotion of Science (JSPS) KAKENHI grant No.
JP20H05853. Y.I. is supported by JSPS KAKENHI grant No. JP18H05458, JP19K14772, program of Leading Initiative for Excellent Young Researchers, MEXT, Japan, and RIKEN iTHEMS Program. A.K., V.T., and Y.I. are also supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan.


\section*{INTRODUCTION}

Nonlocality, which challenges our comprehension and intuition about the
nature, is a key and distinctive feature of quantum world. Three different
types of nonlocal correlations: Bell nonlocality \cite{Bell},
Einstein-Podolsky-Rosen (EPR) steering \cit
{EPR1,EPR2,Review1,Review2,RMP2009} and entanglement \cite{Entanglement}
have opened an epoch of unrelenting exploration of quantum correlations
since they were introduced. The notion of EPR steering was introduced as the
phenomenon that one can remotely steer the (conditional) quantum state owned
by the other party through local measurements on half of an entangled state
\cite{EPR1,EPR2}. According to the hierarchy of nonlocality, EPR steering
stands between entanglement and Bell nonlocality \cite{EPR from info}. The
existence of steering does not imply the violation of any Bell inequality,
while the violation of at least one Bell inequality immediately implies
steering in both directions which is referred as two-way steering \cit
{OneWay1}. Remarkably, there is a situation that a quantum state may be
steerable from Alice to Bob, but not vice versa, which is called the one-way
steering \cite{OneWay1,OneWay2,DME,PRAQin,Caiyin,NC,PRL118,OneWay4,OneWay3}.
This intriguing feature comes from the intrinsically asymmetric with respect
to the two subsystems.

The one-way EPR steering has been observed for both continuous variable (CV)
\cite{OneWay1,OneWay2,DME,PRAQin,Caiyin} and discrete variable (DV) systems
\cite{NC,PRL118,OneWay4,OneWay3}. Gaussian states and Gaussian operations
play a central role in the analysis and implementation of CV quantum
technology \cite{Adesso,xin3}. When considering generally bipartite Gaussian
states and using Gaussian measurements, the steerability between the
submodes is referred as Gaussian steering. What's more interesting, it has
been shown that the direction of Gaussian EPR steering can be manipulated in
noisy environment \cite{PRAQin}. The temporal quantum steering and
spatio-temporal steering have also been investigated \cite{s1,PRA93,s2}.
Very recently, the EPR steering using hybrid CV and DV entangled state \cit
{hybrid steering} and remote generation of Wigner-negativity through EPR
steering \cite{Wigner-neg} are demonstrated. As applications, EPR steering
has been identified as an essential resource for one-sided
device-independent quantum key distribution \cite{1sQKD2,1sQKD3,1sQKD4},
quantum secret sharing \cite{QSS3,QSS2}, secure quantum teleportation \cit
{Tel1,Tel2,Tel3}, securing quantum networking tasks \cite{PRA99} and
subchannel discrimination \cite{Subchannel1,Subchannel2}.

Besides two-mode EPR steering, the multipartite EPR steering has also been
widely investigated since it has potential application in quantum network.
Genuine multipartite steering exists if it can be shown that a steering
nonlocality is necessarily shared among all observers \cite{xin10}. The
criterion for multipartite EPR steering and for genuine multipartite EPR
steering has been developed \cite{xin8,xin10}. Genuine high-order EPR
steering has also been experimentally demonstrated \cite{xin11}.
Multipartite Gaussian EPR steering in Greenberger-Horne-Zeilinger (GHZ)
state \cite{OneWay2,OCdeng} and a four-mode cluster state \cite{DME} have
been demonstrated. The experiments of multipartite Gaussian EPR steering
only demonstrated EPR steering in the case of bipartite splitting at
present, while genuine multipartite EPR steering for Gaussian states has not
been demonstrated. Recently, the distribution of multipartite steering in a
quantum network by separable states has been demonstrated \cit
{xiangyu,wangmh}.

The performance of modern quantum communication sort of relies on the
interaction between quantum system and environment, which is generally
described by the theory of decoherence. Decoherence, caused by the loss and
noise encountered in quantum channel, can destroy fragile quantum properties
of quantum states, and thus the information carried by the quantum states
will be affected. Decoherence effects on entangled state has been well
studied \cite{ESD1,ESD2,esd3,esd4,pra83,esd5}. It has been shown that two
entangled qubits become completely separable in a finite-time under the
influence of vacuum which is the so called entanglement sudden death \cit
{ESD1,ESD2}, and disentanglement occurs in Gaussian multipartite entangled
states when one mode is transmitted in a noisy channel \cite{NonM3}. It is
essential to recover entanglement when entanglement sudden death happens.\ A
number of methods to recover the destroyed entanglement have been
demonstrated, such as the non-Markovian environment \cite{NonM1,NonM2,NonM3
, weak measurement \cite{weakm} and feedback \cite{feedb}. Since EPR
steering is so intriguing and has enormous potential in quantum
communication applications, it is imperative and significant to investigate
the influence on Gaussian EPR steering made by decoherence. In 2015, the
decoherence of steering was theoretically investigated \cite{Josab15}.
Further more, it is unknown whether the methods used to recover entanglement
can be used to recover EPR steering when sudden death of EPR steering
happens.

Here, we experimentally investigate properties of Gaussian EPR steering for
a two-mode squeezed state (TMSS) when it is transmitted in lossy and noisy
quantum communication channels, respectively. Firstly, we investigate the
influence of purity of quantum states on Gaussian EPR steering. We show that
in a lossy channel, the impurity of the state leads to decrease of
steerabilities and two-way steering region, but it never leads to sudden
death of Gaussian EPR steering. Secondly, we find the Gaussian EPR steering
totally disappear at certain noise level in a noisy quantum channel, which
demonstrates sudden death of EPR steering. Thirdly, we successfully revive
the Gaussian EPR steering after its sudden death by establishing a
correlated noisy channel. We also confirm that the steering quantifier
proposed by I. Kogias et al. and the steering criterion proposed by M. Reid
are equivalent for a TMSS with Gaussian measurement. The presented results
provide useful references for applying Gaussian EPR steering in noisy
environment.

\section*{Experimental scheme}

\textbf{The principle}

\begin{figure*}[tbp]
\begin{center}
\includegraphics[width=120mm]{Fig1.eps}\vspace*{-0.3cm}
\end{center}
\caption{The principle of experiment. \textbf{a,} Schematic of EPR steering
in lossy channel. One mode of the two-mode squeezed state (TMSS) $\hat{A}$
is kept by Alice, the other mode of TMSS $\hat{B}$ is distributed to Bob
over a lossy channel which is experimentally mimicked by a beam-splitter,
the vacuum mode couples into system from the other input port of the
beam-splitter. The transmission efficiency of the lossy channel is $\protec
\eta _{L}$. \textbf{b,} Schematic of disappearance and revival of EPR
steering. Mode $\hat{B}$ is distributed to Bob over a noisy channel with an
excess noise higher than the vacuum noise. The transmission efficiency of
the noisy channel is $\protect\eta _{N}$. An ancillary beam with noise
correlated to noisy channel couples with the transmitted mode on a
beam-splitter with transmission coefficient $T$. \textbf{c,} The
experimental setup. A pure -3 dB TMSS at the sideband frequency of 3 MHz is
generated by a nondegenerate optical parametric amplifier (NOPA). Alice and
Bob perform homodyne measurement on their states respectively. An ancillary
beam is used to revive the EPR steering. PBS, polarization beam splitter;
LO, local oscillator; EOM, electro-optic modulator; HD, homodyne.}
\end{figure*}

The state we use in the experiment is a deterministically prepared TMSS,
whose quantum correlations between quadratures are expressed by $\Delta
^{2}\left( \hat{x}_{A}+\hat{x}_{B}\right) =\Delta ^{2}\left( \hat{p}_{A}
\hat{p}_{B}\right) =2e^{-2r}$, where $\hat{x}=\hat{a}+\hat{a}^{\dagger }$
and $\hat{p}=\left( \hat{a}-\hat{a}^{\dagger }\right) /i$ are the amplitude
and phase quadratures of an optical mode respectively, $r$ is the squeezing
parameter ranging from 0 to infinite which correspond to no squeezing and
the ideal squeezing respectively. Under this definition, the noise of the
vacuum state is normalized to 1. All Gaussian properties of a TMSS can be
described by its covariance matri
\begin{equation}
\sigma _{AB}=\left(
\begin{array}{cccc}
\alpha  & 0 & \gamma  & 0 \\
0 & \alpha  & 0 & -\gamma  \\
\gamma  & 0 & \beta  & 0 \\
0 & -\gamma  & 0 & \beta
\end{array
\right) =\left(
\begin{array}{cc}
\mathbf{A} & \mathbf{C} \\
\mathbf{C}^{T} & \mathbf{B
\end{array
\right)
\end{equation
with $\sigma _{ij}=Cov\left( \hat{R}_{i},\hat{R}_{j}\right) =\frac{1}{2
\left\langle \hat{R}_{i}\hat{R}_{j}+\hat{R}_{j}\hat{R}_{i}\right\rangle
-\left\langle \hat{R}_{i}\right\rangle \left\langle \hat{R}_{j}\right\rangle
$, $i,j=1,2,3,4$, where $\hat{R}=(\hat{x}_{A},\hat{p}_{A},\hat{x}_{B},\hat{p
_{B})^{T}$ is a vector composed by the amplitude and phase quadratures of
bipartite optical beams \cite{Adesso}. Thus the TMSS covariance matrix can
be partially expressed as (the cross correlations between different
quadratures of one mode are taken as $0$)
\begin{equation}
\sigma =\left[
\begin{array}{cccc}
\bigtriangleup ^{2}\hat{x}_{A} & 0 & Cov\left( \hat{x}_{A},\hat{x
_{B}\right)  & 0 \\
0 & \bigtriangleup ^{2}\hat{p}_{A} & 0 & Cov\left( \hat{p}_{A},\hat{p
_{B}\right)  \\
Cov\left( \hat{x}_{A},\hat{x}_{B}\right)  & 0 & \bigtriangleup ^{2}\hat{x
_{B} & 0 \\
0 & Cov\left( \hat{p}_{A},\hat{p}_{B}\right)  & 0 & \bigtriangleup ^{2}\hat{
}_{B
\end{array
\right] .
\end{equation
For the theoretical covariance matrix elements of the TMSS we used in the
experiment, $\mathbf{A}=\alpha \mathbf{I}$, $\mathbf{B}=\beta \mathbf{I}$,
\mathbf{C}=\gamma \mathbf{Z}$, $\mathbf{I}$ and $\mathbf{Z}$ are the Pauli
matrice
\begin{equation}
\mathbf{I}=\left(
\begin{array}{cc}
1 & 0 \\
0 &
\end{array
\right) ,\mathbf{Z}=\left(
\begin{array}{cc}
1 & 0 \\
0 & -
\end{array
\right)
\end{equation
and
\begin{eqnarray}
\alpha  &=&\frac{V_{s}+V_{as}}{2}  \notag \\
\beta  &=&\frac{V_{s}+V_{as}}{2} \\
\gamma  &=&\frac{V_{s}-V_{as}}{2}  \notag
\end{eqnarray
in which $V_{s}=e^{-2r}$ and $V_{as}=e^{2r}+\delta $ are the squeezing and
anti-squeezing fluctuation of the initial squeezing state, respectively. The
existence of $\delta $ is caused by classical and uncorrelated noise in
quantum resource and leads to the impurity of quantum resource. To partially
reconstruct all relevant entries of its associated covariance matrix we
perform $6$ different measurements on the output optical modes. These
measurements include the amplitude and phase quadratures of the output
optical modes $\bigtriangleup ^{2}\hat{x}_{A},\bigtriangleup ^{2}\hat{p
_{A},\bigtriangleup ^{2}\hat{x}_{B},\bigtriangleup ^{2}\hat{p}_{B}$, and the
cross correlations $\Delta ^{2}\left( \hat{x}_{A}+\hat{x}_{B}\right) $ and
\Delta ^{2}\left( \hat{p}_{A}-\hat{p}_{B}\right) $. The covariance elements
are calculated via the identities \cite{recontructcm
\begin{align}
Cov\left( \hat{R}_{i},\hat{R}_{j}\right) & =\frac{1}{2}\left[ \Delta
^{2}\left( \hat{R}_{i}+\hat{R}_{j}\right) -\Delta ^{2}\hat{R}_{i}-\Delta ^{2
\hat{R}_{j}\right] ,  \notag \\
Cov\left( \hat{R}_{i},\hat{R}_{j}\right) & =-\frac{1}{2}\left[ \Delta
^{2}\left( \hat{R}_{i}-\hat{R}_{j}\right) -\Delta ^{2}\hat{R}_{i}-\Delta ^{2
\hat{R}_{j}\right] .
\end{align
As an example, the partially reconstructed covariance matrix of the
experimentally prepared pure TMSS is
\begin{equation}
\sigma _{AB}=\left(
\begin{array}{cccc}
1.26 & 0 & -0.79 & 0 \\
0 & 1.28 & 0 & 0.80 \\
-0.79 & 0 & 1.32 & 0 \\
0 & 0.80 & 0 & 1.2
\end{array
\right) .
\end{equation
With these theoretical and experimentally partially reconstructed covariance
matrix, we can quantify the Gaussian EPR steering of quantum states.

The steerability of Bob by Alice ($A\rightarrow B$) for a two-mode Gaussian
state with Gaussian measurement can be quantified by \cite{criterion}
\begin{equation}
\mathcal{G}^{A\rightarrow B}(\sigma _{AB})=\max \bigg\{0,\frac{1}{2}\ln
\frac{\det \mathbf{A}}{\det \sigma _{AB}}\bigg\}.  \label{eqn:parameter}
\end{equation
The quantity $\mathcal{G}^{A\rightarrow B}$ vanishes iff the state described
by $\sigma _{AB}$ is nonsteerable by Gaussian measurements. The steerability
of Alice by Bob ($B\rightarrow A$), quantified by $\mathcal{G}^{B\rightarrow
A}$, can be obtained by replacing the role of A by B.

It is worth mentioning that the first criterion for the demonstration of the
EPR paradox was proposed by M. Reid in 1989 \cite{ReidPRA}, and it was
linked to the concept of steering lately \cite{PRA80}. For a two-mode
Gaussian state, the quantifier as given by Eq. (4) is equivalent to the
Reid's criterion which is based on Heisenberg uncertainty relation, i.e., a
state is steerable from Alice to Bob if the following relation is violated:
\begin{equation}
V_{X_{B}|X_{A}}V_{P_{B}|P_{A}}\geq 1,
\end{equation
where $V_{X_{B}|X_{A}}$ and $V_{P_{B}|P_{A}}$ are the conditional variances
of Bob's measurements conditioned on Alice's results, and
V_{X_{B}|X_{A}}V_{P_{B}|P_{A}}=\det \sigma _{AB}/\det \mathbf{A}$. For
steering from Bob to Alice, the roles of A and B are swapped. This criterion
does not require the assumption of Gaussian states or Gaussian systems to be
valid.

The PPT criterion \cite{ppt} is applied to characterise the entanglement,
which is a sufficient and necessary condition for a Gaussian TMSS. The PPT
value can be determined by
\begin{equation}
\sqrt{\frac{\Gamma -\sqrt{\Gamma ^{2}-4\det \sigma _{AB}}}{2}},
\end{equation
where $\Gamma =\det \mathbf{A}+\det \mathbf{B}-2\det \mathbf{C}$. When the
PPT value is less than $1$, the two modes are entangled.

At first, we investigate the effect of purity of quantum state on EPR
steering. The purity of a quantum state $\varrho $ is $\mu \left( \varrho
\right) =\mathbf{Tr}\varrho ^{2}$ \cite{Adesso}, which varies in the range
of $\frac{1}{D}\leq \mu \leq 1$, where $D=\dim \mathcal{H}$ for a given
Hilbert space $\mathcal{H}$. The minimum is reached by the totally random
mixture, the upper bound is saturated by pure states. In the limit of CV
systems ($D\rightarrow \infty $ ), the minimum purity tends to zero. The
purity of any N-mode Gaussian state is only function of the symplectic
eigenvalues $\upsilon _{i}$ of covariance matrix $\sigma$ \cite{purity2},
which is given by
\begin{equation}
\mu \left( \varrho \right) =\frac{1}{\prod_{i}\upsilon _{i}}=\frac{1}{\sqrt
\det \sigma }}.
\end{equation
Regardless of the number of modes, the purity of a Gaussian state is fully
determined by the global symplectic invariant $\det\sigma $ alone \cit
{purity2}.

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=80mm]{Fig2.eps}
\end{center}
\caption{The influence of purity of TMSS on the Gaussian EPR steering in a
lossy channel. \textbf{a,} Steering parameter quantified by Eq. (4). \textbf
b,} Steering certification according to Eq. (5). The region under the bound
of $1$ (black dotted line) indicates the existence of steerability.
Theoretical results for different purities of $1$, $0.8$, $0.6$, $0.4$ and
0.2$ are considered, which correspond to different $\protect\delta $ of
0,0.5,1.33,2.99$ and $7.97$, respectively. The experimental data (black
circles and squares) for the pure TMSS agree well with the theoretical
prediction (black solid and dash curves). Error bars of experimental data
represent $\pm 1$ standard deviation.}
\end{figure}

In practice, the loss and noise in the generation system will lead to
generation of impure quantum states \cite{esd5,PRA79}. For example, the
phonon noise in the three-color entangled state prepared by a non-degenerate
optical parametric oscillator (OPO) is a type of classical and uncorrelated
noise which leads to an impure three-color entangled state. Scattered light
by thermal phonons inside a second-order nonlinear crystal is the source of
additional phase noise observed in OPOs \cite{esd5,PRA79}. For our TMSS, the
relationship between purity $\mu \left( \varrho \right) $ and $\delta $ is
\begin{equation}
\mu \left( \varrho \right) =\frac{1}{\sqrt{\det \sigma _{AB}}}=\frac{1}
V_{s}V_{as}}=\frac{e^{2r}}{e^{2r}+\delta }.
\end{equation

The scheme of the transmission of Gaussian EPR steering in a lossy channel
for pure and impure TMSSs is shown in Fig. 1a. Mode $\hat{B}$ of the TMSS is
distributed in a lossy channel, the mode after transmission is given by
\hat{B}_{L}=\sqrt{\eta _{L}}\hat{B}+\sqrt{1-\eta _{L}}\hat{v}$, where $\eta
_{L}$ and $\hat{v}$ represent the transmission efficiency of quantum channel
and the vacuum mode induced by loss into the quantum channel, respectively.
The partially reconstructed covariance matrix of the TMSS distributed over
lossy channel in the experiment we hav
\begin{eqnarray}
\alpha _{L} &=&\frac{V_{s}+V_{as}}{2}  \notag \\
\beta _{L} &=&\frac{\eta _{L}\left( V_{s}+V_{as}\right) }{2}+1-\eta _{L} \\
\gamma _{L} &=&\frac{\sqrt{\eta _{L}}\left( V_{s}-V_{as}\right) }{2}  \notag
\end{eqnarray}

Secondly, the transmission of Gaussian EPR steering in a noisy channel is
investigated, as shown in Fig. 1b, where the excess noise higher than the
vacuum noise exists. Mode $\hat{B}$ of the TMSS is distributed in the noisy
channel, the mode after transmission is given by $\hat{B}_{N}=\sqrt{\eta _{N
}\hat{B}+\sqrt{\left( 1-\eta _{N}\right) g}\hat{N}+\sqrt{1-\eta _{N}}\hat{v}
, where $\hat{N}$ and $g$ are the Gaussian noise in the channel with
transmission efficiency of $\eta _{N}$ and the magnitude of noise,
respectively. The partially reconstructed covariance matrix of the TMSS
distributed over noisy channel is
\begin{eqnarray}
\alpha _{N} &=&\frac{V_{s}+V_{as}}{2}  \notag \\
\beta _{N} &=&\frac{\eta _{N}\left( V_{s}+V_{as}\right) }{2}+\left( 1-\eta
_{N}\right) \left( g+1\right)  \\
\gamma _{N} &=&\frac{\sqrt{\eta _{N}}\left( V_{s}-V_{as}\right) }{2}  \notag
\end{eqnarray}

Finally, we demonstrate the revival of Gaussian EPR steering by establishing
a correlated noisy channel, where an auxiliary channel carrying on the noise
correlated with that of the noisy channel is introduced. In today's
communication system, it is relevant to consider channels with correlated
noise (i.e. non-Markovian environment), because the time and space
correlated noise exist naturally \ \cite{correlated noise1,correlated
noise2,correlated noise3}. Furthermore, it has been experimentally found
that correlated noisy channel can be established by two bounded fibers \cit
{bound fiber}, thus the method we use in experiment to built the
non-Markovian environment is applicable and practical significant. As shown
in Fig. 1b. The auxiliary noisy channel based on an ancillary coherent state
$\hat{c}_{an}$, which is modulated to carry the noise information correlated
with that in the noisy quantum channel, saying $\hat{c}_{an}^{\prime }=\hat{
}_{an}+\sqrt{g_{an}}\hat{N}$ ($g_{an}$ is the magnitude of the correlated
noise and is experimentally adjustable). Thus there exists correlation
between the quantum system and environment. The noisy transmitted mode $\hat
B}_{N}$ couples with the auxiliary noisy channel on a revival beam-splitter
with transmission efficiency $T$. One of the output mode we need from the
revival beam-splitter is $\hat{B}_{R}=\sqrt{\eta _{N}T}\hat{B}+\sqrt{\left(
1-\eta _{N}\right) gT}\hat{N}-\sqrt{\left( 1-T\right) g_{an}}\hat{N}+\sqrt
\left( 1-\eta _{N}\right) T}\hat{v}-\sqrt{1-T}\hat{c}_{an}$, when the
adjustable $g_{an}$ and $T$ satisfy
\begin{equation}
\frac{g}{g_{an}}=\frac{1-T}{\left( 1-\eta _{N}\right) T},  \label{condition}
\end{equation
it can be described as
\begin{equation}
\hat{B}_{R}=\sqrt{\eta _{N}T}\hat{B}+\sqrt{\left( 1-\eta _{N}\right) T}\hat{
}-\sqrt{1-T}\hat{c}_{an}.
\end{equation
Thus the excess noise in mode $\hat{B}_{R}$ can be removed totally, and the
steerability between $\hat{B}_{R}$ and $\hat{A}$ can be revived to some
degree. After the revival operation, the covariance matrix of the revived
state is
\begin{eqnarray}
\alpha _{R} &=&\frac{V_{s}+V_{as}}{2}  \notag \\
\beta _{R} &=&\frac{\eta _{N}T\left( V_{s}+V_{as}\right) }{2}+1-\eta _{N}T \\
\gamma _{R} &=&\frac{\sqrt{\eta _{N}T}\left( V_{s}-V_{as}\right) }{2}  \notag
\end{eqnarray}

\textbf{Experimental set-up}

The experimental set-up is shown in Fig. 1c. We use a non-degenerate optical
parametric amplifier (NOPA) to prepare the TMSS. NOPA is pumped by a
continuous-wave intracavity frequency-doubled Nd:YAP-LBO laser with two
output wavelengths at 540 nm and 1080 nm. The NOPA consists of an $\alpha
-cut type-II KTP crystal which is front face coated (work as the input
coupler with transmittance of 21.2\% and 0.04\% at 540 nm and 1080 nm,
respectively) and a concave mirror (work as the output coupler with
transmittance of 0.5\% and 12.5\% at 540 nm and 1080 nm, respectively). Via
the frequency down-conversion process of the pump field at 540 nm inside the
NOPA which operated at deamplification condition [the relative phase between
the pump laser and the injected signal equals to $\left( 2n+1\right) \pi $],
the TMSS at 1080 nm we want can be generated. Experimentally, the TMSS we
prepared is a nearly pure -3.0 dB state at sideband frequency of 3 MHz.

The loss in the quantum channel is experimentally mimicked by a
beam-splitter, which is composed by a half-wave plate and a polarization
beam-splitters (PBS), the vacuum mode couples into system from the other
input port of the beam-splitter. In the noisy environment, to couple the
transmitted mode $\hat{B}$ and the noise-modulated coherent beam (auxiliary
noisy channel $\hat{c}_{an}^{\prime }$), the beam-splitter (revival
beam-splitter) is composed by a PBS, a half-wave plate and another PBS as
shown in Fig. 1c, where the transmission efficiency $\eta _{N}$ and $T$ are
adjusted by setting the half-wave plate. The correlated noises between the
quantum system and environment are generated by same signal generator
experimentally. We reconstruct the covariance matrix of the quantum state
with the homodyne detection systems, and quantify its steering properties.

\section*{RESULTS}

\textbf{The effect of purity on EPR steering}

\begin{figure*}[tbp]
\begin{center}
\includegraphics[width=145mm]{Fig3.eps}\vspace*{-0.3cm}
\end{center}
\caption{Disappearance and revival of EPR steering in noisy channel. \textbf
a,b} Steering parameters quantified by Eq. (4) and (5) in a noisy channel
where the variance of the excess noise is taken as $5$ times of vacuum
noise, respectively. Black dash and solid curves are the theoretical
predictions of $\mathcal{G}^{B_{N}\rightarrow A}$ and $\mathcal{G
^{A\rightarrow B_{N}}$, respectively. Region I, two-way steering; region II,
one-way steering; region III, no steering. Red dash and solid curves are the
theoretical predictions of $\mathcal{G}^{B_{R}\rightarrow A}$ and $\mathcal{
}^{A\rightarrow B_{R}}$ after steering revival, respectively. The Blue dash
and solid curves are the theoretical predictions of $\mathcal{G
^{B_{L}\rightarrow A}$ and $\mathcal{G}^{A\rightarrow B_{L}}$ in a pure
lossy but noiseless channel. \textbf{c,d} Steering parameters quantified by
Eqs. (4) and (5) in a transmission efficiency fixed ($\protect\eta _{N}=0.6
) noisy channel with different noise levels. Black dash and solid curves are
theoretical predictions of $\mathcal{G}^{B_{N}\rightarrow A}$ and $\mathcal{
}^{A\rightarrow B_{N}}$ as a function of the excess noise in the unit of
shot noise level. Red dash and solid curves are the theoretical predictions
of $\mathcal{G}^{B_{R}\rightarrow A}$ and $\mathcal{G}^{A\rightarrow B_{R}}$
after steering revival, respectively. The inset figure shows the detailed
theoretical predictions of $\mathcal{G}^{B_{R}\rightarrow A}$ and $\mathcal{
}^{A\rightarrow B_{R}}$ when the excess noise is relatively low. Error bars
of experimental data represent $\pm 1$ standard deviation and are obtained
based on the statistics of the measured data.}
\end{figure*}

Figure 2 shows the effect of purity of the TMSS on EPR steering in a lossy
channel, where the steerabilities quantified by two criteria given by Eqs.
(4) and (5) are shown in Fig. 2a and 2b, respectively. For a pure $-3$ dB
TMSS, one-way steering appears in a lossy channel, and the critical point of
one-way steering property occurs at $\eta _{L}=0.5$. For comparison, four
different purities of the TMSSs at $0.8$, $0.6$, $0.4$ and $0.2$ are
considered, which show that the maximum steerability and the phenomenon of
one-way EPR steering property happens earlier with the decreasing of purity.
This confirms that the impurity of the TMSS leads to decrease of
steerabilities and two-way steering region. However, the steerability from
Alice to Bob never disappear until the purity tends to zero, which means
that the impurity of the TMSS does not lead to sudden death of Gaussian EPR
steering. The results in Fig. 2a and 2b confirm that two criteria given by
Eqs. (4) and (5) are equivalent for a TMSS with Gaussian measurement in a
lossy channel.

\begin{figure*}[tbp]
\begin{center}
\includegraphics[width=180mm]{Fig4.eps}
\end{center}
\caption{Results for the effect of purity of the TMSS on EPR steering in a
noisy channel. \textbf{a, b,} The dependence of Gaussian EPR steering on
transmission efficiency and excess noise with initial pure and impure
(purity is 0.8) TMSS, respectively. The steerabilities of $\mathcal{G
^{A\rightarrow B_{N}}$ and $\mathcal{G}^{B_{N}\rightarrow A}$ are given by
the light red and blue colors, respectively. \textbf{c,} The effect of noise
on EPR steering of the TMSS with purities of 1 [blue curves] and 0.8 [red
curves], respectively. The excess noise is taken as $2$ times of vacuum
noise.}
\end{figure*}

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=56mm]{Fig5.eps}
\end{center}
\caption{The region plot parameterized by transmission efficiency and excess
noise. The region below the black curve, which is the bound of entanglement
sudden death, represents the existence of the entanglement. The region below
the blue and red curves, which are the bound of the disappearance of
\mathcal{G}^{A\rightarrow B}$ and $\mathcal{G}^{B\rightarrow A}$
respectively, represents the existence of EPR steering. }
\end{figure}

In the experiment, to determine the error bar of the experimental results,
we measure three covariance matrices for each quantum state and obtain the
corresponding EPR steering parameters. By calculating the mean value and the
standard deviation of these three EPR steering parameters, error bars are
obtained. The errors of the experimental data mainly come from the phase
fluctuation of the phase locking system, and we take much efforts to
suppress phase fluctuation of each phase locking system to be $\theta
\approx 1.5^{\circ }$ by adjusting locking parameters.

\textbf{The effect of noise on EPR steering}

As shown in Fig. 3a and 3b, sudden death of Gaussian EPR steering of the
TMSS in a noisy channel is observed, where the variance of the excess noise
is taken as $5$ times of vacuum noise. Specifically, in region I (
0.96<\eta _{N}\leq 1$), the two-way steering exists, both Alice and Bob can
steer with each other. With the decreasing of $\eta _{N}$, when $0.86<\eta
_{N}\leq 0.96$ (labeled as region II), one-way steering appears, i.e. Bob
can steer Alice's state while Alice cannot steer Bob's state. More
seriously, when $\eta _{N}\leq 0.86$, the steerability between Alice and Bob
totally disappeared (labeled as region III), i.e. the sudden death of the
Gaussian EPR steering happens.

We also demonstrate sudden death of Gaussian EPR steering with different
excess noise levels (noise levels are taken as $5,10,15,20$ times of vacuum
noise) at fixed transmission efficiency ($\eta _{N}=0.6$) in Fig. 3c and 3d.
We can see the two-way steering (labeled as region I) and/or one-way
steering (labeled as region II) could still exist when the excess noise is
relatively low. When the variance of the excess noise is higher than 0.47
(labeled as region III), the sudden death of the Gaussian EPR steering
happens. The results in Fig. 3 confirm that two criteria given by Eqs. (4)
and (5) are equivalent for the sudden death and revival of Gaussian EPR
steering of a TMSS with Gaussian measurement.

The effect of purity of initial TMSS on EPR steering in a noisy environment
is also investigated as shown in Fig. 4. Figure 4a and 4b are the
three-dimensional theoretical results of Gaussian EPR steering properties of
initial TMSS with different purities of $1$ and $0.8$, respectively. It is
obvious that the direction of EPR steering is changed with different noise
levels, which is the same with the result in Ref. \cite{PRAQin}. We can see
that the blue surface ($\mathcal{G}^{B_{N}\rightarrow A}$) and light red
surface ($\mathcal{G}^{A\rightarrow B_{N}}$) have crossover, which happens
with excess noise of $0.24$ and $0.49$ for the pure and impure TMSSs,
respectively. Figure 4c shows the effect of purity of the TMSS on EPR
steering in a noisy channel with excess noise of $2$ times of vacuum noise.
It's obvious that the sudden death of EPR steering for the impure state
happens earlier and the maximum steerability is smaller than that of the
pure state.

\noindent\ \ \textbf{The revival of EPR steering}

To revive the disappeared EPR steering, we fix $T=0.9$ for the revival
beam-splitter and adjust the magnitude of correlated noise $g_{an}$. The
parameters $g/g_{an}$ are set to $0.14,0.18,0.28,$ and $0.56$ for different
quantum channel efficiencies of $0.2,0.4,0.6,$ and $0.8$, respectively,
according to Eq. (8). After the revival operation, the disappeared Gaussian
EPR steering is revived as shown in Fig. 3a and 3b. Because the excess noise
in mode $\hat{B}_{R}$ is removed totally in the revival operation, the
Gaussian EPR steering is revived even at the noise level of $20$ times of
vacuum noise , as shown in Fig. 3c and 3d. Generally, $T$ is better to be
close to $1$, because it results in adding an extra linear loss on the
revived mode, this is why the steerabilities between mode $\hat{B}_{R}$ and
\hat{A}$ are not good as that between mode $\hat{B}_{L}$ and $\hat{A}$ [the
pure lossy transmitted regime, presented by blue dash and solid curves in
Fig. 3a].
\begin{figure}[tbp]
\begin{center}
\includegraphics[width=80mm]{Fig6.eps}
\end{center}
\caption{Results for distillation of Gaussian steering with a correlated
noisy channel and measurement-based NLA. Black solid and dash curves are the
theoretical predicted steerabilities in a noisy channel with excess noise 5
times of vacuum noise. Red solid and dash curves are the predicted
steerabilities after revival of Gaussian steering with a correlated noisy
channel. The orange and purple curves are the predicted steerabilities when
the measurement-based NLA based on Bob's measurement results
with gain of 1.2 is applied before and after the revival of Gaussian
steering, respectively. The Blue solid and dash curves are the predicted
steerabilities in a pure lossy but noiseless channel. }
\end{figure}

\section*{DISCUSSION}

Note that the noise in quantum channel and the noise in initial TMSS have
different influences on EPR steerability. With the decreasing of purity of
initial TMSS, the maximum steerability decreases and the phenomenon of
one-way EPR steering property happens earlier than that of\ the pure state
(between $0.5\leq \eta \leq 1$), but the death of the EPR steering never
happens. While with the increasing of noise in the quantum channel, the
two-way EPR steering firstly declines to one-way EPR steering, and the
sudden death of the Gaussian EPR steering finally occurs. Thus we can
distinguish where the excess noise exists (whether in the communication
channel or in the initial TMSS) by the Gaussian EPR steerability easily.

The sudden death of Gaussian EPR steering is different from that of the
Gaussian entanglement, which is characterised by the positive partial
transposition (PPT) criterion \cite{ppt} (see Methods for details). Figure 5
shows the entanglement and steering region parameterized by transmission
efficiency and excess noise, where the purple region represents the two-way
steering region, the cyan region represents the one-way steering region of
\mathcal{G}^{A\rightarrow B}$, the light red region represents the one-way
steering region of $\mathcal{G}^{B\rightarrow A}$, and the yellow region
represents the entanglement without steering region, respectively. Firstly,
the direction of the EPR steering is changed with the increase of excess
noise which is not the case in entanglement sudden death. The sudden death
of EPR steering and entanglement happens when there is excess noise in the
quantum channel. For a pure $-3$ dB TMSS, the direction of EPR steering is
changed when the excess noise is 0.25 times of vacuum noise. This result
comes from the intrinsic asymmetric character of EPR steering. Secondly, it
is obvious that the sudden death of Gaussian EPR steering happens earlier
than that of the entanglement with same excess noise level. Thirdly,
different from entanglement, the sudden death of Gaussian EPR steering in a
noisy channel is directional.

Recovering EPR steering with correlated noisy channel can only remove the
effect of excess noise in the transmission channel, the effect of loss in
the transmission channel cannot be removed. Some other technologies such as
noiseless linear amplification can be used to eliminate the effect of loss
\cite{linear amp1,linear amp2,linear amp3,linear amp4}. Recently, it has
been shown that the measurement-based noiseless linear amplification can be
used to distill Gaussian EPR steering in both lossy and noisy channels,
which provides another feasible way to recover Gaussian EPR steering in
quantum channel \cite{arxivLiu}.

Comparing the present work of EPR steering revival with a correlated noisy
channel and distillation of EPR steering with measurement-based NLA, the
differences are as following. Firstly, the steering revival with a
correlated noisy channel can only remove the effect of excess noise in the
quantum channel, the effect of loss in the transmitted channel cannot be
removed. But the distillation of EPR steering works for both lossy and noisy
channels. Secondly, the revival of Gaussian EPR steering with a correlated
noisy channel can remove the excess noise in the quantum channel totally,
but distillation of EPR steering with measurement-based NLA cannot remove
the excess noise in all range of transmission efficiency [orange curves in
Fig. 6]. Thirdly, the revived steerability with a correlated noisy channel
cannot exceed the initial steerability (that in a pure lossy channel), but
the distilled steerability can exceed the initial steerability.

Combining the methods of correlated noisy channel and the NLA may provide an
option to distill Gaussian EPR steering in lossy and noisy environment. As
shown in Fig. 6, if the measurement-based NLA based on Bob's
measurement results with gain of 1.2 is applied after the revival of
Gaussian EPR steering with a correlated noisy channel, the revived
steerability could be enhanced by using the measurement-based NLA. In this
case, the revival of steering with a correlated noisy channel removes the
effect of excess noise in quantum channel. Followed by the measurement-based
NLA, the effect of loss in quantum channel can be eliminated and the
distillation of Gaussian steering can be obtained.

Our work demonstrates evolution of bipartite Gaussian EPR steering of a TMSS
transmitted in lossy and noisy quantum communication channels. We show the
noise in quantum channel and the noise in the initial state result in
different influences on EPR steerability. The impurity of the state only
decrease the steerabilities and two-way steering region, but never leads to
sudden death of EPR steering. While the excess noise in the quantum channel
leads to sudden death of Gaussian EPR steering. We successfully revive EPR
steering after the sudden death occurs in noisy channel by establishing a
correlated noisy channel. The presented results provide useful references
for understanding the steering properties of quantum state which interacts
with different environments, also have potential applications in asymmetric
quantum information processing exploiting EPR steering as a valuable
resource.

\textbf{ACKNOWLEDGMENTS}

This research was supported by the NSFC (Grant Nos. 11904160, 11834010 and
62005149), the program of Youth Sanjin Scholar, National Key R\&D Program of
China (Grant No. 2016YFA0301402), and the Fund for Shanxi \textquotedblleft
1331 Project\textquotedblright\ Key Subjects Construction.

X. D. and Y. L. contributed equally to this work.