# A First Analysis of the Stability of Takens’ Embedding

Han Lun Yap, Armin Eftekhari, Michael B. Wakin, and Christopher J. Rozell

##  I Introduction

Data in the form of time series pervades many areas of science and
engineering. Traditionally, one might impose a suitable probability model for
the data, yielding tools such as the ARMA model and the Kalman filter [1].
More recently, modeling the time-series data as a one-dimensional measurement
of the possibly high-dimensional states of a nonlinear dynamical system has
also produced many positive results [2, 3].

The theoretical foundation of time-series analysis based on the dynamical
system model hinges on a seminal paper by Takens [4], whose ideas were later
expanded by Sauer et al. [5]. The surprising result of the paper asserts that
when the states of a dynamical system are confined to a low-dimensional
attractor, complete information about the hidden states can be preserved in
the time-series output (see Figure 1). Indeed, many systems of interest do
have this type of structure [6], and a variety of algorithms for tasks such as
time-series prediction exploit Takens’ result [2].

To be precise, suppose we have a dynamical system whose system states $x(t)$
lie on a submanifold $\mathcal{A}\subset\mathbb{R}^{N}$. Given a sampling
interval $T_{s}$, one could define the discretized dynamics of the system
through the flow function $\phi_{T_{s}}:\mathcal{A}\rightarrow\mathcal{A}$
defined as $x(t+T_{s})=\phi_{T_{s}}(x(t))$. The time series $s(t):=h(x(t))$ is
modeled as a one-dimensional observation of the system dynamics through the
observation function $h$. Usually the sampling time $T_{s}$ is fixed and we
will drop the subscript $T_{s}$ from $\phi_{T_{s}}$ to ease notations. Takens
[4] defined the _delay coordinate map_
$F:\mathbb{R}^{N}\rightarrow\mathbb{R}^{M}$ as a mapping of the state vector
$x(t)\in\mathbb{R}^{N}$ of a dynamical system to a point
$F(x(t))\in\mathbb{R}^{M}$ in the _reconstruction space_ by taking $M$
uniformly spaced samples of the time series $s(t)$ and concatenating them into
a single vector, i.e.,
$F(x(t))=[s(t),\;s(t-T_{s}),\;\cdots,\;s(t-(M-1)T_{s})]^{T}$. Takens’ result
states that for _almost every_ smooth observation function $h(\cdot)$, the
delay coordinate map is an _embedding_ (i.e., a one-to-one immersion) of the
state space attractor with dimension $d$ when $M>2d$. As seen in Figure 1,
despite the state being hidden from direct observation, the topology of the
attractor that characterizes the dynamical system can be preserved in the
time-series data when it is arranged into a delay coordinate map.

Figure 1: Takens’ embedding of a state-space attractor $\mathcal{A}$ into
reconstruction space from time series measurements.

In the absence of measurement or system noise, Takens’ result indicates that
the delay coordinate vectors $F(x(t))$ are a proxy for the hidden system
states $x(t)$. However, in the presence of noise, a one-to-one mapping as
guaranteed by Takens’ result may not be sufficient to guarantee the robustness
of any processing performed in the reconstruction space. The underlying
problem is that while Takens’ theorem guarantees the preservation of the
attractor’s _topology_ , it does not guarantee that the _geometry_ of the
attractor is also preserved. To be precise, Takens’ result guarantees that two
points on the attractor do not map to the same point in the reconstruction
space, but there are no guarantees that close points on the attractor remain
close under this mapping (or far points remain far). Consequently, relatively
small imperfections could have arbitrarily large, unwanted effects when the
delay coordinate map is used in applications.

In the signal processing community, recent work has highlighted the importance
of well-conditioned measurement operators to ensure the geometry of a low-
dimensional signal family is preserved. The notion of a _stable embedding_ by
a measurement operator $F$ introduced in compressed sensing (CS) [7] has been
extended to manifold signal families [8, 9, 10, 11]. The results in [8]
guarantee that distances between points on a $D$-dimensional submanifold are
approximately preserved when $F$ is a random, linear orthoprojector and
$M=O(d\log(N))$.111The notation $O(\cdot)$ simply means “on the order of”.
This result is extended in [10, 11] to include other random stable embedding
operators. Moreover, the results in [9, 11] removed the logarithmic dependence
on $M$. Besides signal recovery [12], the stable embedding of a signal family
also ensures that data processing and inference algorithms can perform with
approximately the same guarantees in the measurement space as in the ambient
space [13].

The conditions for Takens’ Embedding Theorem to hold ignore the effects of
noise. As such, time-series analysis in practice requires a careful empirical
determination of the sampling time $T_{s}$ and number of delays $M$, typically
via looking at the first null of mutual information or autocorrelation between
two consecutive samples of the time series [14, 2]. This usually results in
using a number of delay coordinates $M$ larger than the minimum prescribed by
Takens’ theorem [15].

In this paper, we use tools and ideas in CS to provide a first theoretical
justification for the choice of $M$ in noisy conditions. In fact, we show that
under certain conditions on the dynamical system, measurement function $h$,
number of delays $M$ and sampling time $T_{s}$, the delay-coordinate map can
be a stable embedding of the dynamical system’s attractor.

##  II Main Results

For the rest of this section, consider a dynamical system given by the flow
$\phi=\phi_{T_{s}}$ with sampling time $T_{s}$ and whose states are lying on a
Riemannian submanifold $\mathcal{A}$ of $\mathbb{R}^{N}$. We shall also assume
that $\phi$ is a smooth function. Before presenting the stable embedding
results proper, we will need to introduce an important characterization of
manifolds that will be useful for describing certain local and global
properties of the system manifold $\mathcal{A}$.

######  Definition II.1.

Let $\mathcal{A}$ be a compact Riemannian submanifold of $\mathbb{R}^{N}$. The
_reach_ $\tau$ of $\mathcal{A}$ is defined as the largest positive real number
having the following property: The open normal bundle about $\mathcal{A}$ of
radius $a$ is embedded in $\mathbb{R}^{N}$ for all $a<\tau$.

Reach is an equivalent measure to the condition number used in [16, 10, 8].
See the above cites and [11, 17] for more on reach and its implications.

On top of the properties of the system manifold, our results also depend on
the types of measurement functions used to form the time series data. First,
consider a set of basis functions $h_{p}:\mathbb{R}^{N}\rightarrow\mathbb{R}$
for $p=1,\cdots,P$ that will span a $P$-dimensional subspace of all functions
from $\mathbb{R}^{N}$ to $\mathbb{R}$. An element chosen from this class of
measurement functions can be denoted by $h_{\alpha}$ for some
$\alpha\in\mathbb{R}^{P}$ defined as
$h_{\alpha}:=\sum_{p=1}^{P}\alpha_{p}h_{p}$. Next, we define the delay
coordinate map with measurement function $h_{\alpha}$ as
$F_{\alpha}(x):=\sum_{p=1}^{P}\alpha_{p}F_{p}(x)=\left[h_{\alpha}(x),\cdots,h_{\alpha}\left(\phi^{-M+1}\left(x\right)\right)\right]^{T}$,
where $F_{p}$ is defined as the delay-coordinate vector formed with the $p$-th
basis function
$F_{p}(x):=\left[h_{p}(x),\cdots,h_{p}\left(\phi^{-M+1}\left(x\right)\right)\right]^{T}$.
In this paper, we are considering the class of linear functions. Thus, we take
$P=N$ (where $N$ is the ambient dimension) and
$h_{p}(\cdot):=\langle\cdot,e_{p}\rangle$ where $e_{p}$ is the $p$-th
canonical vector.

Additionally, we need a metric to describe the dynamical system in question.
For any state $x\in\mathcal{A}$, we define the $M\times N$ matrix $G_{x}$ as
$G_{x}:=\left(F_{1}(x)\;|\;\cdots\;|\;F_{N}(x)\right),$ and define the
trajectory vector $\widetilde{g}\left(x\right)\in\mathbb{R}^{MN}$ starting at
point $x\in\mathcal{A}$ as
$\widetilde{g}\left(x\right):=\left[x,\phi^{-1}\left(x\right),\cdots,\phi^{-M+1}\left(x\right)\right]^{T}$.
Observe that $G_{x}$ is basically the matrix version of
$\widetilde{g}\left(x\right)$ and that
$F_{\alpha}(x)=\sum_{i=1}^{N}\alpha_{i}F_{i}(x)=G_{x}\alpha$. We also define
the trajectory manifold
$\widetilde{g}\left(\mathcal{A}\right)\subset\mathbb{R}^{MN}$ as the set of
all trajectories in $\mathcal{A}$, i.e.,
$\widetilde{g}\left(\mathcal{A}\right):=\left\\{\widetilde{g}\left(x\right)\,:\,x\in\mathcal{A}\right\\}$.
Notice that since the first $N$ entries of $\widetilde{g}\left(x\right)$ are
distinct for all
$\widetilde{g}\left(x\right)\in\widetilde{g}\left(\mathcal{A}\right)$ and
$\phi$ is a smooth function, $\widetilde{g}\left(\mathcal{A}\right)$ is
diffeomorphic to $\mathcal{A}$. Then for any $G\in\mathbb{R}^{M\times N}$,
define $r(G):=\frac{\|G\|_{F}^{2}}{\|G\|_{2}^{2}}$ which we shall call the
_soft-rank_ of the matrix $G$. In the stable embedding results that follow, we
shall see that the soft-rank of the matrices $G_{x}-G_{y}$ for all pairs
$x,y\in\mathcal{A}$ capture how well information about the dynamical system
lying on $\mathcal{A}$ is captured by the basis measurement functions
$\\{h_{p}\\}$. The soft-rank is a refined metric to replace the rank of a
matrix.

