#  Topological comparison of some dimension reduction methods using persistent
homology on EEG data

Eddy Kwessi111Corresponding author: Department of Mathematics, Trinity
University, 1 Trinity Place, San Antonio, TX 78212, Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we explore how to use topological tools to compare dimension
reduction methods. We first make a brief overview of some of the methods often
used dimension reduction such as Isometric Feature Mapping, Laplacian
Eigenmaps, Fast Independent Component Analysis, Kernel Ridge Regression,
t-distributed Stochastic Neighbor Embedding. We then give a brief overview of
some topological notions used in topological data analysis, such as, barcodes,
persistent homology, and Wasserstein distance. Theoretically, these methods
applied on a data set can be interpreted differently. From EEG data embedded
into a manifold of high dimension, we apply these methods and we compare them
across persistent homologies of dimension 0, 1, and 2, that is, across
connected components, tunnels and holes, shells around voids or cavities. We
find that from three dimension clouds of points, it is not clear how distinct
from each other the methods are, but Wasserstein and Bottleneck distances,
topological tests of hypothesis, and various methods show that the methods
qualitatively and significantly differ across homologies.

##  1 Introduction

In topological data analysis, one is interested in understanding high
dimensional structures from low dimensional ones and how discrete structures
can be aggregated to form a global structure. It can be a difficult task to
even think or believe that high dimensional object exist beyond three
dimensions since we can not visualize objects beyond a three-dimensional
space. However, embedding theorems, for instance Whitney (1936) and Takens
(1981) embedding theorems, clearly show that these high dimension structures
do in fact exist. From a practical point of view, to make inferences on
structures embedded in high dimensional ambient spaces, some kind of
dimensional reduction needs to occur. From a data analysis point of view,
dimension reduction amounts to data compression where a certain amount of
information may be lost. This dimension reduction is part of manifold
learning, which can be understood as a collection of algorithms for recovering
low dimension manifolds embedded into high dimensional ambient spaces, while
preserving meaningful information, see Ma and Fu (2012). The algorithms for
dimension reduction may be classified into linear and nonlinear methods or
parametric or nonparametric methods, where the goal is to select or extract
coarse features from high dimensional data. Among the pioneering linear
methods is the Principal Component Analysis (PCA) introduced by Hotelling
(1933). Its primary goal is to reduce the data to a set of orthogonal linear
projections ordered by decreasing variances. Another linear method is the
multidimensional scaling (MSD) where the data are aggregated using a measure
of proximity, which could be a distance, or a measure of association such as
correlation, or any other method describing how close entities can be. Linear
Discriminant Analysis (LDA) is a linear method similar to PCA consisting of
writing a categorical dependent variable as a linear combination of continuous
independent variables. As such, it is opposite to an Analysis of Variance
(ANOVA) where the dependent variable is continuous and the independent
variables are categorical. The focus of this paper will be on nonlinear
techniques, which as their linear counterparts, aim to extract or select low
dimensional features while preserving important information. Since there are
many such methods, our focus will be on Isometric Feature Mapping (ISOMAP),
Laplacian Eigenmaps, Fast Independent Component Analysis, (Fast-ICA),
T-distributed Stochastic Neighbor Embedding (t-SNE). We will compare them
using Persistent homology (PH). PH is one the many techniques of topological
data analysis (TDA) that can be used to identify features in data that remain
persistent over multiple and different scales. This tool can provide new
insights into seemingly known or unknown data and has the potential to uncover
interesting hidden information embedded within data. For instance, PH has been
used to provide new insights on the topology of deep neural networks, see
Naizait et al. (2020). PH has successfully been used to provide new
perspectives on viral evolution, see Chan et al. (2013). The following
examples of successful applications can be found in Otter et al. (2017)
including but not limited to better understanding of sensor-network coverage,
see de Silva and Ghrist (2007), proteins, see Xia and Wei (2014); Gameiro et
al. (2015), dimensional structure of DNA, see Emmett et al. (2016), see cells
development, Rizvi et al. (2017), robotics, see Bhattacharya et al. (2015);
Pokorny et al. (2016); Vasudevan et al. (2013), signal processing, see Chung
et al. (2009); Guillemard et al. (2013), spread of contagions, see Taylor et
al. (2015), financial networks, see Leibon et al. (2008), applications in
neuroscience, see Giusti et al. (2016); Sizemore et al. (2019), time-series
output of dynamical systems, see Maletić et al. (2016), on EEG Epileply, see
Chung et al. (2023). The approach is the last reference is of particular
interest to us. Indeed that paper, authors considered EEG measured on the
healthy person during sleep. They used the method of false nearest neighbors
to estimate the embedding dimension. From there, persistent barcodes diagrams
are obtained and reveal that topological noise persists at certain dimension
and vanish some others. This paper has a similar approach and is organized as
follows: in Section 2, we review theories behind some dimension reduction
methods, then in Section 3, we give an overview of the essentials of
persistent homology, in Section 4, we discuss how to apply persistent homology
to data, and compare the methods on an EEG dataset, using persistent homology.
Finally in Section 5, we make some concluding remarks.

##  2 Review of a selected dimension reduction methods

Let us note that some of the review methods below are extensively described in
Ma and Fu (2012). To have all of our ideas self-contained, let us re-
introduced a few concepts. In the sequel, $\left\|\cdot\right\|$ is the
euclidian norm in $\mathbb{R}^{d}$, for some $d\geq 3$. In the sequel,
topological spaces $\mathscr{M}$ will considered to be second-countable
Hausdorff, that is, (a) Every pair of distinct points has a corresponding pair
of disjoint neighborhoods. (b) Its topology has a countable basis of open
sets. This assumption is satisfied in most topological spaces and seems
reasonable.

###  2.1 Preliminaries

######  Definition 1.

A topological space $\mathscr{M}$ is called a (topological) manifold if
locally, it resembles a real $n$-dimensional Euclidian space, that is, there
exists $n\in\mathbb{N}$ such that for all $x\in\mathscr{M}$, there exists a
neighborhood $U_{x}$ of $x$ and a homeomorphism $f:U_{x}\to\mathbb{R}^{n}$.
The pair $(U_{x},f)$ is referred to as a chart on $\mathscr{M}$ and $f$ is
called a parametrization at $x$.

######  Definition 2.

Let $\mathscr{M}$ be a manifold. $\mathscr{M}$ is said to be smooth if given
$x\in\mathscr{M}$, the parametrization $f$ at $x$ has smooth or continuous
partial derivatives of any order and can be extended to a smooth function
$F:\mathscr{M}\to\mathbb{R}^{n}$ such that $F\big{|}_{\mathscr{M}\cap
U_{x}}=f$.

######  Definition 3.

Let $\mathscr{M}$ and $\mathscr{N}$ be differentiable manifolds and consider
$\psi:\mathscr{M}\to\mathscr{N}$ a function. $\psi$ is said to be an immersion
if $\psi$ is a differentiable function and its derivative is everywhere
injective. In other words, $\psi:\mathscr{M}\to\mathscr{N}$ is an immersion if
$D_{x}\psi:T_{x}\mathscr{M}\to T_{\psi(x)}\mathscr{N}$ is an injective
function at every point $x$ of $\mathscr{M}$, where $T_{x}\mathscr{M}$
represents the tangent plane to $\mathscr{M}$ at $x$.

######  Definition 4.

Let $\mathscr{M}$ and $\mathscr{N}$ be differentiable manifolds. A function
$\psi:\mathscr{M}\to\mathscr{N}$ is an embedding if $\psi$ is an injective
immersion.

Let us introduce the notion of boundary of topological manifold that will be
important in the sequel.

######  Definition 5.

Consider a Hausdorff topological manifold $\mathscr{M}$ homeomorphic to an
open subsets of the half-euclidian space $\mathbb{R}^{n}_{+}$. Let the
interior $Int(\mathscr{M})$ of $\mathscr{M}$ be the subspace of $\mathscr{M}$
formed by all points $s$ that have a neighborhood homeomorphic to
$\mathbb{R}^{n}$. Then the boundary of $\mathscr{M}$ is defined as complement
of $Int(\mathscr{M})$ in $\mathscr{M}$, that is, $\mathscr{M}\setminus
Int(\mathscr{M})$, which is an $n-1$-dimensional topological manifold.

###  2.2 Isomap

Isometric Feature Mapping (Isomap) was introduced by Tenenbaum et al. (2000).
The data are considered to be a finite sample $\left\\{\bm{v}_{i}\right\\}$
from a smooth manifold $\mathscr{M}$. The two key assumptions are: (a) there
exists of an isometric embedding $\psi:\mathscr{M}\to\mathscr{X}$ where
$\mathscr{X}=\mathbb{R}^{d}$, where the distance on $\mathscr{M}$ is the
geodesic distance or the shortest curve connecting two points, (b) the smooth
manifold $\mathscr{M}$ is a convex region of $\mathbb{R}^{m}$, where $m<<d$.
The implementation phase has three main steps.

  1. 1.

For a fixed integer $K$ and real number $\epsilon>0$, perform an
$\epsilon-K$-nearest neighbor search using the fact that the geodesic distance
$D^{\mathscr{M}}(v_{i},v_{j})$ between two points on $\mathscr{M}$ is the same
(by isometry) as their euclidian distance $\left\|v_{i}-v_{j}\right\|$ in
$\mathbb{R}^{d}$. $K$ is the number of data points selected within a ball of
radius $\epsilon$.

  2. 2.

Having calculated the distance between points as above, the entire data set
can be considered as a weighted graph with vertices
$\bm{v}=\left\\{v_{i}\right\\}$ and edges $\bm{e}=\left\\{e_{ij}\right\\}$,
where $e_{ij}$ connects $v_{i}$ with $v_{j}$ with a distance
$w_{ij}=D^{\mathscr{M}}(v_{i},v_{j})$ considered as an associated weight. The
geodesic distance between two data points $v_{i}$ and $v_{j}$ is estimated as
the graph distance between the two edges, that is, the number of edges in a
shortest path connecting them. We observe that this shortest path is found by
minimizing the sum of the weights of its constituent edges.

  3. 3.

Having calculated the geodesic distances $D^{G}=\left\\{w_{ij}\right\\}$ as
above, we observe that $D^{G}$ is a symmetric matrix, so we can apply the
classical Multidimensional Scaling algorithm (MDS) (see Torgerson (1952)) to
$D^{G}$ by mapping (embedding) them into a feature space $\mathscr{Y}$ of
dimension $d$ while preserving the geodesic distance on $\mathscr{M}$.
$\mathscr{Y}$ is generated by a $d\times m$ matrix whose $i$-th column
represents the coordinates of $v_{i}$ in $\mathscr{Y}$.

###  2.3 Laplacian Eigenmaps

The Laplacian Eigenmaps (LEIM) algorithm was introduced by Belkin and Niyogi
(2002). As above, the data ${\bm{v}}=\left\\{v_{i}\right\\}$ are supposed to
be from a smooth manifold $\mathscr{M}$. It also has three main steps:

  1. 1.

For a fixed integer $K$ and real number $\epsilon>0$, perform an
$\epsilon-K$-nearest neighbor search on symmetric neighborhoods. Note that
given two points $v_{i},v_{j}$, their respective $K$-neighborhood $N_{i}^{K}$
and $N_{j}^{K}$ are symmetric if and only $v_{i}\in
N_{j}^{K}\Longleftrightarrow v_{j}\in N_{i}^{K}$.

  2. 2.

For a given real number $\sigma>0$ and each pair of points $(v_{i},v_{j})$,
calculate the weight
$w_{ij}=e^{-\frac{\left\|u_{i}-v_{j}\right\|^{2}}{2\sigma^{2}}}$ if $v_{i}\in
N_{j}^{K}$ and $w_{ij}=0$ if $v_{i}\notin N_{j}^{K}$. Obtain the adjacency
matrix $\bm{W}=(w_{ij})$. The data now form a weighted graph with vertices
$\bm{v}$, with edges $\bm{e}=\left\\{e_{ij}\right\\}$, and weights
$\bm{W}=\left\\{w_{ij}\right\\}$, where $e_{ij}$ connects $v_{i}$ with $v_{j}$
with distance $w_{ij}$.

  3. 3.