Under the assumptions that there are no periodic orbits of period less than
$M$ on $\mathcal{A}$ and $M\leq N$, it is easy to see that
$\operatorname{rank}\left(G_{x}-G_{y}\right)=M$ for all pairs of
$x,y\in\mathcal{A}$. Having full rank only says that the matrix $G_{x}-G_{y}$
has $M$ non-zero singular values, but it does not say anything about the
distribution of these singular values. On the other hand, the soft-rank is the
ratio of the $\ell_{2}$ and $\ell_{\infty}$ norms of the vector of singular
values, therefore $r(G_{x}-G_{y})=cM$ for some $c\in(\frac{1}{M},\;1]$ if
$G_{x}-G_{y}$ is full-rank. The value of $c$ tells us about the energy
distribution of the singular values, with $c$ tending to 1 when all the
singular values have almost the same values and tending to $\frac{1}{M}$ when
only one of the singular value is significant while the rest are close to
zero. We remark that via linear algebra, we have $c\approx 1$ whenever the
rows of $G_{x}-G_{y}$ have approximately equal length and are approximately
orthogonal.

With all the required notations now defined, we finally present our stable
embedding result for linear measurement functions.

######  Theorem II.1.

Suppose we have a dynamical system described by a flow $\phi$ with sampling
time $T_{s}$ and whose states lie on a $D_{\mathcal{A}}$-dimensional
submanifold $\mathcal{A}$ of $\mathbb{R}^{N}$. Also, suppose that the
trajectory manifold $\widetilde{g}\left(\mathcal{A}\right)$ has volume
$V_{\widetilde{g}\left(\mathcal{A}\right)}$ and reach
$\tau_{\widetilde{g}(\mathcal{A})}$. Let $\alpha\in\mathbb{R}^{N}$ be an
i.i.d. Rademacher $\pm 1$ sequence, and let $\delta$ be a predetermined stable
embedding conditioning. Define the infimum soft-rank
$r(\mathcal{A}):=\inf_{x\neq y\in\mathcal{A}}r(G_{x}-G_{y})$ and assume that
$V_{\widetilde{g}(\mathcal{A})}\geq
O\left(\tau_{\widetilde{g}(\mathcal{A})}^{D_{\mathcal{A}}}\right).$ If

$\displaystyle r(\mathcal{A})\geq$ $\displaystyle
O\left(\delta^{-2}D_{\mathcal{A}}\log\left(\frac{r(\mathcal{A})^{1/2}V_{\widetilde{g}(\mathcal{A})}^{1/D_{\mathcal{A}}}}{\tau_{\widetilde{g}(\mathcal{A})}}\right)\right),$
(1)

then except with a probability of at most
$\left(\frac{r(\mathcal{A})^{D_{\mathcal{A}}}V_{\widetilde{g}(\mathcal{A})}}{\tau_{\widetilde{g}(\mathcal{A})}^{D_{\mathcal{A}}}}\right)^{-C}$
with $C$ an absolute constant, the delay-coordinate map with $M$ delays
$F_{\alpha}$ is a stable embedding of $\widetilde{g}\left(\mathcal{A}\right)$,
i.e., for all $x,y\in\mathcal{A}$,

$\displaystyle(1-\delta)\leq\frac{\left\|F_{\alpha}(x)-F_{\alpha}(y)\right\|_{2}^{2}}{\left\|\widetilde{g}\left(x\right)-\widetilde{g}\left(y\right)\right\|_{2}^{2}}\leq(1+\delta).$
(2)

Due to space limitations, the proof of this theorem will be deferred to a
future publication. We note that the assumption
$V_{\widetilde{g}(\mathcal{A})}\geq
O\left(\tau_{\widetilde{g}(\mathcal{A})}^{D_{\mathcal{A}}}\right)$ is to
maintain a failure probability less than 1.

Our result describes how the geometry and the dynamics on the manifold coupled
with the number of measurements $M$, as seen through the soft-rank of
$G_{x}-G_{y}$ in (1), dictate the stable embedding conditioning of the delay-
coordinate map. Thus, this result is different from the typical Takens’
embedding result in two distinct ways. First, in the typical Takens’ embedding
result, only the topology of the manifold — basically the fact that the system
states lie on a manifold and the dimension of the manifold — is utilized.
Second, the typical result does not say what happens when we have more
measurements beyond the minimum required for an embedding.

The atypical condition (1) demands further elaboration. Suppose the number of
measurements $M$ is less than the ambient space dimension $N$ and there are no
periodic orbits of period less than $M$. First, notice that if the _infimum
soft-rank_ $r(\mathcal{A})$ scales with the number of measurements $M$, then
we will have a typical manifold stable embedding condition whereby a manifold
is stably embedded whenever $M$ scales with the dimension of the manifold and
logarithmically with other characterizations of the manifold (which is exactly
what appears on the right hand side of (1)), and the conditioning $\delta$ of
the embedding improves with increasing $M$. We remark that the improvement of
$\delta$ with $M$ matches what was observed by Casdagli et al. [15], namely
that with more measurements, distortion by the delay-coordinate map is
decreased. Now, for every $x,y\in\mathcal{A}$, the soft-rank $r(G_{x}-G_{y})$
indeed depends on the number of measurements $M$, but it also depends on the
sampling time $T_{s}$ coupled with various geometric properties of the
manifold and the dynamical system. As discussed previously, the soft-rank can
indeed be equal or close to $M$, but only when the rows of the matrix
$G_{x}-G_{y}$ (i.e., the chords $x-y,\phi^{-1}(x)-\phi^{-1}(y),\cdots$) have
almost the same length and are almost orthogonal. Moreover, these requirements
on $G_{x}-G_{y}$ have to hold for all $x,y\in\mathcal{A}$.

As expected, the length-preserving and orthogonality requirements described
above can be very stringent. For example, for a dynamical system whose inverse
flow $\phi^{-1}$ has a large maximal Lyapunov exponent, the requirement that
the chords retain almost the same length will be easily violated for any
decent sampling time $T_{s}$ and for $x,y$ close to one another.222 When the
inverse flow $\phi^{-1}$ has Lyapunov exponent $\lambda$, we have
$\|\phi^{-m}(x)-\phi^{-m}(y)\|_{2}\approx e^{m\lambda}\|x-y\|_{2}$ for $x,y$
close enough. Therefore, the chords going down the matrix $G_{x}-G_{y}$ grow
exponentially in length, severely violating the length-preserving requirement.
This is the effect of irrelevance as described in, e.g., [15]. Suppose $T_{s}$
is small and $x\approx y$. Then, the $M$ consecutive chords forming the $M$
rows of $G_{x}-G_{y}$ can more or less be approximated by vectors on a tangent
plane of $\mathcal{A}$ (say the tangent plane of $x$). Since
$\dim(\mathcal{A})=D_{\mathcal{A}}$ which could potentially be much less than
$M$, this means that the matrix $G_{x}-G_{y}$ will have only $D_{\mathcal{A}}$
significant singular values and thus $r(G_{x}-G_{y})=cM$ with $c\ll 1$. This
is the effect of redundancy. We will need many more measurements $M$ to “break
out” of the tangent plane to achieve $r(G_{x}-G_{y})=cM$ with a decent $c$
(i.e., $c\approx 1$). It is useful to add that if the manifold $\mathcal{A}$
indeed only lie on a low-dimensional subspace of $\mathbb{R}^{N}$, then the
infimum soft-rank will be _upper bounded_ by the dimension of this subspace
and thus will _not_ scale with $M$. This remark is also true when $M$ exceeds
the ambient dimension $N$, in which case the infimum soft-rank will plateau at
$N$. Notice that this imposes the limitation that $M\leq N$ for this result to
be useful which thus restricts the class of dynamical systems that we can
consider. Notably, this result is not adequate for dynamical systems residing
on low-dimensional ambient spaces ($N$ small), e.g., the Lorenz system.
Nonetheless, we intend to show in a future article that we can break this
ambient space dimension “barrier” by considering nonlinear measurement
functions.

Another important observation is that instead of embedding state-space vectors
$x\in\mathcal{A}$, we are in fact embedding trajectory vectors
$\widetilde{g}\left(x\right)\in\widetilde{g}\left(\mathcal{A}\right)\subset\mathbb{R}^{MN}$
as seen in (2). The variables on the right-hand side of (1) are also based on
the geometry of $\widetilde{g}\left(\mathcal{A}\right)$, e.g., the terms
$V_{\widetilde{g}\left(\mathcal{A}\right)}$ and
$\tau_{\widetilde{g}\left(\mathcal{A}\right)}$. Maintaining distance between
trajectories of the dynamical system in the reconstruction space may be
advantageous for some applications. We can translate this into a stable
embedding of the ambient space attractor $\mathcal{A}$ since there there is a
diffeomorphism (hence isomorphism) between the ambient and the trajectory
manifold. Understandably, this translation comes with a degradation of the
stable embedding conditioning $\delta$ previously obtained. While a detailed
mathematical treatment of this translation will be deferred to a future
article, it suffices to note here that the requirements for the minimum
degradation of the stable embedding conditioning caused by the translation
usually run contrary to requirements for the best infimum soft-rank in (1).
This is also somewhat similar to the tradeoff between redundancy and
irrelevancy in time series analysis. Redundancy happens when the sampling time
and/or the Lyapunov exponent (of the inverse flow) is small. In this case, the
trajectory and ambient vectors are scalar multiples of one another, which
means that there is minimum degradation of the stable embedding conditioning
going from trajectory to ambient space. However, the soft-rank of
$G_{x}-G_{y}$ suffers because now the rows of the matrix $G_{x}-G_{y}$ are
very similar to one another which violates the orthogonality condition (as
described before). However when irrelevance occurs (meaning that either the
sampling time or the Lyapunov exponent of the inverse flow becomes large), the
consecutive chords in $G_{x}-G_{y}$ rapidly decorrelate. When this occurs, the
soft-rank of $G_{x}-G_{y}$ may improve due to the diversity but the
conditioning is degraded when we pass from trajectory space to ambient space.