Consider $\bm{\Lambda}=\left\\{\lambda_{ij}\right\\}$ be a diagonal matrix
with $\displaystyle\lambda_{ii}=\sum_{j}w_{ij}$ and define the graph Laplacian
as $\bm{L}=\bm{\Lambda}-\bm{W}$. Then $\bm{L}$ is positive definite so let
$\widehat{\bm{Y}}$ be the $d\times n$ matrix that minimizes
$\displaystyle\sum_{i,j}w_{ij}\left\|\bm{y}_{i}-\bm{y}_{j}\right\|^{2}=tr(\bm{TLY^{T}})$.
Then $\widehat{\bm{Y}}$ can used to embed $\mathscr{M}$ into a $d$-dimensional
space $\mathscr{Y}$, whose $i$-th column represents the coordinates of $v_{i}$
in $\mathscr{Y}$.

###  2.4 Fast ICA

The Fast Independent Component Analysis (Fast-ICA) algorithms were introduced
by Hyvärinen (1999). As above, the data $\bm{v}$ is considered to be from a
smooth manifold $\mathscr{M}$. It is assumed that the data $\bm{v}$ is
represented as an $n\times m$ matrix $(v_{ij})$ that can be flattened into a
$n\times m$ vector. As in Principal Component Analysis (PCA), in Factor
Analysis, Projection Pursuit, or Independent Component Analysis (ICA), by
considering the data as an $n\times m$-dimensional observed random variable,
the goal is to determine a matrix $\bf{W}$ such that $\bm{s=W^{T}v}$, where
$\bm{s}$ is a $n\times m$-dimensional random variables having desirable
properties such as optimal dimension reduction, or other interesting
statistical properties such as minimal variance. Optimally, the components of
$\bm{s}$ should provide source separation (the original data source $\bm{v}$
is assumed corrupted with noise) and feature extraction and be independent of
each other. In a regular ICA, the matrix $\bm{W}$ is found by minimizing the
mutual information, a measure of dependence between given random variables. In
fast ICA algorithms, the matrix $\bm{W}$ is found by using a Newton fixed
point approach, with an objective function taken as the differential entropy
given as
$J_{G}(\bm{W})=\left(\mathbb{E}[G(\bm{W^{T}W})]-\mathbb{E}[G(z)]\right)^{2}$,
where it is assumed that $\bm{W}$ is such that
$\mathbb{E}[\bm{(W^{T}W)}^{2}]=1$, and $z$ is standard normal distribution.
$G$ is a function referred to as the contrast function that include but is not
limited to $G(u)=\alpha^{-1}\log(\cosh(\alpha
u)),G(u)=-\sigma^{-1}e^{-0.5\sigma u^{2}},G(u)=0.25u^{4}$, where
$\alpha\in[1,2]$ and $\sigma\approx 1$. From a dynamical system point of view,
the fixed point is locally asymptotically stable with the exception of
$G(u)=0.25u^{4}$ where stability becomes global. For simplification purposes,
let $g(x)=G^{\prime}(x)$. The key steps are:

  1. 1.

Data preparation: it consists of centering the data $\bm{v}$ with respect to
the column to obtain $\bm{v}^{c}$. That is, $\displaystyle
v^{c}_{ij}=v_{ij}-\frac{1}{m}\sum_{j=1}^{m}v_{ij}$, for $i=1,2,\cdots,n$. The
centered data are then whitened, that is, $\bm{v}^{c}$ is linearly transformed
into $\bm{v}_{w}^{c}$, a matrix of uncorrelated components. This is
accomplished through an eigenvalue decomposition of the covariance matrix
$\bm{C=v^{c}(v^{c})^{T}}$ to obtain two matrices $\bm{V},\bm{E}$, respectively
of eigenvectors and eigenvalues so that $\bm{\mathbb{E}[C]=VEV^{T}}$. The
whitened data are found as $\bm{v^{c}_{w}=E^{-1/2}V^{T}v^{c}}$ and simply
referred to again as $\bm{v}$ for simplicity.

  2. 2.

Component extraction: Let
$F(\bm{W})=\mathbb{E}[\bm{v}g(\bm{W^{T}v})]-\beta\bm{W}$ for a given constant
$\beta=\mathbb{E}[\bm{W_{a}^{T}v}g(\bm{W_{a}^{T}v})]$, where ${W_{a}}$ is the
optimal weight matrix. Applying the Newton scheme
($x_{n+1}=x_{n}-F(x_{n})[F^{\prime}(x_{n})]^{-1}$) to the differentiable
function $J_{G}$, we have

     * •

Select a random starting vector $\bm{W}_{0}$.

     * •

For $n\geq 0$,
$\bm{W}_{n+1}=\mathbb{E}[\bm{v}g(\bm{W}_{n}^{T}\bm{v})]-\mathbb{E}[g^{\prime}(\bm{W}_{n}^{T}\bm{v})]\bm{W}_{n}$.

     * •

We then normalize $\bm{W}_{n+1}$ as
$\displaystyle\frac{\bm{W}_{n+1}}{\left\|\bm{W}_{n+1}\right\|}$.

     * •

We repeat until a suitable convergence level is reached.

     * •

From the last matrix $\bm{W}$ obtained, we let $\bm{s=W^{T}v}$.

###  2.5 Kernel Ridge Regression

The Kernel Ridge Regression (KRR) is constructed as follows: as above, the
data $\bm{v}$ is considered to be from a smooth manifold $\mathscr{M}$ of
dimension say $d$. It is assumed that the data $\bm{v}$ is represented as an
$n\times m$ matrix $\left\\{v_{ij}\right\\}$ that can be flattened into a
$n\times m$ vector. Suppose we are in possession of
$\bm{u}=(u_{1},u_{2},\cdots,u_{n})$ data corresponding to a response variable
and covariates given as $\bm{v}=(\bm{v_{1},v_{2},\cdots,v_{n}})$ where
$\bm{v_{i}}=(v_{ij})^{T}$ for $j=1,2,\cdots,m$. With the Least Square method,
on can find the best linear model between the covariates $\bm{v}=(v_{i})$ and
the response $\bm{u}=(u_{i})$, by minimizing the objective function
$\displaystyle
L(\bm{W})=\frac{1}{2}\sum_{i=1}^{L}(u_{i}-\bm{W^{T}v}_{i})^{2}$, where
$\bm{W}$ is a $1\times n$ vector. Least square methods are notorious for
overfitting. The Ridge regression is a compromise that uses a penalized
objective function such as $\displaystyle
L(\bm{W})=\frac{1}{2}\sum_{i=1}^{L}(u_{i}-\bm{W^{T}v}_{i})^{2}+\frac{\lambda}{2}\left\|\bm{W}\right\|^{2}$.
The solution can be found as $\displaystyle\bm{W}=\left(\lambda
I+\sum_{i=1}^{n}\bm{v_{i}v_{i}}^{T}\right)^{-1}\left(\sum_{i=1}^{n}u_{i}\bm{v}_{i}\right)$.
In case the true nature of the relationship between the response and
covariates is nonlinear, we can replace $\bm{v}_{i}$ with
$\varphi(\bm{v}_{i})$ where $\varphi$ is a nonlinear function
$\mathbb{R}^{m}\to\mathbb{R}$. In particular, if the response is qualitative,
that is, say labels, then we have a classification problem and $\varphi$ is
referred to as feature map. Note that when using $\varphi$, the number of
dimensions of the problem is considerably high. Put
$\Phi=\varphi(\bm{v})=(\varphi(v_{1}),\varphi(v_{2}),\cdots,\varphi(v_{n}))$.
Replacing $\bm{v_{i}}$ with $\varphi(\bm{v_{i}})$, the solution above becomes
$\displaystyle\bm{W}=\left(\lambda
I+\sum_{i=1}^{n}\varphi(\bm{v_{i}})\varphi(\bm{v_{i}})^{T}\right)^{-1}\left(\sum_{i=1}^{n}u_{i}\varphi(\bm{v}_{i})\right)=(\lambda
I+\Phi\Phi^{T})^{-1}\Phi\bm{u}^{T}$. Consider the following identity
$AB^{T}(C+BAB^{T})^{-1}=(A^{-1}+B^{T}C^{-1}B)^{-1}B^{T}C^{-1}$ for given
invertible matrices $A,C$ and a matrix $B$. Applying this with $A=C=I$ and
$B=\Phi$ we have $\bm{W^{T}}=\bm{u}\left[\Phi^{T}(\lambda
I+\Phi^{T}\Phi)^{-1}\right]=\bm{u}\left[(\lambda
I+\Phi^{T}\Phi)^{-1}\Phi^{T}\right]$. Therefore, given a new value
$\bm{v}_{n}$, the predicted value is
$y_{n}=\bm{W}^{T}\Phi(\bm{v}_{n})=\bm{u}(\Phi^{T}\Phi+\lambda
I)^{-1}\Phi^{T}\Phi(\bm{v}_{n})=\bm{u}(K+\lambda I)^{-1}\kappa(\bm{v}_{n})$,
where $\displaystyle
K=K(\bm{v}_{i},\bm{v}_{i})=\Phi^{T}\Phi=\sum_{i=1}^{n}\varphi(\bm{v}_{i})^{T}\varphi(\bm{v}_{i})$
and $\kappa(\bm{v}_{n})=K(\bm{v}_{i},\bm{v}_{n})$. K is referred to as the
kernel, which is the only quantity needed to be calculated, thereby
significantly reducing the computational time and dimensionality of the
problem. In practice, we may use a linear kernel
$K(\bm{x},\bm{y})=\bm{x}^{T}\bm{y}$, or a Gaussian kernel
$K(\bm{x},\bm{y})=e^{-\sigma\left\|\bm{x}-\bm{y}\right\|^{2}}$, for some
$\sigma>0$, where $\left\|\cdot\right\|$ is a norm in $\mathbb{R}^{m}$ and
$\sigma$ is given real constant.

###  2.6 t-SNE

Stochastic Neighbor Embedding (SNE) was proposed by Hinton and Roweis (2002).
t-SNE latter followed and was proposed by van der Maaten and Hinton (2008).
t-distributed stochastic neighbor embedding (t-SNE) is a dimension reduction
method that amounts to assigning data two or three dimensional maps . As
above, we consider the data $\bm{v}=(v_{ij})=(v_{k})$ ($k=1,2,\cdots,N$ with
$N=n\times m$) to be from a smooth manifold $\mathscr{M}$ of high dimension,
say $d$. The main steps of the method are:

  * •

Calculate the asymmetrical probabilities $p_{kl}$ as
$p_{kl}=\frac{e^{-\delta_{kl}}}{\sum_{k\neq l}e^{-\delta_{kl}}}$, where
$\delta_{kl}=\frac{\left\|v_{k}-v_{l}\right\|^{2}}{2\sigma_{i}}$ represents
the dissimilarity between $v_{k}$ and $v_{l}$ and $\sigma_{i}$ is a parameter
selected by the experimenter or by binary search. $p_{kl}$ represents the
conditional probability that datapoint $v_{l}$ is the neighborhood of
datapoint $v_{k}$, if neighbors were selected proportionally to their
probability density under a normal distribution centered at $v_{k}$ and
variance $\sigma_{i}$.

  * •

Assuming that the low dimensional data are
$\bm{u}=(u_{k}),\leavevmode\nobreak\ k=1,2,\cdots,N$, the corresponding
dissimilarity probabilities $q_{kl}$ are calculated under constant variance as
$q_{kl}=\frac{e^{-d_{kl}}}{\sum_{k\neq l}e^{-d_{kl}}}$, where
$d_{kl}=\left\|u_{k}-u_{l}\right\|^{2}$ in the case of SNE and as
$\displaystyle q_{kl}=\frac{(1+d_{kl})^{-1}}{\sum_{k\neq
l}{(1+d_{kl})^{-1}}}$, for t-SNE.

  * •

Then, we minimize the Kullback-Leibler divergence between $p_{kl}$ and
$q_{kl}$ given as $\displaystyle
L=\sum_{k=1}^{N}\sum_{l=1}^{N}p_{kl}\log\left(\frac{p_{kl}}{q_{kl}}\right)$,
using the gradient descent method with a momentum term with the scheme
$\displaystyle\bm{w}^{t}=\bm{w}^{t-1}+\eta\frac{\partial
L}{\partial\bm{u}}+\alpha(t)(\bm{w}^{t-1}-\bm{w}^{t-2})$ for $t=2,3,\cdots,T$
for some given $T$. Note that $\bm{w}^{0}=(u_{1},u_{2},\cdots,u_{N})\sim
N(0,10^{-4}\bm{I})$, where $\bm{I}$ is the $N\times N$ identity matrix, $\eta$
is a constant representing a learning rate, and $\alpha(t)$ is $t$-th momentum
iteration. We note that $\displaystyle\frac{\partial
L}{\partial\bm{u}}=\left(\frac{\partial L}{\partial u_{k}}\right)$ for
$k=1,2,\cdots,N$ where $\displaystyle\frac{\partial L}{\partial
u_{k}}=4\sum_{l=1}^{N}(p_{kl}-q_{kl})(u_{k}-u_{l})(1+d_{kl})^{-1}.$

  * •