Lastly, we remark that our results mirror the notion that “almost every
measurement function provides an embedding” appearing in a typical Takens’
embedding statement (see [5]). This notion is represented by the random vector
of coefficients $\alpha\in\mathbb{R}^{N}$. Because $F_{\alpha}$ is a stable
embedding with high probability on the coefficients $\alpha$, this means that
most measurements functions in the space of functions defined by
$\\{h_{\alpha}=\alpha\;|\;\alpha\in\\{-1,+1\\}^{N}\\}$ can result in a stable
embedding by delay coordinate maps. In fact, the class of measurement
functions considered can be vastly expanded by replacing the Rademacher
sequence $\alpha$ with a Gaussian sequence instead. However due to the lack of
space, this result will be discussed in a future article.

##  III Experiment

We shall demonstrate our theoretical results with a simple example. Consider
the following discrete-time dynamical system with system states at time step
$n$, $x_{n}\in\mathbb{R}^{N}$, defined through $x_{n+1}=\Phi x_{n}$ where the
(linear) flow $\Phi\in\mathbb{R}^{N\times N}$ is a shift matrix, i.e.,

$\displaystyle\Phi=\left(\begin{array}[]{ccccc}0&1&0&\cdots&0\\\
\vdots&&\ddots&\ddots&0\\\ 0&&&0&1\\\ 1&0&\cdots&\cdots&0\end{array}\right).$

Also, suppose the initial condition of the system is given by
$x_{0}=[1,0\cdots,0]^{T}$. This system can be thought of as the dynamics of a
single point object translating down the entries of a $N$-dimensional vector.
Thus, the system states $x_{n}$ lie on a one-dimensional manifold
$\mathcal{A}$ that is parameterized by the object’s location on the vector. We
suppose that we only get to observe a one-dimensional time series derived by
taking a linear projection of the system states, i.e., we get to observe the
time series $s_{n}=\alpha^{T}x_{n}$ where $\alpha\in\mathbb{R}^{N}$ is an
i.i.d. Rademacher sequence. The goal is to study how well the delay coordinate
vectors formed using the time series $s_{n}$ provide a stable embedding of the
trajectory vectors $\widetilde{g}\left(x_{n}\right)$ of the dynamical system.

The discussion following Theorem II.1 suggests that the stable embedding
conditioning is dependent on how well the infimum soft-rank $r(\mathcal{A})$
scales with the number of measurements $M$. For $x\neq y\in\mathcal{A}$,
$G_{x}-G_{y}$ is a circulant matrix since $\Phi$ is a shift matrix. Moreover,
the first row of of $G_{x}-G_{y}$ is a vector of zeros except for a ‘$1$’ and
a ‘$-1$’. With these observations, the following lemma shows that the infimum
soft-rank of this matrix can be calculated analytically.

######  Lemma III.1.

For the discrete-time dynamical system described in the text, we have
$r(\mathcal{A})\geq\frac{M}{2}$.

The proof of this lemma follows from the fact that
$(G_{x}-G_{y})(G_{x}-G_{y})^{T}$ is also a circulant matrix and the formulas
for the eigenvalues of circulant matrices are well-known. This lemma tells us
that for this system, we get a linear scaling of the infimum soft-rank with
$M$. Thus Theorem II.1 will imply that the stable embedding conditioning
$\delta$ will scale with $M^{-1/2}$.

##  IV Conclusion

This paper is a first step towards understanding whether Takens’ embedding can
be stablilized through concatenating more time series measurements. Many
important issues are not yet addressed in this work. These include the study
of nonlinear measurement functions, the translation of the results from
trajectory manifold to state-space manifold, and the extension to fractal
attractors instead of just Riemannian manifold attractors. We intend to
address these issues in a future work.

## References

  * [1] P. J. Brockwell and R. A. Davis, _Introduction to Time Series and Forecasting_. Springer, 2002.
  * [2] H. Kantz and T. Schreiber, _Nonlinear Time Series Analysis_. Cambridge University Press, 2004, vol. 7.
  * [3] T. Asefa, M. Kemblowski, U. Lall, and G. Urroz, “Support Vector Machines for Nonlinear State Space Reconstruction: Application to the Great Salt Lake Time Series,” _Water Resources Research_ , vol. 41, no. 12, 2005. [Online]. Available: http://www.agu.org/pubs/crossref/2005/2004WR003785.shtml
  * [4] F. Takens, “Detecting Strange Attractors in Turbulence,” in _Dynamical Systems and Turbulence, Warwick 1980_ , ser. Lecture Notes in Mathematics, D. Rand and L.-S. Young, Eds. Springer Berlin / Heidelberg, 1981, vol. 898, pp. 366–381.
  * [5] T. D. Sauer, J. A. Yorke, and M. Casdagli, “Embedology,” _J. Statistical Physics_ , vol. 65, no. 3-4, pp. 579–616, Nov. 1991.
  * [6] S. H. Strogatz, _Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering_. Westview Press, 1994.
  * [7] E. J. Candès, “Compressive Sampling,” in _Proc. Int. Congress of Mathematicians_ , vol. 3, 2006, pp. 1433–1452.
  * [8] R. G. Baraniuk and M. B. Wakin, “Random Projections of Smooth Manifolds,” _Foundations of Computational Mathematics_ , vol. 9, no. 1, pp. 51–77, 2009\.
  * [9] K. L. Clarkson, “Tighter Bounds for Random Projections of Manifolds,” in _Proc. 24th Annu. Symp. Computational Geometry_. ACM, 2008, pp. 39–48.
  * [10] H. L. Yap, M. B. Wakin, and C. J. Rozell, “Stable Manifold Embeddings with Structured Random Matrices,” _IEEE J. Selected Topics in Signal Processing_ , vol. 7, no. 4, pp. 720–730, Sep. 2013.
  * [11] A. Eftekhari and M. B. Wakin, “New Analysis of Manifold Embeddings and Signal Recovery from Compressive Measurements,” _arXiv:1306.4748v4_ , 2014.
  * [12] P. Shah and V. Chandrasekaran, “Iterative Projections for Signal Identification on Manifolds: Global Recovery Guarantees,” in _Proc. 49th Annu. Allerton Conf._ Allerton, CA: IEEE, Sep. 2011, pp. 760–767.
  * [13] M. A. Davenport, P. T. Boufounos, M. B. Wakin, and R. G. Baraniuk, “Signal Processing with Compressive Measurements,” _IEEE J. Selected Topics in Signal Processing_ , vol. 4, no. 2, pp. 445–460, Apr. 2010.
  * [14] A. M. Fraser and H. L. Swinney, “Independent Coordinates for Strange Attractors from Mutual Information,” _Physical Review A_ , vol. 33, no. 2, pp. 1134–1140, 1986.
  * [15] M. Casdagli, S. Eubank, J. D. Farmer, and J. F. Gibson, “State Space Reconstruction in the Presence of Noise,” _Physica D: Nonlinear Phenomena_ , vol. 51, no. 1-3, pp. 52–98, 1991.
  * [16] P. Niyogi, S. Smale, and S. Weinberger, “Finding the Homology of Submanifolds with High Confidence from Random Samples,” _Technical Report No. TR-2004-08, University of Chicago_ , pp. 1–23, Mar. 2006.
  * [17] H. Federer, “Curvature measures,” _Transactions of the American Mathematical Society_ , pp. 418–491, 1959.


# What powers the radio emission in TDE AT2019dsg: a long-lived jet or the
disruption itself?

Tatsuya Matsumoto,1,2,3 Tsvi Piran,1 and Julian H. Krolik4
1Racah Institute of Physics, Hebrew University, Jerusalem, 91904, Israel
2Research Center for the Early Universe, Graduate School of Science,
University of Tokyo, Tokyo 113-0033, Japan
3Department of Physics, Graduate School of Science, University of Tokyo, Tokyo
113-0033, Japan
4Physics and Astronomy Department, Johns Hopkins University, Baltimore, MD
21218, USA
E-mail<EMAIL_ADDRESS>Research Fellow

###### Abstract

The tidal disruption event, AT2019dsg, was observed across a broad range of
electromagnetic-wavelengths from radio to X-rays, and it was possibly
accompanied by a high-energy neutrino. We study the nature of the radio-
emitting outflow by analyzing the synchrotron self-absorbed spectra in terms
of the equipartition model. We find that the time evolution of the outflow
radius can be interpreted as either free expansion or deceleration. If the
former, the outflow was launched $\simeq 40$ days before the optical peak; if
the latter, the launch was $\simeq 10$ days after the optical peak. In
addition, the energy in the radio-emitting region increases over time. This
second conclusion is most naturally interpreted by a scenario resembling the
earliest stage of a supernova remnant: as more and more material is swept up,
it is heated by the forward shock at the expense of the outflow’s kinetic
energy. Energy injection from an accreting BH cannot be completely excluded,
but the injection rate is very different from the fallback luminosity,
requiring further physical explanation. If the neutrino association is real,
the scale of energy injection needed is much greater than for the radio
emission, suggesting that the detected neutrino did not arise from the radio
emitting region.

###### keywords:

transients: tidal disruption events

††pubyear: 2021††pagerange: What powers the radio emission in TDE AT2019dsg: a
long-lived jet or the disruption itself?–References

##  1 Introduction

A supermassive black hole (BH) lurking in a center of a galaxy can destroy a
star approaching very close due to its strong gravity (Hills, 1975; Rees,
1988). Tidal disruption events (TDEs) have been observed in many wavelengths:
X-rays (Saxton et al., 2020); UV/optical (Roth et al., 2020; van Velzen et
al., 2020); and radio (Alexander et al., 2020). Understanding these
observations will help us to take a census of the population of supermassive
BHs, revealing their mass and spin distributions. It will also help us
investigate relativistic effects around BHs by shedding light on the jet-
launching mechanism, as well as stellar dynamics in galactic nuclear clusters
(Stone et al., 2020).

However, the emission mechanisms of TDEs are still under debate. This is
mainly because the evolution of bound debris, which finally accretes to the
BH, is not well understood. For instance, although the observed luminosity in
many TDEs is dominated by the UV/optical band, initial predictions estimated
that it should be very weak (Rees, 1988). Two models compete to explain it,
one applicable if debris efficiently circularizes around the BH (Metzger &
Stone, 2016), the other if it does not (Piran et al., 2015). To resolve this
issue, numerical simulations and statistical studies using increasingly large
samples (e.g. van Velzen et al., 2021) are in progress. In addition, a close
examination of individual well-observed events is still useful for
understanding TDEs.

AT2019dsg is one of the most well-observed TDEs. This event was discovered as
a bright UV/optical TDE on 9 April 2019 by the Zwicky Transient Facility (van
Velzen et al., 2021). Follow-up observation revealed rapidly declining X-rays
$\simeq 40$ days after the discovery (Cannizzaro et al., 2021; Stein et al.,
2021) and radio emission lasting more than 500 days (Stein et al., 2021;
Cendes et al., 2021; Mohan et al., 2021). In addition, IceCube reported
detection of a high-energy neutrino in the same direction as the TDE (Stein et
al., 2021), although the probability that this is actually associated with the
TDE is not high. No TDE detected since ASASSN-14li (Miller et al., 2015;
Alexander et al., 2016; Cenko et al., 2016; Holoien et al., 2016; Jiang et
al., 2016; van Velzen et al., 2016; Brown et al., 2017) has such a rich set of
multiwavelength data.

In this paper we focus on the radio, which is useful for studying the dynamics
of outflow. The radio spectra have a peak, which is caused by the synchrotron
self-absorption (SSA); applying the equipartition method to the time-
dependence of the frequency of the peak and its luminosity, we can estimate
both outflow properties and the density of the circumnuclear medium (CNM)
surrounding the BH (Chevalier, 1998; Barniol Duran et al., 2013). Stein et al.
(2021) and Cendes et al. (2021) carried out such an analysis and estimated the
radius and energy in the radio-emitting site as well as the CNM density.
Interestingly, they found that the energy in the emitting region increases
with time. They attributed the increase to energy injection from the BH. We
revisit this possibility and propose an alternative interpretation.

We organize this paper as follows. In §2 we review the equipartition method
and estimate the physical quantities of the radio-emitting outflow. We analyze
the time evolution of the equipartition radius in §3 and the CNM density
profile in §4. In §5 we discuss the implications of the increasing energy of
the radio-emitting region and find that, rather than indicating energy
injection from the BH to the radio outflow, it is more naturally interpreted
as the result of shock propagation in an impulsive event. We discuss the
implications of our result for possible neutrino emission and the origin of
the outflow in §6 and conclude this paper in §7.

##  2 The Equipartition method

Table 1: Radio data and results of our equipartition analysis for AT2019dsg. The time $\Delta t$ is measured since the discovery in the observer frame (9 April 2019), which is different from that in Cendes et al. 2021 who set the origin as 10 days before the discovery. The equipartition radius, energy, Lorentz factor and total number of emitting electrons at $\nu_{\rm p}$, and CNM density are calculated by Eqs. (1) - (5). We adopt two (freely-coasting and decelerating) velocity models given by fiting the time evolution of $R_{\rm eq}$, and the corresponding density slope $k$. The spectra of the first two epochs are not of good quality and they are excluded from our analysis. $\Delta t$ | $F_{\rm p}$ | $\nu_{\rm p}$ | $R_{\rm eq}$ | $E_{\rm eq}$ | $\gamma_{\rm e}$ | $N_{\rm e}$ |  $n$ (freely-coasting) |  $n$ (decelerating)
---|---|---|---|---|---|---|---|---
$[\rm d]$ | [mJy] | [10GHz] | [$10^{16}$cm] | [$10^{48}$erg] |  | [$10^{50}$] | [$\rm cm^{-3}$] | [$\rm cm^{-3}$]
(42) | ($0.47\pm 0.06$) | ($1.29\pm 0.30$) | ($1.04\pm 0.25$) | ($1.19\pm 0.50$) | ($59.05\pm 4.71$) | ($0.83\pm 0.26$) | ($22.29\pm 20.02$) | ($12.68\pm 11.39$)
(45) | ($0.56\pm 0.06$) | ($1.62\pm 0.49$) | ($0.90\pm 0.27$) | ($1.17\pm 0.53$) | ($59.62\pm 4.76$) | ($0.80\pm 0.29$) | ($34.01\pm 38.37$) | ($22.04\pm 24.87$)
50 | $0.74\pm 0.02$ | $1.62\pm 0.15$ | $1.02\pm 0.10$ | $1.63\pm 0.54$ | $60.53\pm 4.84$ | $1.09\pm 0.20$ | $32.01\pm 14.21$ | $24.74\pm 10.98$
70 | $0.78\pm 0.04$ | $1.07\pm 0.12$ | $1.59\pm 0.20$ | $2.63\pm 0.90$ | $60.70\pm 4.86$ | $1.75\pm 0.36$ | $13.82\pm 7.09$ | $16.81\pm 8.62$
72 | $0.74\pm 0.02$ | $1.05\pm 0.07$ | $1.58\pm 0.13$ | $2.53\pm 0.82$ | $60.53\pm 4.84$ | $1.69\pm 0.30$ | $13.35\pm 5.15$ | $16.78\pm 6.47$
120 | $1.30\pm 0.06$ | $0.89\pm 0.08$ | $2.43\pm 0.25$ | $5.85\pm 1.96$ | $62.41\pm 5.04$ | $3.71\pm 0.71$ | $8.55\pm 3.84$ | $18.31\pm 8.22$
151 | $1.11\pm 0.03$ | $0.79\pm 0.07$ | $2.53\pm 0.26$ | $5.43\pm 1.80$ | $61.88\pm 4.98$ | $3.49\pm 0.66$ | $7.03\pm 3.13$ | $18.59\pm 8.28$
178 | $1.38\pm 0.04$ | $0.50\pm 0.05$ | $4.44\pm 0.45$ | $11.17\pm 3.71$ | $62.61\pm 5.05$ | $7.05\pm 1.33$ | $2.67\pm 1.19$ | $8.16\pm 3.64$
290 | $0.89\pm 0.04$ | $0.38\pm 0.04$ | $4.76\pm 0.59$ | $8.70\pm 2.95$ | $61.14\pm 4.91$ | $5.72\pm 1.16$ | $1.69\pm 0.87$ | $7.78\pm 3.98$
551 | $0.38\pm 0.04$ | $0.18\pm 0.05$ | $6.80\pm 1.78$ | $6.69\pm 2.80$ | $58.38\pm 4.64$ | $4.76\pm 1.53$ | $0.44\pm 0.43$ | $3.40\pm 3.30$

We briefly review the equipartition method, which is applicable to a radio-
emitting outflow with a spectral peak caused by the SSA (Pacholczyk, 1970;
Scott & Readhead, 1977; Chevalier, 1998; Barniol Duran et al., 2013). Consider
an outflow producing synchrotron emission with an SSA spectral-peak flux
density $F_{\rm p}$ at a frequency $\nu_{\rm p}$. The equipartition method
gives us the radius and energy of the radio-emitting region under the
assumption that the post-shock thermal energy is distributed equally to
relativistic electrons and magnetic field.

We follow the expressions for the radius and energy of Barniol Duran et al.
(2013), but with corrections.111We take into account an additional dependence
on radius $R$ in the relativistic electrons’ energy, $\propto R^{2(1-p)}$,
which is ignored by Barniol Duran et al. (2013). This adds minor correction
factors in their original expressions. We assume that an arbitrary fraction of
the dissipated kinetic energy is transferred to the electrons
$\varepsilon_{\rm e}$ and magnetic field $\varepsilon_{\rm B}$ (see their
equations 27 and 28 in the Newtonian limit):

$\displaystyle R_{\rm eq}$ $\displaystyle\simeq 1\times 10^{17}{\,\rm
cm\,}\biggl{[}21.8\times(525)^{p-1}\gamma_{\rm
m}^{2-p}\biggl{(}\frac{p+1}{3}\biggl{)}\biggl{]}^{\frac{1}{2p+13}}$
$\displaystyle\Big{[}F_{\rm p,mJy}^{\frac{p+6}{2p+13}}d_{\rm
L,28}^{\frac{2(p+6)}{2p+13}}\nu_{\rm
p,10}^{-1}(1+z)^{-\frac{3p+19}{2p+13}}\Big{]}$
$\displaystyle\,\,\,\,\,\,f_{\rm A}^{-\frac{p+5}{2p+13}}f_{\rm
V}^{-\frac{1}{2p+13}}\Big{\\{}(4\xi\epsilon)^{\frac{1}{2p+13}}\Big{\\}}\ ,$
(1) $\displaystyle E_{\rm eq}$ $\displaystyle\simeq 1.3\times 10^{48}{\,\rm
erg\,}\biggl{(}\frac{2p+13}{17}\biggl{)}$
$\displaystyle\biggl{[}21.8\biggl{(}\frac{p+1}{3}\biggl{)}\biggl{]}^{-\frac{2(p+1)}{2p+13}}\big{[}(525)^{p-1}\gamma_{\rm
m}^{2-p}\big{]}^{\frac{11}{2p+13}}$ $\displaystyle\,\,\,\Big{[}F_{\rm
p,mJy}^{\frac{3p+14}{2p+13}}d_{\rm L,28}^{\frac{2(3p+14)}{2p+13}}\nu_{\rm
p,10}^{-1}(1+z)^{-\frac{5p+27}{2p+13}}\Big{]}f_{\rm
A}^{-\frac{3(p+1)}{2p+13}}f_{\rm V}^{\frac{2(p+1)}{2p+13}}$
$\displaystyle\,\,\,\,\,\,\biggl{\\{}(4\xi)^{\frac{11}{2p+13}}\biggl{[}\frac{11}{2p+13}\epsilon^{-\frac{2(p+1)}{2p+13}}+\frac{2(p+1)}{2p+13}\epsilon^{\frac{11}{2p+13}}\biggl{]}\biggl{\\}}\
.$ (2)