Then we use $\bm{u}=\bm{w}^{T}$ as the low dimensional representation of
$\bm{v}$.

##  3 Persistent Homology

In the sequel, we will introduce the essential ingredients needed to
understand and compute persistent homology.

###  3.1 Simplex complex

######  Definition 6.

A real $d$-simplex $S$ is a topological manifold of dimension $d$ that
represents the convex hull of $d+1$ points. In other words:

$S=\left\\{(t_{0},t_{1},\cdots,t_{d})\in\mathbb{R}^{d}:t_{i}\geq
0\quad\mbox{and $\displaystyle\sum_{i=1}^{d}t_{i}=1$}\right\\}\;.$ (3.1)

Example 1.  A $0$-simplex is a point, a $1$-simplex is an edge, a $2$-simplex
is a triangle, a $3$-simplex is a tetahedron, a $4$-simplex is a pentachoron,
etc.

Figure 1: An illustration of 0, 1, 2, 3, and 4-simplices.

######  Remark 7.

We observe that a $d$-simplex $S$ can also be denoted as

$S=[V_{0},V_{1},\cdots,V_{d}],\quad\mbox{where $V_{i}=\left\\{\mbox{vertices
of $V_{i}$ }\right\\},i=0,1,\cdots,d$}\;.$

We also note that the dimension of $V_{i}$ is $i$.

######  Definition 8.

Given a simplex $S$, a face of $S$ is another simplex $R$ such that
$R\subseteq S$ and such that the vertices of $R$ also the vertices of $S$.

Example 2.  Given a 3-simplex (a tetrahedron), it has 4 different 2-simplex or
2 dimensional faces, each of them with three 1-simplex or 1-dimensional faces,
each with three 0-simplex or 0-dimensional faces.

######  Definition 9.

A simplicial complex $\Sigma$ is a topological space formed by different
simplices not necessarily of the same dimension which have to satisfy the
gluing condition, that is:

  1. 1.

Given $S_{i}\in\Sigma$, its face $R_{i}\in\Sigma$.

  2. 2.

Given $S_{i},S_{j}\in\Sigma$, either $S_{i}\cap S_{j}=\emptyset$ or $S_{i}\cap
S_{j}=R_{i}=R_{j}$, the faces of $S_{i}$ and $S_{j}$ respectively.

Figure 2: Example of a simplicial complex. $J$ is 0-simplex, $A$ and $D$ are
1-simplices, $B,C,G$, and $H$ are 2-simplices, $E$ and $F$ are 3-simplices,
and $I$ is a 4-simplex. We note that $A\cap B$ is a 0-simplex. $B\cap C$ is a
1-simplex and a face of $B$ and $C$ respectively. $E\cap F$ is a 2-simplex and
a face of $E$ and $F$. $G\cap H$ is a 1-simplex and $I\cap H$ is 1-simplex.

It is important to observe that a simplicial complex can be defined very
abstractly. Indeed,

######  Definition 10.

A simplicial complex $\Sigma=\left\\{S:S\subseteq\Omega\right\\}$ is a
collection of non-empty subsets of a set $\Omega$ such that

  1. 1.

For all $\omega\in\Omega$, then $\left\\{\omega\right\\}\in\Sigma$.

  2. 2.

For any set $U$ such that $S\subset U$ for some $S\in\Sigma$, then
$U\in\Sigma$.

Example 3.  Let $\Omega=\left\\{1,2,3,4\right\\}$. We can define the following
simplicial complexes on $\Omega$.
1\.
$\displaystyle\Sigma_{1}=\left\\{\left\\{1\right\\},\left\\{2\right\\},\left\\{3\right\\},\left\\{4\right\\},\left\\{1,2\right\\},\left\\{1,3\right\\},\left\\{2,3\right\\},\left\\{1,2,3\right\\}\right\\}\;.$
2\. $\Sigma_{2}=\mathscr{P}(\Omega)\setminus\left\\{\emptyset\right\\}$, where
$\mathscr{P}(\Omega)$ is the set of all subsets of $\Omega$.

###  3.2 Homology and persistent homology

######  Definition 11.

Let $\Sigma$ be a simplicial complex.
We define the Abelian group generated by the $j$-simplices of $\Sigma$ as
$C_{j}(\Sigma)$.
We define a boundary operator associated with $C_{j}(\Sigma)$ as a
homomorphism

$\partial_{j}:C_{j}(\Sigma)\to C_{j-1}(\Sigma)\;.$

We define the chain complex associated with $\Sigma$ as the collection of
pairs

$C(\Sigma)=\left\\{(C_{j}(\Sigma),\partial_{j}),j\in\mathbb{Z}\right\\}\;.$

Now we can define a homology group associated with a simplicial complex.

######  Definition 12.

Given a simplicial complex $\Sigma$, put $A_{j}(\Sigma):=kern(\partial_{j})$
and $B_{j}(\Sigma):=Im(\partial_{j+1})$. Then the $j$th homology group
$\mathbb{H}_{j}(\Sigma)$of $\Sigma$ is defined as quotient group between
$A_{j}(\Sigma)$ and $B_{j}(\Sigma)$, that is,

$\mathbb{H}_{j}(\Sigma)=\frac{A_{j}(\Sigma)}{B_{j}(\Sigma)}\;.$

What this reveals is the presence of “holes” in a given shape.

######  Remark 13.

It is important to observe that
$\displaystyle\mathbb{H}_{j}(\Sigma)=\frac{<\mbox{$j$-dimensional
cycles}>}{<\mbox{$j$-dimensional boundaries}>}$, where $<U>$ stands for the
span of $U$, and a cycle is simply a shape similar to a loop but without
necessarily a starting point.
Another important remark is that the boundary operator can indeed be defined
as

$\partial_{j}(\Sigma):=\sum_{k=0}^{j}(-1)^{k}[V_{0},\cdots,\widehat{V}_{-k},\cdots,V_{j}]\;,$

where $\widehat{V}_{-k}$ means not counting the vertices of $V_{k}$. This
shows that $\partial_{j}(\Sigma)$ lies in a $j-1$-simplex.
Another remark is that $\partial_{j-1}\circ\partial_{j}=0$ for $0\leq j\leq d$
.

Now that we know that homology reveals the presence of “holes”, we need to
find a way of determining how to count these “holes”.

######  Definition 14.

Given a simplicial complex $\Sigma$, the $j$th Betti number $b_{j}(\Sigma)$ is
the rank of $\mathbb{H}_{j}(\Sigma)$ or

$b_{j}(\Sigma)=dim(A_{j}(\Sigma))-dim(B_{j}(\Sigma))\;.$

In other words, the smallest cardinality of a generating set of the group
$\mathbb{H}_{j}(\Sigma)$.
In fact since the elements of $A_{j}(\Sigma)$ are $j$-dimensional cycles and
that of $B_{j}(\Sigma)$ are $j$-dimensional boundaries, the Betti number
counts the number of independent $j$-cycles not representing the boundary of
any collection of simplices of $\Sigma$.

Example 4.

  1. 1.

$b_{0}$ is the number of connected components of the complex.

  2. 2.

$b_{1}$ is the number of tunnels and holes.

  3. 3.

$b_{2}$ is the number of shells around cavities or voids.

######  Definition 15.

Let $\Sigma$ be simplicial complex and let $N$ be a positive integer. A
filtration of $\Sigma$ is a a nested family
$\Sigma^{F}_{N}:=\left\\{\Sigma_{p},0\leq p\leq N\right\\}$ of sub-complexes
of $\Sigma$ such that

$\Sigma_{0}\subseteq\Sigma_{1}\subseteq\Sigma_{2}\subseteq\cdots\subseteq\Sigma_{N}=\Sigma\;.$

Now let $\mathbb{F}_{2}$ be the field with two elements and let $0\leq p\leq
q\leq N$ be two integers. Since $\Sigma_{p}\subseteq\Sigma_{q}$, the inclusion
map $Incl_{pq}:\Sigma_{p}\to\Sigma_{q}$ induces an $\mathbb{F}_{2}$-linear map
$g_{pq}:\mathbb{H}_{j}(\Sigma_{p})\to\mathbb{H}_{j}(\Sigma_{q})$. We can now
define, for any $0\leq j\leq d$, the $j$-th persistent homology of a
simplicial complex $\Sigma$.

######  Definition 16.

Consider a simplicial complex $\Sigma$ with filtration $\Sigma^{F}_{N}$, for
some positive integer $N$. The $j$-th persistent homology
$\mathbb{H}_{j}^{p\to q}(\Sigma)$ of $\Sigma$ is defined as the pair:

$\mathbb{H}_{j}^{p\to
q}(\Sigma,\mathbb{F}_{2}):=(\left\\{\mathbb{H}_{j}(\Sigma_{p}),0\leq p\leq
N\right\\},\left\\{g_{pq},0\leq p\leq q\leq N\right\\})\;.$

In a sense, the $j$-th persistent homology provides a more refined information
than the homology of the simplicial complex in that it informs us of the
changes of features such as connected components, tunnels and holes, shells
around voids through the filtration process. It can be visualized using a
“barcode” or a persistent diagram. The following definition is borrowed from
Ghrist (2008):

######  Definition 17.

Consider a simplicial complex $\Sigma$, a positive integer $N$, and two
integers $0\leq p\leq q\leq N$. The barcode of the $j$-th persistent homology
$\mathbb{H}_{j}^{p\to q}(\Sigma,\mathbb{F}_{2})$ of $\Sigma$ is a graphical
representation of $\mathbb{H}_{j}^{p\to q}(\Sigma,\mathbb{F}_{2})$ as a
collection of horizontal line segments in a plane whose horizontal axis
corresponds to a parameter and whose vertical axis represents an arbitrary
ordering of homology generators.

We finish this section with the introduction of the Wasserstein and Bottleneck
distances, used for the comparison of persistent diagrams.

######  Definition 18.

Let $p>1$ be a real number. Given two persistent diagrams $X$ and $Y$, the
$p$-th Wasserstein distance $W_{p}(X,Y)$ between $X$ and $Y$ is defined as

$W_{p}(X,Y):=\inf_{\eta:X\to Y}\sum_{x\in
X}\left\|x-\eta(x)\right\|_{\infty}^{p}\;,$

where $\eta$ is a perfect matching between the intervals of $X$ and $Y$.
The Bottleneck distance is obtained when $p=\infty$, that is, it is given as

$W_{\infty}(X,Y):=\inf_{\eta:X\to Y}\sup_{x\in
X}\left\|x-\eta(x)\right\|_{\infty}\;.$

##  4 Applications to Data

In the presence of data, simplicial complexes will be replaced by sets of data
indexed by a parameter, therefore transforming these sets into parametrized
topological entities. On these parametrized topological entities, the notions
of persistent homology introduced above can be computed, especially the Betti
number, in the form of “barcode”. To see how this could be done, let us
consider the following definitions:

######  Definition 19.

For a given collection of points $\left\\{s_{\delta}\right\\}$ in a manifold
$\mathscr{M}$ of dimension $n$, its Čech complex $C_{\delta}$ is a simplicial
complex formed by $d$-simplices obtained from a sub-collection
$\left\\{x_{\delta,k},0\leq k\leq d,0\leq d\leq n\right\\}$ of points such
that taken pairwise, their $\delta/2$-ball neighborhoods have a point in
common.

######  Definition 20.

For a given collection of points $\left\\{s_{\delta}\right\\}$ in a manifold
$\mathscr{M}$ of dimension $n$, its Rips complex $R_{\delta}$ is a simplicial
complex formed by $d$-simplices obtained from a sub-collection
$\left\\{x_{\delta,k},0\leq k\leq d,0\leq d\leq n\right\\}$ of points which
are pairwise within a distance of $\delta$.

######  Remark 21.