Here the observables are the power-law index of the electron distribution $p$
(given by the spectral slope) and the luminosity distance corresponding to the
redshift $z$, $d_{\rm L}$, as well as $F_{\rm p,mJy}(=F_{\rm p}/\rm mJy)$ and
$\nu_{\rm p}$. We use the notation $Q_{x}=Q/10^{x}$ in cgs units unless
otherwise specified. $\gamma_{\rm m}$ is the minimal Lorentz factor of
relativistic electrons. In the standard prescription to estimate $\gamma_{\rm
m}$ (e.g. Sari et al., 1998), it becomes less than unity for a non-
relativistic outflow with a velocity $v\lesssim v_{\rm DN}\simeq
0.1\,c\,[(p-2)\varepsilon_{\rm e,-1}/(p-1)]^{-1/2}$ (the so-called deep-
Newtonian phase, Huang & Cheng, 2003; Sironi & Giannios, 2013).222When the
outflow expands at larger velocity than $v\gtrsim v_{\rm DN}$, $\gamma_{\rm
m}$ depends on the outflow velocity and we have to determine $R_{\rm eq}$ and
$v$ at the same time. In this phase, we set $\gamma_{\rm m}=2$ because only
relativistic electrons produce synchrotron emission. The area and volume-
filling fractions of the outflow are defined so that the emission area and
volume are given by $f_{\rm A}4\pi R^{2}$ and $f_{\rm V}\pi R^{3}$,
respectively. We also define $\epsilon\equiv\big{[}(\varepsilon_{\rm
B}/\varepsilon_{\rm e})/(2(p+1)/11)\big{]}$ and $\xi\equiv 1+\varepsilon_{\rm
e}^{-1}$ (see below for the meanings of them).

The term “equipartition energy” is often used just for the sum of equally-
distributed energy of the relativistic electrons and the accompanying magnetic
field. However, the actual energy could be larger, and this possible increase
is expressed by the terms in curly brackets in Eqs. (1) and (2) (besides a
factor of 4 arising from a correction to the isotropic number of radiating
electrons in the non-relativistic case). The factor $\epsilon$ describes a
deviation of the ratio of magnetic to electron energy from $2(p+1)/11$, the
value for the minimal energy case. The factor $\xi$ represents additional
energy stored in the (hot) protons existing in the shocked region. With these
corrections, the equipartition energy in Eq. (2) is the sum of the energy of
(hot) protons, relativistic electrons, and magnetic field in the emitting
site. Note that the outflow’s total energy (including kinetic energy, for
example) could be much larger, as we discuss later.

Importantly the factors $\xi$ and $\epsilon$ as well as the geometrical
parameters $f_{\rm A}$ and $f_{\rm V}$ hardly influence $R_{\rm eq}$, so this
argument enables us to determine the emitting radius with only small
uncertainty. The most significant dependence of $R_{\rm eq}$ comes from the
covering fraction of the emitting region, $f_{\rm A}$, and it is roughly given
by $R_{\rm eq}\propto f_{\rm A}^{-0.5}$.

The Lorentz factor and number of radiating electrons at $\nu_{\rm p}$ are
given by (equations 14 and 15 in Barniol Duran et al., 2013)

$\displaystyle\gamma_{\rm e}$ $\displaystyle\simeq 525\,F_{\rm p,mJy}d_{\rm
L,28}^{2}\nu_{\rm p,10}^{-2}(1+z)^{-3}f_{\rm A}^{-1}R_{\rm eq,17}^{-2}\ ,$ (3)
$\displaystyle N_{\rm e}$ $\displaystyle\simeq 4.0\times 10^{54}\,F_{\rm
p,mJy}^{3}d_{\rm L,28}^{6}\nu_{\rm p,10}^{-5}(1+z)^{-8}f_{\rm A}^{-2}R_{\rm
eq,17}^{-4}\ ,$ (4)

respectively. The (pre-shock) CNM density is estimated by dividing the total
number of electrons by the volume swept up by the outflow:

$\displaystyle n$ $\displaystyle=\frac{N_{\rm e}\big{(}\frac{\gamma_{\rm
m}}{\gamma_{\rm e}}\big{)}^{1-p}{\rm max}\big{[}\big{(}\frac{v_{\rm
DN}}{v}\big{)}^{2},\,1\big{]}}{\frac{\Omega}{3-k}R_{\rm eq}^{3}}$ (5)
$\displaystyle\simeq 4100{\,\rm cm^{-3}\,}\big{[}525\gamma_{\rm
m}^{-1}\big{]}^{p-1}F_{\rm p,mJy}^{p+2}d_{\rm L,28}^{2(p+2)}\nu_{\rm
p,10}^{-2p-3}$ $\displaystyle\,\,\,\,\,\,\times(1+z)^{-3p-5}f_{\rm
A}^{-p-1}R_{\rm eq,17}^{-2p-5}\Omega^{-1}(3-k){\rm
max}\biggl{[}\Big{(}\frac{v_{\rm DN}}{v}\Big{)}^{2},\,1\biggl{]}\ .$

Here three correction factors arise. First, most relativistic electrons have
Lorentz factor $\sim\gamma_{\rm m}$, and their number is $(\gamma_{\rm
m}/\gamma_{\rm e})^{1-p}$ times larger than $N_{\rm e}$. Second, when the
outflow is in the deep-Newtonian phase, only a fraction $(v/v_{\rm DN})^{2}$
of electrons are accelerated to relativistic energy (Sironi & Giannios, 2013;
Matsumoto & Piran, 2021a). Third, the density is assumed to follow $n\propto
R^{-k}$ ($k<3$), and the total electron number in the volume swept-up by the
outflow with a constant solid angle $\Omega$ is given by $\int_{0}^{R}\Omega
r^{2}n(r)dr=\Omega n(R)R^{3}/(3-k)$. The index $k$ is known only after the
density profile is calculated without this factor. When the density profile is
not described by a power-law function, we can approximate the integral as
$\simeq\Omega n(R)R^{3}$, which is justified as long as $d\log n/d\log r<3$.

We apply the equipartition method to AT2019dsg using the peak fluxes and
frequencies in Table 1 taken from Stein et al. (2021) and Cendes et al.
(2021). The other observables are $z=0.051$, the corresponding luminosity
distance $d_{\rm L}\simeq 230\,\rm Mpc$, and the power-law index $p=2.7\pm
0.2$. We adopt the same parameters as those of Cendes et al. (2021), who
assumed that the outflow is spherical ($\Omega=4\pi$), and the emitting region
is a shell with the width of $dR=0.1\,R_{\rm eq}$. These assumptions give
$f_{\rm A}=1$ and $f_{\rm V}=0.36$. The other parameters are $\varepsilon_{\rm
e}=0.1$, $\varepsilon_{\rm B}=0.02$, and $\gamma_{\rm m}=2$ (assuming
$v<v_{\rm DN}$, which is justified later). The results are shown in Table 1,
which are basically consistent with those of Cendes et al. (2021) except for
the CNM density; in that estimate, they neglected the correction factor
$(v/v_{\rm DN})^{2}$. For the outflow velocity, we use two (freely-coasting
and decelerating) models motivated by fitting the time evolution of $R_{\rm
eq}$ (see the next section). We also adopt the density slope $k$ obtained by
fitting the calculated profile.

##  3 Launch time of the radio-outflow

Figure 1: (Bottom) Equipartition radius at each observation epoch. The red and
blue curves denote the best-fit curves of freely-coasting and decelerating
models, respectively, with key parameters of $v_{\rm sh}\simeq 0.064\,c$ and
$t_{l}\simeq-18\,\rm d$ (freely-coasting) and $\alpha=0.63$ and $t_{l}=26\,\rm
d$ (decelerating). The quality of spectrum at 551 days is not so good and we
use a gray dot to highlight this. The orange shaded region shows the peak time
of optical emission. (Top) Outflow velocity derived by time derivative of the
best-fit models of $R_{\rm eq}$.

Fitting the time evolution of the equipartition radius, we derive the velocity
and launching time of the radio-emitting outflow (hereafter radio-outflow).
Fig. 1 depicts the equipartition radius at each observation time. Since the
quality of the data at 42 and 45 days after the discovery is too poor to
determine the spectral peak, we do not include these observations in the
following calculation. The frequency of the peak is also somewhat difficult to
determine in the $\Delta t=551$ d data, but we include it in analysis with a
sizable error bar.

We consider two evolution models. The first is a freely-coasting outflow:
$R_{\rm eq}=v_{\rm sh}(\Delta t-t_{l})$, where $\Delta t$ and $t_{l}$ are the
observation time and outflow-launching time, respectively, and both of them
are measured since the discovery. The best-fit parameters are $v_{\rm
sh}=0.064_{-0.005}^{+0.005}\,c\simeq 20000_{-1500}^{+1500}\,\rm km\,s^{-1}$
and $t_{l}=-18_{-11}^{+9}\,\rm d$ (the error represents $1\sigma$), giving a
reduced chi-square of $\chi_{\rm r}^{2}\simeq 2.2$. The second is a
decelerating (or accelerating) outflow with $R_{\rm eq}=A(\Delta
t-t_{l})^{\alpha}$. The best-fit parameters are $t_{l}=26_{-18}^{+11}\,\rm d$
and $\alpha=0.63_{-0.12}^{+0.16}$ with a reduced $\chi_{\rm r}^{2}\simeq
1.8$.; hence, deceleration is favored, as found by Cendes et al. (2021).
However, the distinction in $\chi_{\rm r}^{2}$ depends entirely on the $\Delta
t=551$ d observation, which has a large uncertainty.

It should be noted that for both models the best-fit launching time $t_{l}$ is
independent of the equipartition parameters $\epsilon$, $\xi$, $f_{\rm A}$,
and $f_{\rm V}$ (Krolik et al., 2016) while the other parameters $v_{\rm sh}$
and $A$ depend on them. These best-fit curves are shown in Fig. 1.

The top panel of Fig. 1 depicts the velocity of each model as a function of
time, deriving it from $v_{\rm sh}\equiv dR_{\rm eq}/d(\Delta t)$. Strictly
speaking, this quantity represents not the outflow velocity but the shock
velocity. The actual outflow velocity is a little smaller than the shock
velocity; their ratio is given by the shock jump condition (e.g. Landau &
Lifshitz, 1987): $\zeta\equiv v_{\rm sh}/v=(\hat{\gamma}+1)/2$, where
$\hat{\gamma}$ is the adiabatic index of the shocked material. For
$\hat{\gamma}=5/3$, this factor is order unity, $\zeta=4/3$, but this
correction is important for inferring the origin of the outflow.

The outflow-launching time $t_{l}$ can also constrain the origin of the radio
outflow. We define it using the peak time of the optical emission as a
reference time and taking the time difference between the outflow’s launch and
the optical peak. The optical luminosity peaked at around 27 April 2019, or
$18_{-10}^{+8}$ days after discovery (van Velzen et al., 2021). For the
freely-coasting model, the optical peak was $36_{-10}^{+8}$ days after the
outflow-launch. On the other hand, for the decelerating model the optical peak
precedes the outflow-launch by $8^{+10}_{-8}$ days (Here we used the central
values of $t_{l}$. If we add errors of $t_{l}$ assuming a Gaussian
distribution with $\sigma\simeq 10$ and $15\,\rm d$ the launching time has an
error of $\simeq 13$ and $17\,\rm d$ for each model). We will discuss the
implication of these time differences for the origin of the radio outflow in
§6.

##  4 The CNM density profile

Radio observations of TDEs provide possibly the only way to infer the CNM
density distribution around distant galaxies (e.g. Alexander et al., 2016;
Krolik et al., 2016). Fig. 2 depicts the profiles reconstructed by our
equipartition analysis. The two thick colored curves represent the predictions
made by our two outflow models (note that the estimate of the CNM density
depends on the outflow velocity, see Eq. 5). The velocities are derived from
the best-fit models of the equipartition radius, $v=\zeta^{-1}dR_{\rm
eq}/d(\Delta t)$. Fitting the density profiles with a power-law function, we
find the best fits for $n\propto R^{-1.99}$ and $R^{-0.83}$ for the freely-
coasting and decelerating models, respectively. Our density values are
slightly larger than those of Cendes et al. (2021), mainly because we consider
the deep-Newtonian correction, $(v/v_{\rm DN})^{2}$.

The CNM profiles of our Galactic center and the host galaxies of two other
radio TDEs, ASASSN-14li (Alexander et al., 2016; Krolik et al., 2016) and CNSS
J0019+00 (Anderson et al., 2020) are shown in Fig. 2. The latter two are
derived by the same procedure that we used for AT2019dsg. Their equipartition
radii evolve at a constant velocity of $v_{\rm sh}\simeq 0.03\,c$. The CNM
densities in all the TDE host galaxies are larger than in Sgr A*. The CNM
density profile of AT2019dsg is shallower than in the other TDEs; this
contrast is related to the different time evolution of its equipartition
energy.

Figure 2: CNM density profiles reconstructed from the equipartition analysis.
Red and blue line shows the best-fit power-law functions for the freely-
coasting ($n\propto R^{-1.99}$) and decelerating ($n\propto R^{-0.83}$) models
for AT2019dsg, respectively. Profiles of Milky Way ($n\propto R^{-1}$,
Baganoff et al., 2003; Gillessen et al., 2019) and other TDEs ($n\propto
R^{-2.5}$) are also presented. Black dash-dotted lines show the positions
where the enclosed mass of $M(R)\simeq 4\pi
m_{\rm{p}}nR^{3}=10^{-4}-10^{-2}\,\mathrm{M}_{\odot}$.

##  5 Energy injection from the BH?

Figure 3: (Top) Equipartition energy at each observation epoch. The red and
blue curves represent the best-fit power-law functions with the same outflow-
launching time as the freely-coasting ($t_{l}=-18\,\rm d$) and decelerating
($26\,\rm d$) outflow models, respectively. (Bottom) Ratio of equipartition
energy to the kinetic energy.

We now address the question raised in the Introduction: Is energy continually
injected into the outflow, as inferred by Stein et al. (2021) and Cendes et
al. (2021)? To answer this question, we fit the time evolution of $E_{\rm eq}$
by a power-law function $E_{\rm eq}\propto(\Delta t-t_{l})^{\beta}$ using the
same $t_{l}$ obtained from the fits to $R_{\rm eq}$. The best-fit power-law
indices are $\beta=0.71_{-0.15}^{+0.14}$ and $0.52_{-0.11}^{+0.10}$ for the
freely-coasting and deceleration models, respectively (see top panel in Fig.
3). Stein et al. (2021) and Cendes et al. (2021) found that the equipartition
energy increases as $E_{\rm eq}\propto\tilde{t}$ and $\tilde{t}^{1.5}$,
respectively, where $\tilde{t}$ is the time measured since their own origins.
They found faster energy increase because they used data points only up to 178
days, but we agree that the energy within the radio-producing region increases
as time goes on.

Stein et al. (2021) and Cendes et al. (2021) suggested that the increase is
caused by a continuous energy injection to the shocked region from the
accreting BH. In the following paragraphs, we show that that the observations
do not necessarily imply this. In fact, although we cannot rule it out, we
demonstrate that this interpretation is unlikely. Consider, first, the freely-
coasting (constant velocity) model. Like a supernova remnant, the first stage
of the interaction between the outflow and the surrounding matter involves a
combination of two shocks: a forward shock propagating into the CNM and a
reverse shock propagating into the outflow. Their post-shock regions are
separated by a contact discontinuity. As long as the mass of the outflow is
significantly larger than the shocked CNM mass, the system effectively does
not slow down. As more CNM is swept up and heated by the forward shock, more
energy is transferred to the shocked region at the expense of the very large
kinetic energy of the outflow.333Energy loss by electron synchrotron cooling
is negligible because the cooling timescale is much longer than the dynamical
time, $t_{\rm cool}\simeq 4500{\,\rm d\,}(B/{\rm G})^{-2}(\gamma_{\rm
m}/2)^{-1}$, where the magnetic field is estimated by the equiaprtition
analysis, $B\simeq 0.1-1\,\rm G$ (see also Cendes et al., 2021). This energy
transfer appears as an increase of the energy of the relativistic electrons
and magnetic field, and hence as an increase of the equipartition energy. This
phase of outflow-CNM interaction most naturally explains the observation of
AT2019dsg.

More precisely, the velocity of the shocked region is determined by energy
conservation:

$\displaystyle E_{\rm kin}(>v)=\big{[}M_{\rm ej}(>v)+M(R_{\rm
eq})\big{]}v^{2}/2\ ,$ (6)

where we suppose that the outflow is ejected with a spread of velocities, with
kinetic energy and mass distributions $E_{\rm kin}(>v)$ and $M_{\rm ej}(>v)$.
The CNM mass swept-up by an isotropic outflow is approximated by

$\displaystyle M(R)\simeq 4\pi m_{\rm{p}}nR^{3}\ .$ (7)

In the freely-coasting phase, as long as the outflow mass dominates the total
mass, $M_{\rm ej}(>v)\gg M(R_{\rm eq})\simeq 3\times
10^{-3}\,\mathrm{M}_{\odot}$, the velocity is almost constant. Therefore, the
velocity structure should be very steep in this case (for instance, $E_{\rm
kin}(>v)\propto v^{-l}$ with $l\gg 1$) and the outflow may be characterized by
a single velocity. On the other hand, the equipartition energy is given by

$\displaystyle E_{\rm eq}\sim m_{\rm{p}}nf_{\rm V}\pi R_{\rm
eq}^{3}v^{2}/2\propto M(R_{\rm eq})v^{2}\ .$ (8)

Therefore, without any external energy injection, the equipartition energy
naturally increases as more and more CNM is swept up and accumulates in the
shocked region. To realize an energy increase with
$\beta=0.71_{-0.15}^{+0.14}$ in the freely-coasting model, the density profile
should have $k=2.29_{-0.14}^{+0.15}$,444For a given evolution of $R_{\rm
eq}\propto t^{\alpha}$ and $E_{\rm eq}\propto t^{\beta}$, we find these
exponents are related to the slope of density profile by
$k=(5\alpha-\beta-2)/\alpha$. that is consistent with our result,
$k=1.99_{-0.26}^{+0.28}$. It should be noted that in Eq. (8), we assumed that
the total energy in the emission region is proportional to the dissipated
kinetic energy of the outflow. The bottom panel of Fig. 3 shows that the ratio
of $E_{\rm eq}/[M(R_{\rm eq})v^{2}/2]$ supports this assumption.