1\. It is worth noting that in practice, Rips complexes are easier to compute
than Čech complexes, because the exact definition of the distance on
$\mathscr{M}$ may not be known.
2\. More importantly, from a data analysis point of view, Rips complexes are
good approximations (estimators) of Čech complexes. Indeed, a result from de
Silva and Ghrist (2007) shows that given $\delta>0$, there exits a chain of
inclusions $R_{\delta}\hookrightarrow C_{\delta/\sqrt{2}}\hookrightarrow
R_{\delta/\sqrt{2}}$. 3\. Though Rips complexes and barcodes seems like
challenging objects to wrap one’s head around, there is an ever growing list
of algorithms from various languages that can be used for their visualization.
All the analysis below has been done using R, in particular the TDA package in
R.

###  4.1 Randomly generated data

We generated 100 data points sampled randomly in the square
$[-5,5]\times[-5,5]$. In Figure 3 and Figure 4 below, we illustrate the Rips
and barcode changes through a filtration.

Figure 3: Example of the evolution of Rips Complexes
$\left\\{R_{\delta}\right\\}$ through a filtration with parameter $\delta$. As
we move from left to right, it shows how samples points (blue dots), first
form 0-simplices, then 1-simplices, and so on. In particular, it shows how
connected components progressively evolve to form different types of holes.
Figure 4: Example of the evolution of barcodes through a filtration with
parameter $\delta$ for the same data as above. As we move from left to right,
from top to bottom, it shows appearance and disappearance of lines
($\mathbb{H}_{0}$) and holes ($\mathbb{H}_{1}$) as the parameter $\delta$
changes. It shows that certain lines and holes persist through the filtration
process.

###  4.2 EEG epilepsy data

####  4.2.1 Data Description