Although there is no motivation to consider external energy injection in the
freely-coasting model, we can still consider such a scenario. In particular,
suppose the BH launches a relativistic jet that increases the kinetic energy
of the outflow as $E_{\rm kin}\propto t^{m}$ (here we define $t\equiv\Delta
t-t_{l}$). Depending on whether the outflow mass or swept-up CNM mass
dominates the total mass, we obtain the time evolution of $R_{\rm eq}$:

$\displaystyle R_{\rm eq}\propto\begin{cases}t^{\frac{m+2}{2}}&:M_{\rm
ej}>M(R_{\rm eq})\ ,\\\ t^{\frac{m+2}{5-k}}&:M_{\rm ej}<M(R_{\rm eq})\
,\end{cases}$ (9)

where we consider a power-law CNM distribution $n\propto R^{-k}$ and
approximate $R_{\rm eq}\sim v_{\rm sh}t$. In our fit with $\alpha=1$, the
former case corresponds to a “genuine” freely-coasting outflow,
$m=2(\alpha-1)=0$, which we have proposed above. The outflow mass must be
larger than $\gtrsim 3\times 10^{-3}\,\mathrm{M}_{\odot}$. The latter case
represents an injection scenario, where the outflow would have decelerated
without the energy supply. We find here that
$m=\alpha(5-k)-2=1.01_{-0.28}^{+0.26}$, similar to the value of
$\beta=0.71_{-0.15}^{+0.14}$. In this case, the outflow mass is constrained to
be $\lesssim 3\times 10^{-4}\,\mathrm{M}_{\odot}$. However, as we have pointed
out, there is no a priori reason why the engine should inject energy at this
rate in order to keep the outflow velocity constant. These results are
summarized in Table 2.

In the decelerating model, the swept-up CNM mass must be larger than the
outflow mass, $M(R)\gtrsim M_{\rm ej}$. Eq. (8) implies $E_{\rm kin}\propto
E_{\rm eq}$ and hence the energy in the shocked region increases. Again there
are two possibilities for this increase. One is that it can arise from the
kinetic energy of the outflow. Slower materials catch up with and inject
energy to the shocked region (this scenario is similar to the scenario for the
freely-coasting model). The CNM mass is comparable to the outflow mass $M_{\rm
ej}(>v)$ in this case. Assuming the kinetic energy distribution $E_{\rm
kin}(>v)\propto v^{-l}$ we obtain the evolution of the radius as $R\propto
t^{\frac{l+2}{5-k+l}}$. The slope of the energy distribution is given by
$l=(2+k\alpha-5\alpha)/(\alpha-1)\simeq 1.7$ and the energy distribution is
inferred to be $E_{\rm kin}(>v)\simeq 10^{48}{\,\rm erg\,}(v/0.08\,c)^{-1.7}$
or equivalently $M_{\rm ej}(>v)\simeq 3\times
10^{-4}{\,\mathrm{M}_{\odot}}(v/0.08\,c)^{-3.7}$. At the last epoch of
observation, about $\simeq 9\times 10^{48}\,\rm erg$ has been injected in
total. Again this scenario does not require any energy injection from the BH
engine.

Another is continuous energy injection from the BH. In this case, $M_{\rm
ej}\ll M(R)$ and the time evolution of the radius is given by Eq. (9). We find
that the injected energy evolves as $E_{\rm kin}\simeq 10^{48}{\,\rm
erg\,}(t/24\,{\rm d})^{0.6}$ and the injection luminosity is given by $\simeq
2\times 10^{41}{\,\rm erg\,s^{-1}\,}(t/24\,{\rm d})^{-0.4}$, where the
normalization corresponds to the first observation $\Delta t=50\,\rm d$. This
luminosity is much smaller and decays rather slower than the fallback
luminosity to the BH, $\dot{M}_{\rm fb}c^{2}\sim 10^{47}{\,\rm
erg\,s^{-1}\,}(t/24{\,\rm d})^{-5/3}R_{*,0}M_{*,0}^{1/3}M_{\rm BH,6}^{-2/3}$
(Rees, 1988; Phinney, 1989), where $R_{*}$, $M_{*}$, and $M_{\rm BH}$ are the
disrupted stellar radius and mass, and the BH mass, respectively, normalized
to solar values. We cannot exclude this possibility, but we do not see a
simple physical explanation for this injection. Most importantly, it is much
smaller than the injection rate inferred from the neutrino association (Stein
et al., 2021; Winter & Lunardini, 2021, see below).

Table 2: Possible types of energy injection for freely-coasting and decelerating radio-outflow models of AT2019dsg. The indices $\alpha$, $\beta$, and $k$ are obtained by fitting radius ($R_{\rm eq}\propto t^{\alpha}$), energy ($E_{\rm eq}\propto t^{\beta}$), and density profile ($n\propto R^{-k}$), respectively. For each injection type, the inferred injection profiles are shown. These indexes are given by $m=2(\alpha-1)$ for $M_{\rm ej}>M(R)$, $m=\alpha(5-k)-2$ for $M_{\rm ej}<M(R)$, and $l=(2+k\alpha-5\alpha)/(\alpha-1)$ (see text). Model | Index | Type of energy injection
---|---|---
| $\alpha$ | $\beta$ | $k$ |  $E_{\rm kin}\propto t^{m}$ ($M_{\rm ej}>M(R)$) |  $E_{\rm kin}\propto t^{m}$ ($M_{\rm ej}<M(R)$) | $E_{\rm kin}\propto v^{-l}$
Freely-coasting | 1 | $0.71_{-0.15}^{+0.14}$ | $1.99_{-0.26}^{+0.28}$ |  $m=0$ (Natural) |  $m\simeq 1.0$ (Fine-tuning) | —
Decelerating | $0.63_{-0.12}^{+0.16}$ | $0.52_{-0.11}^{+0.10}$ | $0.83_{-0.26}^{+0.28}$ | — |  $m\simeq 0.6$ (Possible) |  $l\simeq 1.7$ (Possible)

##  6 Discussion

AT2019dsg coincides in position with the high-energy ($\simeq 0.2\,\rm PeV$)
neutrino IceCube-101001A detected on 1 October 2019 (175 days after the
discovery, Stein et al., 2021). Stein et al. (2021) suggested that the
neutrino was produced by a proton accelerated by a radio-emitting outflow
whose energy is as large as $\sim 10^{53}\,\rm erg$. This is much larger than
the energy estimates based on the radio observations. This situation does not
change even if we consider a conical outflow like a jet. Winter & Lunardini
(2021) proposed that the neutrino was produced by an on-axis jet which
continuously injects energy to the radio-emitting site. Similarly, their jet
luminosity of $\simeq 3\times 10^{44}\,\rm erg\,s^{-1}$ and injected energy of
$\simeq 5\times 10^{51}\,\rm erg$ are three orders of magnitude larger than
those implied by our estimates of the energy in the emitting region. The large
discrepancy in energy could be resolved with very small equipartition
parameters $\varepsilon_{\rm e}$ and $\varepsilon_{\rm B}$ are $\lesssim
10^{-4}$. However, the radio outflow would have been accelerated to
relativistic velocity with these injection rates, which are inconsistent with
the inferred emission radii form our equipartition analysis. Moreover, such a
relativistic outflow was not seen by the VLBI observation (Mohan et al.,
2021). Other scenarios where the neutrino production site is not related to
the radio-emitting outflow (Hayasaki & Yamazaki, 2019; Liu et al., 2020;
Murase et al., 2020) are not constrained by our results.

Turning now to the origin of the radio outflow, there are two main channels to
launch Newtonian outflows from TDEs. One is the unbound debris ejected at the
moment of the disruption. The other is a disk wind launched when debris
fallback forms an accretion disk. Because these outflows are launched at
different times, it may be possible to discriminate between them on the basis
of our inferred launching times, $\simeq 40\,\rm d$ before and $\simeq
10{\,\rm d}$ after the optical peak, for the freely-coasting and decelerating
models, respectively. The distribution of the reduced chi-square is shown in
Figs. 5 and 5 for each model.

Figure 4: Distribution of the reduced chi-square for the freely-coasting model
of $R_{\rm eq}=v(t-t_{l})$. The yellow star denotes the best-fit parameter set
of $v=v_{\rm sh}/\zeta\simeq 0.05\,c\simeq 15000\rm\,km\,s^{-1}$ and
$t_{l}\simeq-18\,\rm d$.

Figure 5: The same as Fig. 5 but for the deceleration model of $R_{\rm
eq}\propto(t-t_{l})^{\alpha}$. The best-fit parameters are $\alpha\simeq 0.63$
and $t_{l}\simeq 26\,\rm d$.

These launch times can constrain the outflow’s origin when placed in the
context of an optical emission model. According to the apocenter-shock model
(e.g. Piran et al., 2015), the returning stream passes the pericenter having
suffered little dissipation and collides with later-arriving material near the
orbital apocenter. The optical emission arises from the apocenter shock, and
it peaks when the most bound debris returns to the apocenter after its
fallback time:555Radiation diffusion from the shocked region can lengthen the
timescale.

$\displaystyle t_{\rm fb}\simeq 41{\,\rm d\,}R_{*,0}^{3/2}M_{*,0}^{-1}M_{\rm
BH,6}^{1/2}\ .$ (10)

Thus, the unbound debris is launched $\sim t_{\rm fb}$ before the optical
peak. Ryu et al. (2020a) estimated the timescale $t_{\rm fb}=29_{-4}^{+5}\,\rm
d$ within the apocenter-shock model,666They took into account the order-unity
corrections due to stellar structure and relativistic effects in Eq. (10).
which is consistent with our result $\simeq 40\,\rm d$ (freely-coasting
model). However, the expected solid angle of the unbound debris is $\Omega\sim
0.1-1$ (Krolik et al., 2016; Yalinewich et al., 2019), for which the implied
velocity is $\gtrsim 3$ times larger than the $v\simeq 15000\,\rm km\,s^{-1}$
estimated for $\Omega=4\pi$ (see Eq. 1 and the dependence on $f_{\rm
A}\propto\Omega$). Such a velocity is too high for the unbound debris even if
we consider a fast tail in the velocity distribution (Ryu et al., 2020b;
Matsumoto & Piran, 2021a). Note that the launching time is independent of the
equipartition parameters.

The other optical emission model is the reprocessing model (e.g. Loeb & Ulmer,
1997; Strubbe & Quataert, 2009; Metzger & Stone, 2016) in which an accretion
disk forms when the fallback matter returns to the pericenter.777However, a
naive estimate of the outflow mass suggests that this model requires more mass
than is typically involved in TDEs (Matsumoto & Piran, 2021b). The accretion
rate is larger than the Eddington rate and a massive disk wind is launched,
blocking X-rays from the disk and reprocessing them to optical photons. Such a
massive disk wind can be the radio outflow for the freely-coasting model.
Neither optical emission model predicts outflows launched after the optical
peak, as found in the decelerating model.

##  7 Summary

TDE AT2019dsg presents rich data not only in the optical/UV and X-rays but
also in the radio bands, thereby giving us a good opportunity to study the
dynamics of the radio-emitting outflow. In this work, we revisit the radio
data and results of the equipartition method carried out by Stein et al.
(2021); Cendes et al. (2021). This method gives us the radius and energy
within the radio-emitting site at each observation time when the SSA spectral
peak is available. For the equipartition radius, we fit the evolution with
linear and power-law functions representing freely-coasting and decelerating
outflows, respectively. Assuming isotropic geometry, the former model gives an
outflow velocity $\simeq 1.5\times 10^{4}\,\rm km\,s^{-1}$. The latter model
requires $\simeq 2.5\times 10^{4}{\,\rm km\,s^{-1}\,}(t/24{\,\rm d})^{-0.4}$,
and the inferred CNM density distribution is shallower than that of the former
one due to the different velocity evolution.

The equipartition energy increases with time. This was interpreted as due to
an energy injection from the BH by Stein et al. (2021) and Cendes et al.
(2021). In both the freely-coasting and decelerating models, if there is
continual energy injection, it requires some sort of fine tuning. For the
freely-coasting model, the CNM density profile must be matched to the rate of
energy injection from the central engine in order to keep the velocity
constant. For the decelerating model, the injection rate of $\dot{E}_{\rm
kin}\propto t^{-0.4}$ is much different from the mass fallback rate to the BH.
Moreover the total energy is very small compared to the overall energy budget
of the system.

On the other hand, the increase in energy can be naturally interpreted as
arising from the kinetic energy of the outflow. In the freely-coasting model,
the outflow mass is larger than the external mass and hence the velocity
automatically remains essentially constant. As the outflow sweeps up more and
more matter, the energy within the radio-emitting region increases. In this
case, the external density $n\propto R^{-2}$. In the decelerating model, the
outflow mass is small so that sweeping up material slows it down; the
deceleration then results in originally-slower matter catching up, delivering
additional energy.

To summarize, we find that the basic picture in which radio emission is
produced by a shock running into external material naturally leads to
continuing energy transfer from material driving the shock into the shocked
region. Thus, in any shock-driven system in which the cooling time is longer
than the dynamical time, and in particular here and in other radio TDEs, this
is the simplest and most natural explanation for the gradual growth of the
energy in the radio emitting region.

## acknowledgments

We thank Chi-Ho Chan, Assaf Horesh, and Christopher Irwin for fruitful
discussions and helpful comments. T.M. thanks the Yukawa Institute for
Theoretical Physics at Kyoto University. Discussions during the YITP workshop
YITP-T-21-05 on “Extreme Outflows in Astrophysical Transients” were useful for
this work. This work is supported in part by JSPS Postdoctral Fellowship,
Kakenhi No. 19J00214 (T.M.) and by ERC advanced grant “TReX” (T.P.).

## data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Alexander et al. (2016) Alexander K. D., Berger E., Guillochon J., Zauderer B. A., Williams P. K. G., 2016, ApJ, 819, L25
  * Alexander et al. (2020) Alexander K. D., van Velzen S., Horesh A., Zauderer B. A., 2020, Space Sci. Rev., 216, 81
  * Anderson et al. (2020) Anderson M. M., et al., 2020, ApJ, 903, 116
  * Baganoff et al. (2003) Baganoff F. K., et al., 2003, ApJ, 591, 891
  * Barniol Duran et al. (2013) Barniol Duran R., Nakar E., Piran T., 2013, ApJ, 772, 78
  * Brown et al. (2017) Brown J. S., Holoien T. W. S., Auchettl K., Stanek K. Z., Kochanek C. S., Shappee B. J., Prieto J. L., Grupe D., 2017, MNRAS, 466, 4904
  * Cannizzaro et al. (2021) Cannizzaro G., et al., 2021, MNRAS, 504, 792
  * Cendes et al. (2021) Cendes Y., Alexander K. D., Berger E., Eftekhari T., Williams P. K. G., Chornock R., 2021, arXiv e-prints, p. arXiv:2103.06299
  * Cenko et al. (2016) Cenko S. B., et al., 2016, ApJ, 818, L32
  * Chevalier (1998) Chevalier R. A., 1998, ApJ, 499, 810
  * Gillessen et al. (2019) Gillessen S., et al., 2019, ApJ, 871, 126
  * Hayasaki & Yamazaki (2019) Hayasaki K., Yamazaki R., 2019, ApJ, 886, 114
  * Hills (1975) Hills J. G., 1975, Nature, 254, 295
  * Holoien et al. (2016) Holoien T. W. S., et al., 2016, MNRAS, 455, 2918
  * Huang & Cheng (2003) Huang Y. F., Cheng K. S., 2003, MNRAS, 341, 263
  * Jiang et al. (2016) Jiang N., Dou L., Wang T., Yang C., Lyu J., Zhou H., 2016, ApJ, 828, L14
  * Krolik et al. (2016) Krolik J., Piran T., Svirski G., Cheng R. M., 2016, ApJ, 827, 127
  * Landau & Lifshitz (1987) Landau L. D., Lifshitz E. M., 1987, Fluid Mechanics
  * Liu et al. (2020) Liu R.-Y., Xi S.-Q., Wang X.-Y., 2020, Phys. Rev. D, 102, 083028
  * Loeb & Ulmer (1997) Loeb A., Ulmer A., 1997, ApJ, 489, 573
  * Matsumoto & Piran (2021a) Matsumoto T., Piran T., 2021a, arXiv e-prints, p. arXiv:2107.06289
  * Matsumoto & Piran (2021b) Matsumoto T., Piran T., 2021b, MNRAS, 502, 3385
  * Metzger & Stone (2016) Metzger B. D., Stone N. C., 2016, MNRAS, 461, 948
  * Miller et al. (2015) Miller J. M., et al., 2015, Nature, 526, 542
  * Mohan et al. (2021) Mohan P., An T., Zhang Y., Yang J., Yang X., Wang A., 2021, arXiv e-prints, p. arXiv:2106.15799
  * Murase et al. (2020) Murase K., Kimura S. S., Zhang B. T., Oikonomou F., Petropoulou M., 2020, ApJ, 902, 108
  * Pacholczyk (1970) Pacholczyk A. G., 1970, Radio astrophysics. Nonthermal processes in galactic and extragalactic sources
  * Phinney (1989) Phinney E. S., 1989, in Morris M., ed., IAU Symposium Vol. 136, The Center of the Galaxy. p. 543
  * Piran et al. (2015) Piran T., Svirski G., Krolik J., Cheng R. M., Shiokawa H., 2015, ApJ, 806, 164
  * Rees (1988) Rees M. J., 1988, Nature, 333, 523
  * Roth et al. (2020) Roth N., Rossi E. M., Krolik J., Piran T., Mockler B., Kasen D., 2020, Space Sci. Rev., 216, 114
  * Ryu et al. (2020a) Ryu T., Krolik J., Piran T., 2020a, ApJ, 904, 73
  * Ryu et al. (2020b) Ryu T., Krolik J., Piran T., Noble S. C., 2020b, ApJ, 904, 99
  * Sari et al. (1998) Sari R., Piran T., Narayan R., 1998, ApJ, 497, L17
  * Saxton et al. (2020) Saxton R., Komossa S., Auchettl K., Jonker P. G., 2020, Space Sci. Rev., 216, 85
  * Scott & Readhead (1977) Scott M. A., Readhead A. C. S., 1977, MNRAS, 180, 539
  * Sironi & Giannios (2013) Sironi L., Giannios D., 2013, ApJ, 778, 107
  * Stein et al. (2021) Stein R., et al., 2021, Nature Astronomy, 5, 510
  * Stone et al. (2020) Stone N. C., Vasiliev E., Kesden M., Rossi E. M., Perets H. B., Amaro-Seoane P., 2020, Space Sci. Rev., 216, 35
  * Strubbe & Quataert (2009) Strubbe L. E., Quataert E., 2009, MNRAS, 400, 2070
  * Winter & Lunardini (2021) Winter W., Lunardini C., 2021, Nature Astronomy, 5, 472
  * Yalinewich et al. (2019) Yalinewich A., Steinberg E., Piran T., Krolik J. H., 2019, MNRAS, 487, 4083
  * van Velzen et al. (2016) van Velzen S., et al., 2016, Science, 351, 62
  * van Velzen et al. (2020) van Velzen S., Holoien T. W. S., Onori F., Hung T., Arcavi I., 2020, Space Sci. Rev., 216, 124
  * van Velzen et al. (2021) van Velzen S., et al., 2021, ApJ, 908, 4