The main purpose of the manuscript it to analyze EEG data. We will consider a
publicly available (at
http://www.meb.unibonn.de/epileptologie/science/physik/eegdata.html) epilepsy
data set called here EDATA for simplicity. The data consist of five sets A, B,
C, D, and E. Each containing 100 single-channel EEG segments of 23.6 seconds,
each of which was selected after visual inspection for artifacts (such as
acoustic and electrical shielding, separation of earth ground for laboratory,
interconnectivity of devices on the same phase and ground centrally and
locally) and has passed a weak stationarity criterion. Sets A and B were
obtained from surface EEG recordings of five healthy subjects with eyes open
and closed, respectively. Data were obtained in seizure-free intervals from
five patients in the epileptogenic zone for set D and from the hippocampal
formation of the opposite hemisphere of the brain for set C. Set E contains
seizure activity, selected from all recording sites exhibiting ictal activity.
Sets A and B have been recorded extracranially, whereas sets C, D, and E have
been recorded intracranially. All EEG signals were recorded with the same
128-channel amplifier system, using an average common reference omitting
electrodes containing pathological activity (C,D, and E) or strong eye
movement artifacts (A and B). After 12 bit analog-to-digital conversion, the
data were written continuously onto the disk of a data acquisition computer
system at a sampling rate of 173.61 Hz. Band-pass filtersettings were 0.53–40
Hz (12 dB/oct.)

####  4.2.2 Data analysis

The approach is to first embed the data into a manifold of high dimension.
This was already done in Kwessi and Edwards (2021). The dimension $d=12$ was
found using the method of false nearest neighbors. Depending on the set used,
the size of the data can be very large: for example ($4097\times 100\times
5=2,048,500$) making it very challenging to analyze holistically. In Kwessi
and Edwards (2021), we proposed to construct a complex structure (using all
100 channels for all 5 groups) whose volume changes per group. We would like
to analyze the data further from a persistent homology point of view. This
would mean analyzing 500 different persistent diagrams and making an
inference. We note that simplicial complexes of this data sets are very large
(2 Millions+). Fortunately, we can use the Wasserstein distance to compare
persistent diagrams. To clarify, we will use each of the dimension reduction
method introduced earlier, then proceed with construction of persistent
diagrams. We will then compare them by method and by sets (A, B, C, D, and E).

Single-channel Analysis:

Suppose we select at random one channel among the 100 from set D. Figure 5
below represents a 3 dimensional representation of the embedded data using
Takens embedding method (Tak), plotted using the 3 first three delayed
coordinates $x=x(t),y=x(t-\rho),z=x(t-2\rho)$ where $\rho=1\Delta t$, with
$\Delta t=\frac{1}{fs}=5.76$ ms in Figure 5 (a), then the first three
coordinates in the case of Kernel Ridge Regression (KRRi) in Figure 5 (b),
Isomap (iso.$i$) in Figure 5 (c), Laplacian Eigen Maps (LEIMi) in Figure 5
(d), Fast Independent Component Analysis (ICAi) in Figure 5 (e), and
t-Distributed Stochastic Neighbor Embedding (t-SNEi) Figure 5 (f). From these
three dimensional scatter plots, we can visually observe that the t-SNE plot
(Figure 5 (f)) is relatively different from the other five since it seems to
have more larger voids. How different is difficult to tell from the naked-eye.
Figure 6 represents their corresponding barcodes. It is much clearer looking
at the the persistent diagram for t-SNE (Figure 6 (f)) that it is very
different from the other five, when looking at $\mathbb{H}_{0},\mathbb{H}_{1}$
and $\mathbb{H}_{2}$. Now, a visual comparison is not enough to really assert
a significant difference. Using the Bottleneck distance, we calculate the
distance between the respective persistent diagrams, for $\mathbb{H}_{0}$ and
$\mathbb{H}_{1}$ in Table 1 (a) and $\mathbb{H}_{2}$ in Table 1 (b) below. We
observe from the first table that the Bottleneck distance at $\mathbb{H}_{0}$
and $\mathbb{H}_{2}$ for t-SNE are almost twice as large as for the other
methods. They are comparable to that of LEIM at $\mathbb{H}_{1}$. The other
methods have comparable Bottleneck distances at
$\mathbb{H}_{0},\mathbb{H}_{1}$, and $\mathbb{H}_{2}$, confirming what we
already suspected visually in Figure 5 and Figure 6.

Figure 5: Scatterplots for a Takens projection method (a), KRR method (b),
Isomap (c), LEIM (d), ICA (e), and t-SNE (f). Figure 6: Barcodes for a Takens
projection method(a), KRR method (b), Isomap (c), LEIM (d), ICA (e), and t-SNE
(f).

(a) | (b)
---|---
| $\mathbb{H}_{0}$ | Tak | Iso | KRR | ICA | LEIM | TSNE
---|---|---|---|---|---|---
Tak |  |  |  |  |  |
Iso | 0.0945019 |  |  |  |  |
KRR | 0.0957546 | 0.0200035 |  |  |  |
ICA | 0.0982795 | 0.0157002 | 0.0071899 |  |  |
LEIM | 0.1678820 | 0.1182656 | 0.1247918 | 0.1205499 |  |
TSNE | 0.2238167 | 0.1730406 | 0.1817924 | 0.1759454 | 0.1162392 |
|  $\mathbb{H}_{1}$/$\mathbb{H}_{2}$ | Tak | Iso | KRR | ICA | LEIM | TSNE
---|---|---|---|---|---|---
Tak |  | 0.0363205 | 0.0301992 | 0.0292631 | 0.0291247 | 0.0551774
Iso | 0.0340282 |  | 0.0330687 | 0.0290406 | 0.0236890 | 0.0598517
KRR | 0.0317261 | 0.0279460 |  | 0.0207599 | 0.0212138 | 0.0647935
ICA | 0.0310771 | 0.0270919 | 0.0208086 |  | 0.0242277 | 0.0611090
LEIM | 0.0607389 | 0.0725585 | 0.0702695 | 0.0682761 |  | 0.0542615
TSNE | 0.0757815 | 0.0959521 | 0.0864587 | 0.0861522 | 0.0785030 |

Table 1: Bottleneck distance between the persistent diagrams above at
$\mathbb{H}_{0}$ (black), at $\mathbb{H}_{1}$ (blue), and at $\mathbb{H}_{2}$
(red) .

The analysis above was done using a single channel, selected at random from
the set D. It seems to suggest that the t-SNE method is different from the
other five dimension reduction methods discussed above. Strictly speaking, non
zero Bottleneck distances are indication of structural topological
differences. What they do not say however is, if the differences observed are
significant. To address the issue of significance, we will perform a pairwise
permutation test. Practically, from set $j$ and channel $i$, we will obtain a
persistent diagram $\mathscr{D}^{(j)}_{i}\sim\mathscr{P}^{(j)}$ where
$j\in\left\\{1,2,3,4,5\right\\},i\in\left\\{1,2,\cdots,15\right\\}$, and
$\mathscr{P}^{(j)}$ is the true underlying distribution of persistent
diagrams, see Mileyko et al. (2011) for the existence of these distributions.
We will conduct a pairwise permutation test with null hypothesis
$H_{0}:\mathscr{P}^{(j)}=\mathscr{P}^{(j^{\prime})}$ and alternative
hypothesis $H_{1}:\mathscr{P}^{(j)}\neq\mathscr{P}^{(j^{\prime})}$. We will
use landscape functions (see Berry et al. (2020)) to obtain test statistics.
The p.values obtained were found to be very small, suggesting that the
differences above are indeed all significative across
$\mathbb{H}_{0},\mathbb{H}_{1}$, and $\mathbb{H}_{2}$.

Multiple-channel Analysis:

(a) Within set analysis

In each set, we make a random selection of 15 channels, and we compare the
Bottleneck distances obtained. This means having 15 tables of distance such as
Table 1 (b) above. There will be consistency if the cell value $k(i,j)$ in
Table $k$, where $k\in\left\\{1,2,\cdots,15\right\\}$ and
$i,j\in\left\\{1,2,3,4,5\right\\}$ is barely different from $k^{\prime}(i,j)$
of Table $k^{\prime}$. Large differences will be an indication of topological
differences between the methods within the sets. In Figure 7 below, the
$y$-axis represents Bottleneck distances and the $x$-axis represents channels
indices. The red color is indicative of the Bottleneck distance between
persistent diagrams on $\mathbb{H}_{1}$ and the blue color on $\mathbb{H}_{2}$
from data generated from each of the methods above. We see that overall, while
there are small fluctuations from channels to channels on $\mathbb{H}_{1}$,
the largest fluctuations actually occur on $\mathbb{H}_{2}$. A deeper analysis
reveals that in fact, the large fluctuations are due to large distance between
t-SNE and the other five methods. This confirms the earlier observations
(refer to Figure 6 and Table 1 above) that persistent diagrams are really
different on $\mathbb{H}_{2}$. Topologically, this means that shells around
cavities or voids that persist are not the same when using different dimension
reduction methods. However, the small fluctuations on $\mathbb{H}_{1}$ do not
mean that tunnels and holes that persist are the same. Rather, they do
indicate is that they may not be all very different.

(A) | (B) | (C)
---|---|---
|  |
(D) | (E) | (C)
|  |

Figure 7: Bottleneck distances between the persistent diagrams for 15 channels
within each set A, B, C, D, and E on $\mathbb{H}_{1}$ and $\mathbb{H}_{2}$ for
each of the methods introduced above.

(b) Between set analysis

To analyze the data of Bottleneck distances between sets, we need a summary
statistics for each set from the data above. It is clear from Figure 7 that
the mean would not be a great summary statistics on $\mathbb{H}_{1}$, as there
seems to be too many outliers. We will use the median instead and perform a
pairwise Wilcoxon-Mann-Whitney test. Table 2 below shows the p.value on
$\mathbb{H}_{1}$ and and $\mathbb{H}_{2}$. The take-away is that the last row
of the table suggests that set E is statistically topologically different from
others on $\mathbb{H}_{1}$, at significance level $0.05$. In a way, this is a
confirmation of the results obtained in Kwessi and Edwards (2021) where set E
(seizure) was already shown to be statistically different from other sets.

$\mathbb{H}_{1}$/$\mathbb{H}_{2}$ | A | B | C | D | E
---|---|---|---|---|---
A |  | 0.1975936 | 0.3049497 | 0.2467548 | 0.7432987
B | 0.3202554 |  | 0.3835209 | 0.5066311 | 0.1707835
C | 0.0832231 | 0.1322987 |  | 0.8356690 | 0.7088614
D | 0.2012797 | 0.6292608 | 0.6292608 |  | 0.5067258
E | 0.0049325 | 0.0157855 | 0.0157855 | 0.0114901 |
Table 2: P.values of Wilcoxon-Mann-Whitney tests between sets of median
Bottleneck distances.

##  5 Concluding remarks

In this paper, we have revisited the mathematical descriptions of six
dimension reduction methods. We have given a brief introduction to the very
vast topic of persistent homology. We discussed how to apply persistent
homology to data. In the presence of data (say in three dimension) obtained
either by projecting the data from high dimension into smaller dimension (as
in Takens) or by performing some sort of dimension reduction, it is not always
clear what we see or how different one method is compared to another. From
their mathematical description, they seem to represent different objects.
Further, obtaining theoretically a clear discrimination procedure between
these procedures seems a daunting if not an outright impossible task. Topology
may offer a solution by looking at persistent artifacts through filtration.
From Figure 5, it seems clear that the methods were different but Figure 6
offers a different perspective. In the end, through calculation of Bottleneck
distances and hypothesis tests, we can safely conclude that the methods are
different topologically speaking, in that, the connected components, the
tunnels and holes, the shells around cavities or voids do not match perfectly.
Since these methods are indiscriminately used in many applications, the
message is that replication of results from one method to the next may not be
guaranteed in the grand scheme of things. It does not however render them
useless. In fact, our analysis is limited to one data set, meaning that
another data set may yield different conclusions. Further, due to cost in
calculation, we were limited to only a handful of samples. More, Wasserstein
distance for $p<\infty$ are extremely costly in time to calculate on a regular
computer. Even for $p=\infty$, the Bottleneck distance is also very costly in
time to calculate, especially for $\mathbb{H}_{0}$. This explain why at some
point, we did not provide the comparison on $\mathbb{H}_{0}$. We can infer
from this analysis that topological persistent homologies do change
dramatically at seizure, a finding already obtained in previous analyses, see
Kwessi and Edwards (2021). This suggests that looking at changes in homology
landscapes could be a predictor of seizure. Given that some EEG epilepsy data
are known to contain some deterministic chaos, it be might worthwhile to study
whether persistent homology can also be used for better understanding of
chaotic data in dynamical systems.

## References

  * Belkin and Niyogi (2002) M. Belkin and P. Niyogi.  Laplacian eigenmaps and spectral techniques for embedding and clustering.  In T. G. Dietterich, S. Becker, and Z. Ghahramani, editors, _Advances in Neural Information Processing Systems 14_ , pages 585–591. Cambridge, MA: MIT Press, 2002.
  * Berry et al. (2020) E. Berry, Y.-C. Chen, J. Cisewski-Kehe, and B.T. Fasy.  Functional summaries of persistence diagrams.  _Journal of Applied and Computational Topology_ , 4:211–262, 2020.
  * Bhattacharya et al. (2015) S. Bhattacharya, R. Ghrist, and V. Kumar.  Persistent homology for path planning in uncertain environments.  _IEEETrans Robot_ , 31:578–590, 2015.
  * Chan et al. (2013) J. M. Chan, G. Carlsson, and R. Rabadan.  Topology of viral evolution.  _PNAS_ , 110(46):18566–18571, 2013.
  * Chung et al. (2009) M. K. Chung, P. Bubenik, and P. T. Kim.  Persistence diagrams of cortical surface data.  In J. L. Prince, D. L. Pham, and K. J. Myers, editors, _Information processing in medical imaging. Lecture notes in computer science_ , volume 5636 of _Springer, Berlin_ , pages 386–397, 2009.
  * Chung et al. (2023) M. K. Chung, C. G. Ramos, J. Paiva, F. B. Mathis, V. Prabharakaren, V. A. Nair, E. Meyerand, B. P. Hermann, J. R. Binder, and A. F. Struck.  Unified topological inference for brainnetworks in temporal lobe epilepsy using thewasserstein distance.  _ArXiv_ , 2023.
  * de Silva and Ghrist (2007) V. G. de Silva and R. Ghrist.  Coverage in sensor networks via persistent homology.  _Algebraic Geom Topol_ , 7:339–358, 2007.  doi: DOI10.1140/epjds/s13688-017-0109-5.
  * Emmett et al. (2016) K. Emmett, N. Schweinhart, and R. Rabadán.  Multiscale topology of chromatin folding.  In _Proceedings of the 9th EAIinternational conference on bio-inspired information and communications technologies_ , BICT’15. ICST, pages 177–180, 2016.
  * Gameiro et al. (2015) M. Gameiro, Y. Hiraoka, S. Izumi, Kramár M., K. Mischaikow, and K. Nanda.  A topological measurement of proteincompressibility.  _Jpn J Ind Appl Math_ , 32:1–1–7, 2015.
  * Ghrist (2008) R. Ghrist.  Barcodes: The persistent topology of data.  _Bull. Amer. Math. Soc._ , 45:61–75, 2008.
  * Giusti et al. (2016) C. Giusti, R. Ghrist, and D. Bassett.  Two’s company and three (or more) is a simplex.  _J Comput Neurosci_ , 41:1–14, 2016.
  * Guillemard et al. (2013) M. Guillemard, H. Boche, G. Kutyniok, and F. Philipp.  Persistence diagrams of cortical surface data.  In _10th international conference on sampling theory and applications._ , pages 309–312, 2013.
  * Hinton and Roweis (2002) G. E. Hinton and S. Roweis.  Stochastic neighbor embedding.  In S. Becker, S. Thrun, and K. Obermayer, editors, _Advances in Neural Information Processing Systems_ , volume 15. MIT Press, 2002.
  * Hotelling (1933) H. Hotelling.  Analysis of a complex of statistical variables into principal components.  _Journal of Educationbal Psychology_ , 24(417–441):498–520, 1933.
  * Hyvärinen (1999) A. Hyvärinen.  Fast and robust fixed-point algorithms for independent component analysis.  _IEEE Transactions on Neural Networks_ , 13(4–5):411–430, 1999.
  * Kwessi and Edwards (2021) E. Kwessi and L. Edwards.  Analysis of eeg time series data using complex structurization.  _Neural Computations_ , 33(7), 2021.
  * Leibon et al. (2008) G. Leibon, S. Pauls, D. Rockmore, and R. Savell.  Topological structures in the equities market network.  _Proc Natl AcadSci USA_ , 105:20589–20594, 2008.
  * Ma and Fu (2012) Y. Ma and Y. Fu.  _Manifold Learning: Theory and Applications_.  CRC Press, 2012.
  * Maletić et al. (2016) S. Maletić, Y. Zhao, and M. Rajković.  Persistent topological features of dynamical systems.  _Chaos_ , 26(5):053105, 2016.  doi: 10.1063/1.4949472.PMID:27249945.
  * Mileyko et al. (2011) Y. Mileyko, S. Mukherjee, , and J. Harer.  Probability measures on the space of persistence diagrams.  _Inverse Problems_ , 27:124007, 2011.
  * Naizait et al. (2020) G. Naizait, A. Zhitnikov, and L.-H. Lim.  Topology of deep neural networks.  _Journal of Machine Learning Research_ , 21:1–40, 2020\.
  * Otter et al. (2017) N. Otter, M. A. Porter, U. Tillman, P. Grindrod, and H. A. Harrington.  A roadmap for the computationof persistent homology.  _EPJ Data Science_ , 6(7), 2017.  doi: DOI10.1140/epjds/s13688-017-0109-5.
  * Pokorny et al. (2016) F.T. Pokorny, M. Hawasly, and S. Ramamoorthy.  Topological trajectory classification with filtrations of simplicialcomplexes and persistent homology.  _Int J Robot Res_ , 35:204–223, 2016.
  * Rizvi et al. (2017) A. Rizvi, P. Camara, E. Kandror, T. Roberts, I. Schieren, T. Maniatis, and R. Rabadán.  Single-cell topological rna-seqanalysis reveals insights into cellular differentiation and development.  _Nat Biotechnol_ , 35:551–560, 2017.
  * Sizemore et al. (2019) A. E. Sizemore, J. E. Phillips-Cremins, R. Ghrist, and D. S. Bassett.  The importance of the whole: Topological data analysis for the network neuroscientist.  _Network Neuroscience_ , 3:656––673, 2019.
  * Takens (1981) F. Takens.  Detecting strange attractors in turbulence dynamical systems and turbulence (lecture notes in mathematics), vol. 898, 1981.
  * Taylor et al. (2015) D. Taylor, F. Klimm, H. A. Harrington, M. Kramár, K. Mischaikow, M. A. Porter, and P. J. Mucha.  Topological data analysis ofcontagion maps for examining spreading processes on networks.  _Nat Commun_ , 6, 2015.  Article ID 7723.
  * Tenenbaum et al. (2000) J.B. Tenenbaum, V. de Siva, and J. C. Langford.  A global geometric frameworkfor nonlinear dimensionality reduction.  _Science_ , 290:2319–2323, 2000.
  * Torgerson (1952) W. S. Torgerson.  Multidimensional scaling: I. theory and method.  _Psychometrika_ , 17:410–419, 1952.
  * van der Maaten and Hinton (2008) L. J. P. van der Maaten and G. E. Hinton.  Visualizing data using t-sne.  _Journal of Machine Learning Research_ , 9:2579–2605, 2008\.
  * Vasudevan et al. (2013) R. Vasudevan, A. Ames, and R. Bajcsy.  Persistent homology for automatic determination of human-data based costof bipedal walking.  _Nonlinear Anal Hybrid Syst_ , 7:101–115, 2013.
  * Whitney (1936) H. Whitney.  Differentiable manifolds.  _Ann. Math_ , 37(3):645–680, 1936.
  * Xia and Wei (2014) K. Xia and G.-W. Wei.  Persistent homology analysis of protein structure, flexibility, and folding.  _Int J Numer MethodsBiomed Eng_ , 30:814–844, 2014.


# A localized kinematic structure detected in atomic carbon emission spatially
coincident with a proposed protoplanet in the HD 163296 disk

Felipe Alarcón  Department of Astronomy, University of Michigan, 323 West
Hall, 1085 S University, Ave., Ann Arbor, MI 48109, USA  Edwin A. Bergin
Department of Astronomy, University of Michigan, 323 West Hall, 1085 S
University, Ave., Ann Arbor, MI 48109, USA  Richard Teague Department of
Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Center for Astrophysics — Harvard &
Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

(Accepted November 29, 2022)

###### Abstract

Over the last five years, studies of the kinematics in protoplanetary disks
have led to the discovery of new protoplanet candidates and several structures
linked to possible planet-disk interactions. We detect a localized kinematic
bipolar structure in the HD 163296 disk present inside the deepest dust gap at
48 au from atomic carbon line emission. HD 163296’s stellar jet and molecular
winds have been described in detail in the literature; however, the kinematic
anomaly in C I emission is not associated with either of them. Further, the
velocity of the kinematic structure points indicates a component fast enough
to differentiate it from the Keplerian profile of the disk; and its atomic
nature hints at a localized UV source strong enough to dissociate CO and
launch a C I outflow, or a strong polar flow from the upper layers of the
disk. By discarding the stellar jet and previously observed molecular winds,
we explore different sources for this kinematic feature in C I emission that
could be associated with a protoplanet inflow/outflow, or disk winds.

protoplanetary disks — planet–disk interactions — radio lines: planetary
systems

††journal: ApJL††facilities: ALMA††software: astropy, (Astropy Collaboration
et al., 2013, 2018), CASA (McMullin et al., 2007), bettermoments (Teague &
Foreman-Mackey, 2018), eddy (Teague, 2019), disksurf (Teague et al., 2021),
GoFish (Teague, 2019)

##  1 Introduction

Over the past few years, the analysis of gas kinematic features traced by the
emission of millimetric rotational lines has been proposed as a useful method
to detect putative forming-planets in circumstellar disks (Pinte et al.,
2022). Embedded planets in protoplanetary disk generate kinematic footprints
in the Keplerian pattern traced by molecular emission (Perez et al., 2015).
Such footprints, either in pressure profiles or velocity deviations, have been
observed in a handful of cases, thanks to the Atacama Large Millimeter and
sub-millimeter Array (ALMA), proving the usefulness of the method (Teague et
al., 2018; Pinte et al., 2018, 2022; Izquierdo et al., 2022a).

Much of the focus on disk kinematics have been related to the observation of
planetary wakes. However, there are other possible kinematic features due to
planets that could be observed. This includes meridional flows, and the inflow
of material towards the accreting protoplanet (Fung & Chiang, 2016; Rabago &
Zhu, 2021). Meridional flows and surface infall has been isolated in HD 163296
and HD 169142 (Teague et al., 2019; Yu et al., 2021), but other protoplanet-
induced features remain undetected at present. There is the possibility of
molecular outflows and/or jets powered by the circumplanetary accretion disk
and its internal angular momentum transport (Pudritz et al., 2007a). Such
outflows linked to accreting protoplanets have been proposed in the
literature, although without observed detections so far (Wolk & Beck, 1990;
Quillen & Trilling, 1998; Fendt, 2003; Lubow & Martin, 2012; Gressel et al.,
2013).

HD 163296 disk stands out as one of the most comprehensively studied disks to
date, having been extensively observed and modeled in detail. It has a rich
dust millimeter continuum structure showing gaps, rings and a dust crescent
that has been associated with an unseen giant planet (Isella et al., 2018).
Several models have proposed that massive planets are responsible for carving
such structures (Liu et al., 2018; Zhang et al., 2018). Moreover, studies of
its gas kinematics structure have been able to infer the presence of planets
through Keplerian deviations (Teague et al., 2018, 2021; Pinte et al., 2018;
Izquierdo et al., 2022b; Calcino et al., 2022). Thus, the HD 163296 disk is
one of the best candidates to study planet formation in action.

In this paper, we present the detection of a strong kinematic deviation in the
HD 163296 disk traced by C I line emission which has not been observed in
other molecular tracers (Thi et al., 2004; Law et al., 2021; Teague et al.,
2021; Guzmán et al., 2021). We argue that such a kinematic feature can be
tentatively associated with a protoplanet candidate or a disk wind that has
not been observed before. We organize the present paper as follows. In Section
2 we describe the datasets used from the ALMA Archive and the methodology used
for the gas kinematics analysis. We present the results of our data analysis
in Section 3. In Section 4 we discuss the possible sources of the kinematic
structure observed in C I emission. Finally, we summarize the findings of this
paper in Sect. 5

Figure 1: Velocity centroid map of the C I 3P1-3P0 line emission using
bettermoments with a Gaussian fit. Overlaid are the contours of the 13CO J=2-1
centroid map from the MAPS data (Law et al., 2021) at 0$\farcs$15 resolution.
We show the region of the deepest gap in HD 163296 with dashed lines and the
minor axis of the disk with a black line. The figure shows that the C I line
emission has a strong redshifted component, $\sim$ 2 km s-1 velocity offset,
located at $\sim$0$\farcs$5 to the SW of the star, while the 13CO velocity
centroid contours show a Keplerian profile. It is noteworthy that at that
radius the deepest dust continuum gap of the HD 163296 disk is located (Isella
et al., 2018).  Figure 2: Integrated spectra of the CO J=2-1 and the C I
3P1-3P0 line emission integrated over a 60 deg wedge, 60 au long, centered at
48 au on the SW orientation of the minor axis of the disk. The dim dashed
lines show the respective noise levels of each integrated spectrum. Given that
the integrated spectrum is centered on the minor axis, a Keplerian disk should
produce a symmetric structure around the stellar velocity, such as the CO
J=2-1 emission; while the C I line emission has an important redshifted
component located in the deepest dust gap of HD 163296 at 48 au. We mirrored
the blueshifted C I emission with a dotted line to emphasize the redshifted
emission excess in the wedge.

.

Figure 3: Kinematic deviations between the velocity centroid map of the C I
3P1-3P0 line emission and the CO J=2-1 line from the MAPS survey (Öberg et
al., 2021) at matching spectral and spatial resolutions. The background is
showing the DSHARP 1.3 mm dust continuum emission of the HD 163296 disk
clipped at 65 $\mu$Jy beam-1 (Andrews et al., 2018). The CO J=2-1 velocity
centroid was used as a proxy of the Keplerian rotational profile of the disk
to isolate the kinematic anomaly observed in CI. The Figure shows a bipolar
flow centered at the 48 au gap and at the same azimuthal location than the
putative planet proposed by Isella et al. (2018).

##  2 Observations and Data Reduction

We used archival ALMA observations of the HD 163296 disk from the
2013.1.00527.S and 2015.1.01137.S projects. The 2013.1.00527.S project
observed HD 163296 for 13 minutes with 37 antennas, while the 2015.1.01137.S
project targeted the HD 163296 disk with an integration time of 4 minutes. We
combined the datasets before imaging the C I 3P1-3P0 line observation in Band
8 with a rest frequency 492.161 GHz using the CASA software (McMullin et al.,
2007). To properly combine them, both datasets were phase-centered in the
celestial coordinates J2000 17°56′21′′ -21°57′22.545′′. The dataset was self-
calibrated in two rounds with a solution interval of 30 s. The continuum was
subtracted in the UV space with the uvcontsub routine and the images were
reconstructed using the tclean algorithm in the CASA software. Different
Briggs weighting values were explored during the imaging process to find an
optimal trade-off between sensitivity and spatial resolution. The final chosen
weighting value was equal to 0.5. The channel spacing of the final cube was
600 m s-1 to match the minimum independent spectral resolution of the initial
datasets. The final RMS noise value was 33 mJy beam-1 for a beam size of
0$\farcs$65$\times$0$\farcs$40 with a position angle of 79$\fdg$8 for each
channel. Using those beam dimensions and signal-to-noise ratio, we obtained an
astrometry precision $\sim$ 0$\farcs$19, equivalent to 20 au in the HD 163296
disk.

After the datasets were combined, we analyzed the kinematics features in the C
I line emission using the bettermoments module (Teague & Foreman-Mackey,
2018), eddy (Teague, 2019), disksurf (Teague et al., 2021) and GoFish (Teague,
2019) python packages. By using those tools, we calculated the velocity-
centroid maps and spectra from the disk using the C I line emission. The
velocity-centroid map in Figure 1 was calculated by doing a Gaussian fitting
inside the bettermoments module. For completeness, the velocity-centroid maps
obtained from the Gaussian fit were compared to the ones from the quadratic
fitting, showing the same localized kinematic feature. Given the spatial
resolution of the observations, it was not possible to recover the emission
height of the C I line. Nevertheless, we were able to compare the C I
velocity-centroid map with the CO velocity-centroid from the MAPS data (Teague
et al., 2021) by subtracting them.

##  3 Results

Figure 1 shows the contours of the velocity-centroid maps of the 13CO J=2-1
line emission overlaid on the velocity-centroid map of the C I line emission.
In the Figure, the CO contours show the systemic rotational profile in the
disk, while the C I velocity-centroid map has a clear redshifted lobe located
$\sim$ 50 au to the SW of the star. The redshifted component has a velocity
offset $>$ 2 km s-1, which suggests that is not linked to the Keplerian
rotation of the disk. The redshifted lobe could be either a fast radially
outwards component or vertical motion in the disk.

We provide the C I and CO J=2-1 integrated spectrum of a 60° wedge with a 60
au radial extension centered in the gap at 48 au on the SW side of the minor
axis of the disk in Fig. 2. Given that the wedge is centered on the minor
axis, a Keplerian disk should be symmetric around the systemic velocity of the
star, 5.76 km s-1 (Teague et al., 2021). As expected, the CO spectrum is
symmetric around the stellar systemic velocity, while the C I emission has a
bright redshifted component when compared with the blueshifted Keplerian
component. This asymmetry hints at an extra velocity component in the disk
which is reflected in the velocity-centroid maps shown in Fig. 1.

We show the subtraction of the C I velocity-centroid map with the one from the
CO J=2-1 MAPS data in Figure 3 imaged at the same spatial and spectral
resolution. By doing that we isolate the gas motion that is not associated
with molecular tracers or the Keplerian disk. In the Figure, we also overlay
the 1.3 millimeter dust continuum emission from the DSHARP survey (Andrews et
al., 2018; Isella et al., 2018) to show the location of the bipolar blob when
compared with the dust continuum emission structure. In addition, the cross
represents the proposed location of a putative planet as suggested by Isella
et al. (2018). A planet at this location would be able to create the dust
continuum crescent observed within the gap at 48 au. This putative planet
appears to be spatially coincident with the C I kinematic anomaly, both
radially and azimuthally. There is an extra strong velocity deviation in the
region close to the star. In that case, however, it is harder to isolate the
emission from the inner disk given the beam size.

Finally, we note that this kinematic feature is not detected in other species
as this is not found in the emission maps of any other species with resolved
data. For this purpose we searched the emission images of CO, C2H, CN, HCO+
and numerous other species/transitions observed in the MAPS program (Law et
al., 2021; Guzmán et al., 2021) for evidence of a similar anomaly. No evidence
is present in these deep datasets. In additional support of its localized
nature, no feature from HD 163296 at that velocity is observed in large-scale
emission spectra from organic molecules with IRAM (Thi et al., 2004).

##  4 Discussion

Given its brightness, HD 163296 is one of the most well-characterized
protoplanetary disks to date. This system is associated with a well-
characterized stellar jet that is also aligned along the NE-SW minor axis,
coinciding with the localized kinematic anomaly (Xie et al., 2021). However,
the velocities of the stellar jet have values ranging over 100 km s-1; in
contrast, the C I feature is spread over only a couple km s-1. More
importantly, the jet observed to the SW of the disk is blueshifted (Günther et
al., 2013; Ellerbroek et al., 2014; Xie et al., 2021; Kirwan et al., 2022),
which is expected given the orientation of the HD 163296 disk on the sky,
i.e., the NE side being the closest to us (Teague et al., 2021). On the
contrary, the C I feature at we observe at $\sim$50 au toward the SW is
redshifted. At this location, such a feature could be interpreted as a
vertical inflow/outflow, or a radial flow going outwards. A connection between
the C I kinematic feature and the slower molecular winds observed by Klaassen
et al. (2013) and Booth et al. (2021) is also unlikely as they are also
blueshifted, with a velocity offset from the star $\sim$ -18 km s-1. At the
location closest to the disk in the SW direction, the wind observed in CO
emission is blue-shifted to $\sim$ -13 km s-1; excess C I emission is at
$\sim$8-10 km s-1. Thus, the source of the C I feature is kinematically
decoupled from either the stellar jet or the observed CO molecular wind given
its velocity and its localized nature. Another potential origin for this
feature could be the detection of a planetary wake. However, the kinematic
deviations associated with wakes are of order a few hundreds of m s-1 for
planets less massive than 1 Jupiter mass, and the Doppler flip would be
azimuthally inverted to the one we observed (Pérez et al., 2018), i.e.,
negative and positive velocity deviations are azimuthally switched given the
disk rotation.

The atomic nature of the kinematic feature requires a source of UV photons,
strong enough to photodissociate CO. If the kinematic feature is tracing C I
flows from the upper layers of the disk, it requires a physical source to
launch outflows at that speed that far from the star, which can be linked to
planets. Conversely, if the feature traces C I emission in the midplane, the
UV photons need to be produced locally deep in the disk since CO self-shields
UV photons coming from the star and the Interstellar Medium (see Appendix A).
Such a source can be linked to a circumplanetary disk actively accreting in
the midplane. We are not able to discard a stellar wind for the blue lobe, as
it is closer to the star and also blueshifted; however, if the blueshifted
lobe traces the disk wind, that still would not explain the source of the
localized redshifted lobe at $\sim$ 50 au.

Given that the C I velocity deviation is $>$ 2 km s-1, i.e. $\delta
v_{0}/v_{K}\sim 0.33$, for the anomaly to be associated with a planet, it has
to be around 10 MJup for $\alpha\sim 10^{-3}$ (see Figure 7 in Rabago & Zhu
(2021)). This estimate can be compared to previously reported values in the
literature using observations of the disk at different wavelengths, e.g.,
Zhang et al. (2018) provide an estimate of a 2.18 MJup planet with
$\alpha=10^{-3}$ using dust millimeter continuum. Further, Mesa et al. (2019)
put an upper limit of 7 MJup for any possible companion with direct imaging
observations.

###  4.1 Potential Solution #1: Protoplanet Infall

Vertical infall of material is expected in the presence of a massive planet,
reaching levels even comparable to the local sound speed (Fung & Chiang,
2016). Tracing the infall through C I line emission might be possible if it
were dominated by disk surface gas directly exposed to stellar ultraviolet
(UV) radiation, as the photodissociation of molecules becomes important. In
the case of an embedded massive planet, it could carve a gap deep enough so
the local pressure gradient dominates the gas flow even in the presence of
disk winds (Weber et al., 2022). Moreover, the depleted gap will increase the
penetration of the stellar UV flux into the disk (Facchini et al., 2018;
Alarcón et al., 2020). In that case, we expect the C I abundance to increase
in layers closer to the midplane, such that neutral atomic carbon might be a
tracer of protoplanet infall. Figure 4 shows the C I abundance in a
thermochemical model of HD 163296 including its dust substructure. The Figure
shows that C I remains abundant at lower layers in the disk, closer to the
midplane within the gap near 50 au. Moreover, the contours showing the CO
abundance illustrate that as that region is more exposed to UV radiation, the
photodissociation rate increases. If we assume the feature is tracing the
surface of the disk, we can also put further constraints on the putative
planet’s mass by using the local sound speed in the disk. The sound speed in
our thermochemical model at the C I layer is $c_{s}=1.5$ km s-1, so the
meridional speed of the C I inflow agrees with the models of Fung & Chiang
(2016) for 1 to 4 Jupiter mass planets, i.e., $v_{z}\sim 1-2c_{s}$, at the
location of the planet.

Our model shows that C I is clearly present in the gap, but the effects of
photodissociation via stellar photons do not destroy enough CO in the midplane
to raise the C I abundance to values significant enough to be observable. From
thermochemical models of HD 163296 with rac2d (Du & Bergin, 2014a), the
predicted C I emission for the 3P1-3P0 line matching the beam size of the
observations is 40 mJy beam-1, which is not detectable given the noise value.

The peak intensity at the kinematic anomaly allows us to provide estimates of
the vertical gas flow. By using the redshifted velocity, $v\sim$ 2 km s${-1}$
and the disk inclination, $i=$ 46.7 deg (Öberg et al., 2021), we can estimate
a lower limit for the infall velocity, $v_{z}=v\cos{i}=1.44$ km s-1. Then, we
can estimate the infall timescale using $t_{\mathrm{in}}=h/v_{z}$, where $h$,
is the emission height of C I. From Figure 4, we use $h\approx 10$ au as a
tentative C I emission height, where most of the C I is concentrated. Then,
the estimated infall timescale is $t_{\mathrm{in}}=33$ yr. With the inferred C
I column density and the area of the redshifted anomaly we can estimate the
mass accretion rate by using the following relationship:

$\dot{M}=m_{\rm H}A\frac{N_{CI}}{X_{CI}}\frac{1}{t_{\mathrm{in}}},$ (1)

where $A$ is the area of the redshifted anomaly, $m_{H}$ the hydrogen mass,
$N_{CI}$ the inferred C I column density, $X_{CI}$ the C I abundance, and
$t_{\mathrm{in}}$ the infall timescale. Using the peak intensity at the
velocity anomaly, $I_{\nu}=0.09$ Jy beam-1, the inferred C I column density is
$N_{CI}=$ 3.5 x 1017 cm-2 for an excitation temperature Texc = 45 K (see
Appendix B). Taking into account a transverse area, $A=\pi 20^{2}$ au2, and
$X_{CI}=10^{-5}$, the accretion rate is: $\dot{M}=2\times 10^{-4}\ M_{\rm
Jup}\ \mathrm{yr}^{-1}$, which is within the expected range of accretion rates
for growing planets from the literature (Kley et al., 2001; Szulágyi et al.,
2014; Tanigawa & Tanaka, 2016).

Figure 4: C I abundance in a thermochemical model of HD 163296 including the
gaps and ring from dust continuum emission (Isella et al., 2018) ran with
rac2d (Du & Bergin, 2014b). Contours are drawn for the CO abundance in the
disk. The Figure shows that inside the gap at 48 au, particularly towards the
outer edge, the C I abundance increases towards the midplane due to the
transparency to UV in that region. Thus, it is possible that C I line emission
can trace vertical protoplanet infall in the disk.

An additional point in favor of this scenario (and Solution #2) is that the
dust crescent is located in what would be the L5 Lagrangian point of the
planet. It is expected that a massive planet, such as the one proposed to
carve the dust gap at 48 au, can cause a dust pile-up in the L4 and L5
Lagrangian points (Montesinos et al., 2020), which has also been observed in
the LkCa 15 disk (Long et al., 2022). Therefore, the azimuthal and radial
location of the kinematic C I feature agrees with the presence of multiple
substructures in the disk.

###  4.2 Potential Solution #2: Protoplanet Outflow or Radial De-cretion

If the C I emission is arising from within deeper layers in the disk, there
might be a local UV source generated by the protoplanetary accretion. This
would answer one puzzling aspect of this anomaly - why it is not seen in
molecular emission (see § 3). If a UV source is present to dissociate CO, then
other species, which are more easily photodestroyed, would be as well. In the
presence of optically thin C18O emission, it is possible that a localized
emission gap caused by CO photodissociation would be expected, although there
will be an interplay between temperature and optical depth. However, the C18O
J=2-1 emission is still optically thick at 50 au in the HD 163296 disk, so it
does not trace the midplane at 50 au.

If the C I kinematic feature is not associated with infall onto a protoplanet,
then it could be a localized source of outflow from an accreting protoplanet,
and the kinematic feature is tracing an outflow behind the disk. It is clear
that protoplanets are born with accretion disks and emission from H$\alpha$
has been detected towards PDS 70 (Haffert et al., 2019), which is associated
with accretion onto the young protoplanets. Recently, a faint point source
emission of H$\alpha$ has been observed in the HD 163296 disk (Huélamo et al.,
2022) at the minor axis in the SW of the disk, although at a closer radius,
$r=$ 171 mas. Given the presence of accretion disks, subsequent outflows or
possible jets are expected in order to produce the loss of angular momentum
(Pudritz et al., 2007b; Pascucci et al., 2022). Given that the planet-disk
interactions are strong enough, the C I emission feature can also be tracing
the radial de-cretion flows in the midplane towards the gap’s edges (Fung &
Chiang, 2016; Rabago & Zhu, 2021; Calcino et al., 2022).

Here the gas flow rate from Section 4.1 would reflect the outflow mass loss
rate. Since we are tracing the backside of the disk, the diameter of the
emitting area is constrained by the width of the dust gap, which is 20 au
(Huang et al., 2018), so we assume the area of the anomaly to be, $A=\pi$102
au2. This emitting area is a factor of four smaller than the one assumed in
our Potential Solution #1, so the inferred accretion rate is also a factor of
four smaller. We can then use the relationship $\dot{M}_{\rm outflow}\approx
0.1\dot{M}_{\rm acc}$ (Pudritz et al., 2007b), leading to an estimated
accretion rate of $\dot{M}=2.5\times 10^{-4}$ $M_{\rm Jup}\ \mathrm{yr}^{-1}$.
Even though this value is still within the possible solutions for planetary
accretion rates, it would be for a very strong accretor. This solution would
also require the presence of a localized UV field to account for the lack of
detection of a molecular counterpart.

###  4.3 Potential Solution #3: Disk Wind

Another possibility is that the C I feature is associated with disk winds.
Disk winds are associated with the jets and magnetic activity of the central
star or radiation-driven in the outer regions of disks. In detriment of this
scenario, the stellar jet is blueshifted to the SW (Xie et al., 2021), which
is the closest part of the disk. For the observed molecular wind in the disk,
the velocity offset being at least a few tens of km s-1 (Booth et al., 2021).
Additionally, theoretical and observational evidence have shown that disk
winds associated with magnetic fields are launched from the inner disk, $<$ 10
au for Class II disks (Goodson et al., 1999; Matt & Pudritz, 2005; Pascucci et
al., 2022). There is a possibility that the feature is tracing an FUV-driven
wind on the far side of the disk, although it is expected that this would be
observed as a diffused extended emission as well. Additionally, the bipolar
nature of the feature offset from the central star makes this scenario less
likely. Nevertheless, the blueshifted lobe is close to the stellar location,
without additional high-resolution data we cannot discard that the blueshifted
lobe is tracing a stellar wind. However, this scenario does not provide an
origin for the offset of the bipolar structure and the localized redshifted
lobe at $\sim$50 au.

Further, observations with the APEX telescope did not show strong C I emission
in HD 163296 (Kama et al., 2016). This points out that the observed kinematic
feature is local rather than a large-scale emission, which is usually the case
for disk winds. Nevertheless, molecular winds launched from the disk surface
at a few km s-1 are still possible.

##  5 Summary

In this work, we have presented the detection of a strong kinematic deviation
in C I line emission in the HD 163296 protoplanetary disk. We proposed that
such structure is linked to either the presence of a strong UV source in the
midplane of the disk, dissociating surrounding CO, which would explain the
strong feature in C I ; or the feature is coming from the disk atmosphere and
is tracing inflow of gas into the midplane. The fact that is not observed in
CO or other molecular tracers, such as the ones observed by the MAPS program,
suggests that this kinematic component is mostly atomic in nature and not
molecular. Its atomic nature constrains the expected local physical and
chemical conditions at the source of this feature. The evidence for this
feature is given below:

  1. 1.

It is located in the strongly depleted inner dust gap, which has been proposed
to be carved by a massive Jovian planet.

  2. 2.

The projection of kinematic structure is parallel to the disk minor axis.

  3. 3.

The structure is observed in two independent datasets (see Appendix C).

  4. 4.

The kinematic structure appears to be atomic in composition, requiring a
localized UV source.

  5. 5.

The dust continuum crescent at the inner dust gap in HD 163296 is located in
one of the dust pile-up Lagrangian points of the putative planet.

  6. 6.

Planetary wakes from hydrodynamic simulations are at the level of a few
hundreds of m s-1. Both the molecular wind and the stellar jet’s component to
the SW of the disk are blueshifted, while the strongest emission in the C I
kinematic structure is redshifted. The velocity of the molecular wind in HD
163296 ranges between -20 and -18 km s-1 offset from the stellar velocity
(Booth et al., 2021), while the blue-shifted component of the stellar jet has
a velocity, $v_{\rm jet}=280$ km s-1 (Xie et al., 2021). Therefore, we discard
other sources of velocity structure.

Deeper observations of C I 3P1-3P0 line emission at higher spatial resolution
and precision are needed to confirm the link between the C I kinematic feature
in the HD 163296 disk and a putative protoplanet. If a protoplanet nature is
confirmed, the C I kinematic feature would be the first detection of an atomic
carbon protoplanetary footprint, introducing a novel technique of protoplanet
detection, while providing further constraints in the planet formation process
and composition of the planet-feeding gas. If the source of the C I feature is
coming from either a disk wind or planetary accretion, both of which require a
local UV source, it would produce a distinctive feature in atomic forbidden
lines, some of which are potentially detectable with JWST. Additionally, a
more robust velocity-centroid map of the HD 163296 disk in C I emission with
ALMA would provide further evidence to disentangle the different possible
scenarios.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00527.S.
, ADS/JAO.ALMA#2015.1.01137.S. ALMA is a partnership of ESO (representing its
member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST
and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and
NAOJ. The National Radio Astronomy Observatory is a facility of the National
Science Foundation operated under cooperative agreement by Associated
Universities, Inc. The National Radio Astronomy Observatory is a facility of
the National Science Foundation operated under cooperative agreement by
Associated Universities, Inc. We thank the anonymous referee for the
insightful reading of the manuscript, and for providing useful and helpful
comments to improve the quality of this work.

##  Appendix A CO Self-Shielding

One of the proposed solutions for the C I kinematic feature is due to
protoplanet infall from the upper layers of the disk. However, since the
feature is not observed in molecular tracers, it requires that the inflow is
actually coming from the very top layers of the disk, or the atomic layer of
the disk is pushed down due to dust depletion. Another important factor to
consider, in particular for C I is CO self-shielding (Visser et al., 2009). As
CO self-shields itself it avoids its own destruction due to high-energy
photons. We show the CO self-shielding factor in Figure 5. It shows that in a
99% dust-depleted gap, the CO self-shielding layer is pushed down. However,
that push is not significant enough to put the C I layer inside the Hill’s
sphere of a 3 Jupiter mass planet at that location. if the inflow scenario is
the source of the kinematic asymmetry, it points to a chemical mixing between
the upper layers of the disk and planet-feeding gas. If the source of the C I
emission is in the upper layers, it points to an additional UV source in the
midplane. This requires a powerful source of high-energy photons, which could
be an accretion disk around a protoplanet.

Figure 5: CO self-shielding for the HD 163296 disk model with a dust-depleted
gap at 48 au. CO self-shields itself even inside the gap. However, the content
depletion inside the gap allows C I formation through CO destruction at the
deeper layer. The model shows that CO self-shields $\sim$ 99% of the CO
photodissociation photons 10 au above the midplane at 48 au.

##  Appendix B Calculation of the C I Column Density

We use the peak intensity to calculate the C I column density in the kinematic
anomaly, in order to estimate the mass flow. From radiative transfer equations
for continuum subtracted line emission, the observed intensity is given by the
following:

$J_{\nu_{0},\mathrm{obs}}=J_{\nu_{0}}(T_{\rm
exc})\Big{(}1-\mathrm{exp}(-\tau_{\nu_{\rm peak}})\Big{)},$ (B1)

where $J_{\nu_{0}}$ is the mean intensity at the line peak, $T_{\rm exc}$ the
excitation temperature, which we set to 45 K, and $\tau_{\nu_{\rm peak}}$ the
optical depth at the line peak.

The optical depth at the line peak is given by:

$\tau_{\nu_{\rm
peak}}=\frac{c^{2}}{8\pi\nu^{2}}\Big{(}\mathrm{exp}\Big{(}\frac{h\nu}{k_{B}T}\Big{)}-1\Big{)}A_{ul}\phi_{\nu_{\rm
peak}}N_{u},$ (B2)

with $c$ the speed of light, $h$ the Planck’s constant, $k_{B}$ the
Boltzmann’s constant, $A_{ul}$ the Einstein coefficient for the respective
line, $\phi_{\nu_{\rm peak}}$ the value of the line profile at the peak of the
line, and $N_{u}$ the column density of the targeted species at the excited
state.

The value of the line profile at the peak is given by:

$\phi_{\nu_{\rm peak}}=\frac{1}{\sqrt{2\pi}\sigma},$ (B3)

with:

$\sigma^{2}=\frac{c^{2}}{2\nu^{2}\Big{(}\frac{2k_{B}T}{m_{CI}}+v_{\rm
turb}^{2}\Big{)}},$ (B4)

where $m_{CI}$ is the mass of C I , $v_{\rm turb}$ a turbulent velocity
component, which we set to 200 m s-1, and $T$ the temperature.

Assuming $J_{\nu}=4\pi B_{\nu}$, the blackbody emission, and combining it with
Eqs. B1, B2 and B3; the excited column density of C I is recovered. To
estimate the total C I column density, we used the partition function as
follows:

$N_{u}=N_{CI}\frac{g_{u}}{Q}\exp{\Big{(}\frac{-E_{u}}{k_{B}T_{\rm
exc}}\Big{)}},$ (B5)

with $g_{u}$ the degeneracy of the upper level, $Q$ the partition function,
$E_{u}$ the energy of the upper level, and $T_{\rm exc}$ the excitation
temperature. Above 30 K, the brightness of the C I 3P1-3P0 line is relatively
insensitive to temperature. Excitation temperatures between 40 K and 100 K
produce column densities, NCI, in the range 3.2-3.9 1017 cm-2, i.e, only a 22%
change. Thus, adopting a value of 45 K for the column density retrieval does
not produce significant changes in the column density given the other sources
of uncertainty. All the specific constant values for the C I 3P1-3P0 line are
taken from the Leiden Atomic and Molecular Database, LAMDA (Schöier et al.,
2005).

##  Appendix C Velocity Centroid in different Datasets

For completeness, we also show the velocity centroids in Figure 6. Both
datasets show the kinematic feature, although the statistical significance of
the feature for each dataset independently is not high enough to assure the
presence of the feature.

Figure 6: Velocity centroid map of the C I 3P1-3P0 line for the 2013.1.00527.S
and 2015.1.01137.S. datasets on the left and right respectively. The kinematic
feature in the SW of the disk is observed in both datasets independently.

## References

  * Alarcón et al. (2020) Alarcón, F., Teague, R., Zhang, K., Bergin, E. A., & Barraza-Alfaro, M. 2020, ApJ, 905, 68, doi: 10.3847/1538-4357/abc1d6
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Booth et al. (2021) Booth, A. S., Tabone, B., Ilee, J. D., et al. 2021, ApJS, 257, 16, doi: 10.3847/1538-4365/ac1ad4
  * Calcino et al. (2022) Calcino, J., Hilder, T., Price, D. J., et al. 2022, ApJ, 929, L25, doi: 10.3847/2041-8213/ac64a7
  * Du & Bergin (2014a) Du, F., & Bergin, E. A. 2014a, ApJ, 792, 2, doi: 10.1088/0004-637X/792/1/2
  * Du & Bergin (2014b) —. 2014b, ApJ, 792, 2, doi: 10.1088/0004-637X/792/1/2
  * Ellerbroek et al. (2014) Ellerbroek, L. E., Podio, L., Dougados, C., et al. 2014, A&A, 563, A87, doi: 10.1051/0004-6361/201323092
  * Facchini et al. (2018) Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104, doi: 10.1051/0004-6361/201731390
  * Fendt (2003) Fendt, C. 2003, A&A, 411, 623, doi: 10.1051/0004-6361:20034154
  * Fung & Chiang (2016) Fung, J., & Chiang, E. 2016, ApJ, 832, 105, doi: 10.3847/0004-637X/832/2/105
  * Goodson et al. (1999) Goodson, A. P., Böhm, K.-H., & Winglee, R. M. 1999, ApJ, 524, 142, doi: 10.1086/307779
  * Gressel et al. (2013) Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59, doi: 10.1088/0004-637X/779/1/59
  * Günther et al. (2013) Günther, H. M., Schneider, P. C., & Li, Z. Y. 2013, A&A, 552, A142, doi: 10.1051/0004-6361/201220772
  * Guzmán et al. (2021) Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, ApJS, 257, 6, doi: 10.3847/1538-4365/ac1440
  * Haffert et al. (2019) Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nature Astronomy, 3, 749, doi: 10.1038/s41550-019-0780-5
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42, doi: 10.3847/2041-8213/aaf740
  * Huélamo et al. (2022) Huélamo, N., Chauvin, G., Mendigutía, I., et al. 2022, arXiv e-prints, arXiv:2210.02212.  https://arxiv.org/abs/2210.02212
  * Isella et al. (2018) Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49, doi: 10.3847/2041-8213/aaf747
  * Izquierdo et al. (2022a) Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022a, ApJ, 928, 2, doi: 10.3847/1538-4357/ac474d
  * Izquierdo et al. (2022b) —. 2022b, ApJ, 928, 2, doi: 10.3847/1538-4357/ac474d
  * Kama et al. (2016) Kama, M., Bruderer, S., Carney, M., et al. 2016, A&A, 588, A108, doi: 10.1051/0004-6361/201526791
  * Kirwan et al. (2022) Kirwan, A., Murphy, A., Schneider, P. C., et al. 2022, arXiv e-prints, arXiv:2205.06569.  https://arxiv.org/abs/2205.06569
  * Klaassen et al. (2013) Klaassen, P. D., Juhasz, A., Mathews, G. S., et al. 2013, A&A, 555, A73, doi: 10.1051/0004-6361/201321129
  * Kley et al. (2001) Kley, W., D’Angelo, G., & Henning, T. 2001, ApJ, 547, 457, doi: 10.1086/318345
  * Law et al. (2021) Law, C. J., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 3, doi: 10.3847/1538-4365/ac1434
  * Liu et al. (2018) Liu, S.-F., Jin, S., Li, S., Isella, A., & Li, H. 2018, ApJ, 857, 87, doi: 10.3847/1538-4357/aab718
  * Long et al. (2022) Long, F., Andrews, S. M., Zhang, S., et al. 2022, ApJ, 937, L1, doi: 10.3847/2041-8213/ac8b10
  * Lubow & Martin (2012) Lubow, S. H., & Martin, R. G. 2012, ApJ, 749, L37, doi: 10.1088/2041-8205/749/2/L37
  * Matt & Pudritz (2005) Matt, S., & Pudritz, R. E. 2005, ApJ, 632, L135, doi: 10.1086/498066
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Mesa et al. (2019) Mesa, D., Langlois, M., Garufi, A., et al. 2019, MNRAS, 488, 37, doi: 10.1093/mnras/stz1662
  * Montesinos et al. (2020) Montesinos, M., Garrido-Deutelmoser, J., Olofsson, J., et al. 2020, A&A, 642, A224, doi: 10.1051/0004-6361/202038758
  * Öberg et al. (2021) Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1, doi: 10.3847/1538-4365/ac1432
  * Pascucci et al. (2022) Pascucci, I., Cabrit, S., Edwards, S., et al. 2022, arXiv e-prints, arXiv:2203.10068.  https://arxiv.org/abs/2203.10068
  * Pérez et al. (2018) Pérez, S., Casassus, S., & Benítez-Llambay, P. 2018, MNRAS, 480, L12, doi: 10.1093/mnrasl/sly109
  * Perez et al. (2015) Perez, S., Dunhill, A., Casassus, S., et al. 2015, ApJ, 811, L5, doi: 10.1088/2041-8205/811/1/L5
  * Pinte et al. (2022) Pinte, C., Teague, R., Flaherty, K., et al. 2022, arXiv e-prints, arXiv:2203.09528.  https://arxiv.org/abs/2203.09528
  * Pinte et al. (2018) Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13, doi: 10.3847/2041-8213/aac6dc
  * Pudritz et al. (2007a) Pudritz, R. E., Ouyed, R., Fendt, C., & Brandenburg, A. 2007a, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 277.  https://arxiv.org/abs/astro-ph/0603592
  * Pudritz et al. (2007b) Pudritz, R. E., Ouyed, R., Fendt, C., & Brandenburg, A. 2007b, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 277.  https://arxiv.org/abs/astro-ph/0603592
  * Quillen & Trilling (1998) Quillen, A. C., & Trilling, D. E. 1998, ApJ, 508, 707, doi: 10.1086/306421
  * Rabago & Zhu (2021) Rabago, I., & Zhu, Z. 2021, MNRAS, 502, 5325, doi: 10.1093/mnras/stab447
  * Schöier et al. (2005) Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369, doi: 10.1051/0004-6361:20041729
  * Szulágyi et al. (2014) Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65, doi: 10.1088/0004-637X/782/2/65
  * Tanigawa & Tanaka (2016) Tanigawa, T., & Tanaka, H. 2016, ApJ, 823, 48, doi: 10.3847/0004-637X/823/1/48
  * Teague (2019) Teague, R. 2019, The Journal of Open Source Software, 4, 1220, doi: 10.21105/joss.01220
  * Teague (2019) Teague, R. 2019, The Journal of Open Source Software, 4, 1632, doi: 10.21105/joss.01632
  * Teague et al. (2019) Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378, doi: 10.1038/s41586-019-1642-0
  * Teague et al. (2018) Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12, doi: 10.3847/2041-8213/aac6d7
  * Teague & Foreman-Mackey (2018) Teague, R., & Foreman-Mackey, D. 2018, Research Notes of the American Astronomical Society, 2, 173, doi: 10.3847/2515-5172/aae265
  * Teague et al. (2021) Teague, R., Law, C. J., Huang, J., & Meng, F. 2021, Journal of Open Source Software, 6, 3827, doi: 10.21105/joss.03827
  * Teague et al. (2021) Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18, doi: 10.3847/1538-4365/ac1438
  * Thi et al. (2004) Thi, W. F., van Zadelhoff, G. J., & van Dishoeck, E. F. 2004, A&A, 425, 955, doi: 10.1051/0004-6361:200400026
  * Visser et al. (2009) Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323, doi: 10.1051/0004-6361/200912129
  * Weber et al. (2022) Weber, M. L., Ercolano, B., Picogna, G., & Rab, C. 2022, MNRAS, 517, 3598, doi: 10.1093/mnras/stac2954
  * Wolk & Beck (1990) Wolk, S. J., & Beck, S. C. 1990, PASP, 102, 745, doi: 10.1086/132697
  * Xie et al. (2021) Xie, C., Haffert, S. Y., de Boer, J., et al. 2021, A&A, 650, L6, doi: 10.1051/0004-6361/202140602
  * Yu et al. (2021) Yu, H., Teague, R., Bae, J., & Öberg, K. 2021, ApJ, 920, L33, doi: 10.3847/2041-8213/ac283e
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47, doi: 10.3847/2041-8213/aaf744