#  SkyCURTAINs: Model agnostic search for Stellar Streams with Gaia data

Debajyoti Sengupta,1 Stephen Mulligan,1 David Shih,2 John Andrew Raine,1 and
Tobias Golling1
1Département de physique nucléaire et corpusculaire, University of Geneva,
Switzerland
2NHETC, Dept. of Physics and Astronomy, Rutgers, Piscataway, NJ 08854, USA
Contact e-mail<EMAIL_ADDRESS>
###### Abstract

We present SkyCURTAINs, a data driven and model agnostic method to search for
stellar streams in the Milky Way galaxy using data from the _Gaia_ telescope.
SkyCURTAINs is a weakly supervised machine learning algorithm that builds a
background enriched template in the signal region by leveraging the
correlation of the source’s characterising features with their proper motion
in the sky. This allows for a more representative template of the background
in the signal region, and reduces the false positives in the search for
stellar streams. The minimal model assumptions in the SkyCURTAINs method allow
for a flexible and efficient search for various kinds of anomalies such as
streams, globular clusters, or dwarf galaxies directly from the data. We test
the performance of SkyCURTAINs on the GD-1 stream and show that it is able to
recover the stream with a purity of 75.4% which is an improvement of over 10%
over existing machine learning based methods while retaining a signal
efficiency of 37.9%.

###### keywords:

Galaxy: Stellar Content – Galaxy: Structure – Stars: Kinematics and Dynamics –
Methods: Weakly Supervised Machine Learning

††pubyear: 2024

##  1 Introduction

When smaller gravitationally bound systems such as globular clusters or
satellite dwarf galaxies, are disrupted by their host galaxy, the stars in
these systems are tidally stripped off. This results in a stream of stars,
named stellar streams, which, over time trace out the orbit of the progenitor
system. Since the interactions between these large-scale gravitationally bound
systems occur over a very long timescale, real time observations of these
events are impossible. Stellar streams are therefore an excellent alternative
probe into the merger history of these systems (Johnston, 1998; Helmi & White,
1999; Carlberg, 2017; Vera-Casanova et al., 2022; Belokurov et al., 2006).
Moreover, the orbits of these streams are sensitive to the gravitational
potential of the host galaxy, and thus can be used to constrain the mass
distribution in it (Johnston et al., 1999; Ibata et al., 2001; Koposov et al.,
2010; Sanders & Binney, 2013; Banik & Bovy, 2019). Over time, due to
gravitational interaction with the surrounding matter, the shape of these
streams change, and the density perturbations therein, such as gaps and spurs,
can also provide insights into the dark matter distribution in the galaxy
(Carlberg et al., 2012; Varghese et al., 2011; Sanders et al., 2016; Bonaca et
al., 2019, 2020) and its properties (Purcell et al., 2012; Necib et al.,
2019). The study of stellar streams is thus crucial to understanding the
formation and evolution of galaxies, and the content thereof.

The _Gaia_ mission (Gaia Collaboration et al., 2018) has provided an
unprecedented dataset of stars in the Milky Way, with accurate astrometric and
photometric measurements. This wealth of data has allowed the development of
several techniques to detect stellar streams (Malhan & Ibata, 2018; Malhan et
al., 2018; Yuan et al., 2018; Meingast, Stefan & Alves, João, 2019; Borsato et
al., 2019; Meingast, Stefan et al., 2019; Ibata et al., 2021). In general,
these methods leverage the astrophysics of stellar streams, such as their
grouping in chemical composition and kinematics, to identify the stream
candidates. For instance, the Streamfinder algorithm (Malhan & Ibata, 2018;
Malhan et al., 2018) assumes a specific model for the gravitational potential
of the Milky Way galaxy, and searches for stars occupying the same
hyperdimensional tubes through a six-dimensional positional and velocity
space.

More recently, several machine learning techniques have been employed to
detect stellar streams. Particularly, Via Machinae (Shih et al., 2021, 2023),
and CWoLa (Pettee et al., 2023) are fully data-driven and have very minimal
model assumptions about the streams. These techniques were originally
introduced in the context of High Energy Physics to find localised
overdensities in the feature space. In the case of kinematically cold stellar
streams, the member stars are expected to produce localised overdensities in
the proper motion feature. One can define a signal region (SR) based on the
proper motion, where there is an increased population of a stellar stream
stars, and side bands (SB1, SB2) on either side of the SR, where the stream
members are not expected to be present (or at a far lower rate, compared to
the SR).

Via Machinae (1.0, 2.0) based on ANODE (Nachman & Shih, 2020), consists of
conditional generative models, that learns the probability distribution of the
kinematic, photometric, and astrometric features of the stars in the SR and
SB, and constructs the likelihood ratio
$\frac{P_{\mathrm{SR}}(\mathbf{x})}{P_{\mathrm{bg}}(\mathbf{x})}$ to tag
anomalous stars. Here, $P_{\mathrm{SR}}$ is the probability distribution of
the stars in the signal region and $P_{\mathrm{bg}}$ is the conditionally
interpolated background density, and $\mathbf{x}$ is the feature vector of the
stars, over which the densities are defined. Thereafter, a line finding
algorithm is used to filter out the stream. The Via Machinae method has been
shown to be very effective at finding streams, however, is computationally
expensive, as it requires training two generative models on the data.

CWoLa, originally introduced in Metodiev et al. (2017), is a computationally
lightweight method that uses a weakly supervised learning approach to detect
streams. Given the SR and SB, CWoLa trains a classifier to distinguish between
the two regions, and then uses the classifier to tag the stars in the SR. The
classifier effectively learns the likelihood ratio
$\frac{P_{\mathrm{SR}}(\mathbf{x})}{P_{\mathrm{SB}}(\mathbf{x})}$. Here,
$P_{\mathrm{SR}}$ is the probability distribution of the stars in the signal
region and $P_{\mathrm{SB}}$ is the probability distribution of the stars in
the side bands. However, the performance of CWoLa is dependent on the choice
of features used in the training. If the selected features are correlated with
the proper motion, the classifier may be biased and produce false excesses in
the SR, even in the absence of a stream.

It is possible to circumvent this bias if a suitable template of the
background is constructed to be used in the CWoLa method. We propose
SkyCURTAINs, that constructs a background-enriched template of the stars in
the SR in a data driven manner. SkyCURTAINs is based on CurtainsF4F, a method
originally developed for anomaly detection in High Energy Physics introduced
in (Raine et al., 2023; Sengupta et al., 2023). CurtainsF4F is a data-driven
weakly supervised strategy that extends the CWoLa method to mitigate the
problem of correlation of discriminatory features with the proper motion
feature. We leverage the correlation of the features with the proper motion to
generate a template in the signal region using the sidebands. This alleviates
the need to sample data from the SB for CWoLa, and results in a template that
is more representative of the background in the SR. One can then use the CWoLa
method to tag the stars in the SR by training a classifier on the template of
the SR data, followed by a line finding algorithm to identify the stream.

Constructing a background enriched template significantly reduces false
positives, which is a big advantage of the SkyCURTAINs method over the
standalone CWoLa method. As we will see in section 3 SkyCURTAINs has a modular
design, and its data efficiency in training allows for an efficient scaling of
the method to larger number of patches.

##  2 Dataset

We demonstrate the SkyCURTAINs method on the _Gaia_ Data Release 2 (GDR2)
(Gaia Collaboration et al., 2018) dataset. GDR2 contains detailed astrometric
and photometric information for over 1.3 billion sources in the Milky Way
galaxy. The dataset characterises the source by the right ascension ($\alpha$)
and declination ($\delta$), the parallax ($\varpi$), proper motions in right
ascension ($\mu_{\alpha}$) and declination ($\mu_{\delta}$), the apparent
magnitude ($G$), and the colour information in the form of the GBP and GRP
bands ($G_{\text{BP}}-G_{\text{RP}}$). The newer _Gaia_ Data Release 3 (GDR3)
comes with improved measurements on radial velocities, but as the SkyCURTAINs
method does not utilise this information, we use the GDR2 dataset. This allows
for a direct comparison with the Via Machinae and CWoLa methods, which were
developed using the GDR2 dataset.

We use the GD-1 stream as the main stream candidate to benchmark and validate
the SkyCURTAINs method. The GD-1 stream is a long and dense stream in the
Milky Way galaxy, discovered in 2006 by Grillmair & Dionatos (2006). The
SkyCURTAINs method uses the stellar membership of the GD-1 stream in (Price-
Whelan & Bonaca, 2018), hereafter referred to as PWB18, as the ground truth to
validate the method. These studies use selections in position, proper motion,
colour, and magnitude space to identify the stars that are members of the GD-1
stream. Although these membership labels are likely not complete and can not
be considered as ground truth, they nonetheless provide a crucial reference
for the validation of the SkyCURTAINs method.

Following (Shih et al., 2021, 2023; Pettee et al., 2023), we choose to divide
the GDR2 dataset into overlapping circular patches of $15^{\circ}$ radius. We
re-center each patch and use patch local coordinates ($\phi$, $\lambda$) and
the corresponding proper-motion ($\mu_{\phi}$, $\mu_{\lambda}$), such that
each patch is centred at $(\alpha_{0},\delta_{0})$ = $(0^{\circ},0^{\circ})$.
This is done to ensure that each patch has an approximately Euclidean distance
metric.

Figure 1: Marginal distribution of the features used in the SkyCURTAINs
method. The stars identified by PWB18 study, are shown in red, and the
background stars are shown in blue. All distributions are normalised to 1, to
help visualise the qualitative differences between the features of background
and stream like stars. The stars in the plot are from the patch with
coordinates centred at $\alpha$ = $146.9^{\circ}$, and $\delta$ =
$35.6^{\circ}$

SkyCURTAINs uses six features associated with each star: kinematic
($\mu_{\phi}^{*}$ = $\mu_{\phi}\cos\lambda$, $\mu_{\lambda}$), spatial
($\phi$, $\lambda$), and photometric ($G$, $G_{\text{BP}}-G_{\text{RP}}$). The
marginal distribution of these features is shown in Figure 1. The entire patch
is used to train the CurtainsF4F model. However, additional fiducial cuts are
applied to the data for the downstream tasks. These include:

  * •

Kinematic cut: $|\mu_{\phi}^{*}|>2$ mas/yr OR $|\mu_{\lambda}|>2$ mas/yr.

  * •

Photometric cut: $0.5\leq G_{\text{BP}}-G_{\text{RP}}\leq 1$, and $G<20.2$.

The kinematic cuts are applied to reject distant stars that produce an
overdensity in the proper motion at $\sim 0$ mas/yr and reduce the sensitivity
of the model to overdensities produced by stellar streams. The cut on
magnitude removes stars that are too dim and ensures we have a uniform
coverage of stars from the _Gaia_ dataset. The cut on colour isolates older,
low-metallicity stars, which are more likely to be stellar stream members.

##  3 SkyCURTAINs Method

SkyCURTAINs is a two stage approach to find stellar streams in a model
agnostic manner. The first stage is CurtainsF4F followed by a CWoLa step. This
stage is used to infer a threshold to select the candidate signal stars for
the second stage. The first stage flags all overdensities as anomalous. But as
we are looking for stellar streams, we need to identify line-like structures
in the candidate signal stars’ population. This is done by the second stage of
the SkyCURTAINs method, which uses the Hough transform (Hough, 1962) for line
detection. Details of training and implementation of the two stages are
discussed in the following sections.

###  3.1 CurtainsF4F

CurtainsF4F constructs a background-enriched template in the signal region by
learning a conditional transformation of the features from the sidebands to
the signal region, as a function of the proper motion. CurtainsF4F uses a
maximum likelihood loss on the transported data and the target data using the
Flows for Flows method introduced in Golling et al. (2023) to learn this
transformation.

A normalising flow (Papamakarios et al., 2021) is a model that learns a
bijective transformation $f_{\phi}:z\rightarrow x$ between a base distribution
to a target distribution under maximum likelihood, where $z\sim p_{\theta}$
and $x\sim P_{X}$. The usual choice for the base distribution is a standard
normal distribution. The loss function for training this normalizing flow
$f_{\phi}$ is given by the change of variables formula

$\log P_{\theta,\phi}(x)=\log
p_{\theta}(f_{\phi}^{-1}(x))-\log\left|\det(J_{f_{\phi}^{-1}(x)})\right|,$

where $J$ is the Jacobian of $f_{\phi}$. In the conditional case this extends
to

$\log P_{\theta,\phi}(x|c)=\log
p_{\theta}\left(f_{\phi}^{-1}(x|c)\right)-\log\left|\det(J_{f_{\phi}^{-1}(x|c)})\right|,$
(1)

where $c$ are the conditional properties, and $\phi$ are the learnable
parameters of the normalizing flow $f$, and $\theta$ are the parameters of the
base distribution.

In anomaly detection methods such as (Nachman & Shih, 2020), one learns the
distribution $p_{\phi}(x|c)$ for $c$ in the sideband regions, and then queries
the conditional normalizing flow for $c$ in the signal region to obtain a
data-driven model for the background template there. This “automatic"
interpolation of the conditional density is simple and effective, however
empirically it was found in (Nachman & Shih, 2020) that for accurate
interpolation into the sideband one needed to train $p_{\phi}(x|c)$ on the
entire complement of the signal region. For a sliding window search, it is
computationally expensive to train a separate flow on the complement of every
signal region. CurtainsF4F improves on this situation by training a second
conditional flow to learn a transformation between left and right sideband
data. This flow is found to interpolate much better as the transformations to
be learnt are much simpler and this simplicity acts as an implicit
regularisation when interpolating to the signal region. One can get an
accurate background template in the signal region with just training on
narrower sidebands instead of the entire complement of the signal region. The
procedure of CurtainsF4F also allows one to train a single base flow to learn
$p_{\phi}(x|c)$ for the entire data, and then sampling from this in narrow
sidebands one can train the top flow to interpolate into any signal region.
Thus, the expensive step of training the base flow need only be done once, and
then the cheap step of training the top flow can be repeated with much less
computational cost.

$\displaystyle\max_{\gamma}\mathop{\mathbb{E}}_{x_{1},x_{2}\sim
P_{\mathrm{SB}}}\left[\log P_{\theta,\phi,\gamma}(x_{1})\right]=$
$\displaystyle\max_{\gamma}\left(\mathop{\mathbb{E}}_{x_{1},x_{2}\sim
P_{\mathrm{SB}}}\left[\log
P_{\theta,\phi}(x_{2})-\log\left|\det(J_{f_{\gamma}^{-1}(x_{1}|c_{x_{1}},c_{x_{2}})})\right|\right]\right);$
(2a) $\displaystyle\max_{\phi}\mathop{\mathbb{E}}_{x\sim P_{SB}}\left[\log
P_{\theta,\phi}(x)\right]=$
$\displaystyle\max_{\phi}\left(\mathop{\mathbb{E}}_{x\sim
P_{\mathrm{SB}}}\left[\log
p_{\theta}(f_{\phi}^{-1}(x|c_{x}))-\log\left|\det(J_{f_{\phi}^{-1}(x|c_{x})})\right|\right]\right).$
(2b)

The optimisation problem for CurtainsF4F is shown in Equation 2. Equation 2a
pertains to the optimisation of the top flow $f_{\gamma}$, while Equation 2b
pertains to the optimisation of the base flow $f_{\phi}$ that defines the log-
likelihood term in the first equation. where $p_{\mathrm{SB}}$ denotes the
distribution of features in the sidebands.
$f_{\gamma}^{-1}(x_{1}|c_{x_{1}},c_{x_{2}})=x_{2}$ is the conditional top
flow, where $x_{1}$ and $x_{2}$ are drawn from the sidebands, and $c_{x_{1}}$
and $c_{x_{2}}$ are the conditional properties of $x_{1}$ and $x_{2}$
respectively. $f_{\phi}^{-1}(x|c_{x})=z$ is the conditional base flow, where
$z\sim p_{\theta}(z)$ is drawn from a standard normal distribution. The
correspondence between the top normalizing flow and the base distributions in
Flows for Flows is shown in Figure 2.

Figure 2: The Flows for Flows architecture for a conditional model. Data $x$
($y$) are drawn from the initial distribution with conditional values $c_{x}$
($c_{y}$) and transformed to new values $c_{y}$ ($c_{x})$ in a cINN
$f_{\gamma(c_{x},c_{y})}$ conditioned on $c_{x}$ and $c_{y}$. The probability
of the transformed data points are evaluated using a second normalizing flow
for the base distribution $f_{\phi^{\prime}(c_{y})}$ ($f_{\phi(c_{x})}$).
$Z_{1},\mathrm{and}\,Z_{2}$ denote analytically known distributions, for
example, a standard normal. In the case where $x$ and $y$ are drawn from the
same underlying distribution $p(x,c)$, the same base distribution $f_{\phi}$
can be used.

CurtainsF4F is trained in both directions. The forward pass transforms data
from low to higher target values of proper motion, whereas the inverse pass
transforms data from high to lower target values. Data are drawn from both SBs
and target proper motion values are randomly assigned to each data point using
all proper motion values in the batch. Data are passed through the network in
a forward or inverse pass, depending on whether the proper motion is larger or
smaller than their initial proper motion. The network is conditioned on a
function of initial and target proper motion, with the two values ordered in
ascending order. This function could be, for example, difference between the
two, or simply both values concatenated. The probability term is evaluated
using a single base distribution trained on the data from SB1 and SB2. The
loss for the batch is calculated from the average of the probabilities
calculated from the forward and inverse passes. A schematic overview is shown
in Figure 3.

Figure 3: A schematic overview of the training procedure for CurtainsF4F with
an event where the target proper motion (c) value is greater than the input
value. A single conditional normalizing flow is used for the base
distribution, conditioned on the target c value $c_{\textit{target}}$, to
determine $p_{\theta}(z|c_{\textit{target}})$. The top normalizing flow is
conditioned on a function of the input ($c_{\textit{input}}$) and target
($c_{\textit{target}}$) proper motion values. For the case where
$c_{\textit{target}}<c_{\textit{input}}$, an inverse pass of the network is
used, and the conditioning property is calculated as
$f(c_{\textit{target}},c_{\textit{input}})$.

The base flow is trained on the sideband data with a standard normal
distribution as the target prior. It is conditioned on the proper motion. The
top flow is trained between data drawn from the sidebands. The transformation
is conditioned on the concatenated tuple of the initial and target proper
motion. Depending on which proper motion is chosen as the conditioning
feature, the downstream task of finding the stream might be affected. This is
because the stream candidate stars may have a non-trivial correlation, and
therefore may produce overdensities of different shapes in the two proper
motions respectively. In this work, we use $\mu_{\lambda}$ as the conditional
feature. The rest of features, i.e. [$\phi$, $\lambda$, $G$,
$G_{\text{BP}}-G_{\text{RP}}$, $\mu_{\phi}\cos\lambda$ ] are used to
characterise the template.

One important aspect of the CurtainsF4F method is the definition of signal and
sideband regions. In Figure 4, we show the distribution of $\mu_{\lambda}$ of
the background stars and GD-1 stream stars in the sidebands (SB1, SB2) and
signal region (SR). Unlike in the Via Machinae method, where the sideband
region was the complementary region of the chosen SR, SkyCURTAINs defines the
sideband region to be typically of 2-6 mas/yr. Since the top flow only needs
to learn a small (but not necessarily trivial) transformation of the sideband
data to generate a template in the SR, we find this width of the sideband to
be sufficient. This also cuts down the total training time compared to Via
Machinae, as despite both methods consisting of two generative models,
SkyCURTAINs’s second generative model effectively learns on narrower
sidebands. To demonstrate the efficacy of the method on the GD-1 stream, we
define the signal region as the interval in which the signal is contained. In
an actual analysis, where the location of the signal is not known a priori,
one would need to scan multiple values of $\mu_{\lambda}$ with the CurtainsF4F
method. Here, the modularity of the CurtainsF4F method comes into play, as the
base flow can be trained on the entire patch of the sky and frozen.
Thereafter, top flows can be trained on individual regions of interest. This
significantly reduces the computational cost of training the model, by
allowing the base flow to be trained once and reused for multiple regions of
interest.

Once the CurtainsF4F model is trained, the background-enriched template is
constructed by transforming the data from the sidebands to the signal region,
conditioned on the tuple of the initial and target proper motion. The target
proper motions in the signal region are sampled from a kernel density fit on
the proper motion distribution. We train an ensemble of $10$ multi layer
perceptron based classifiers111Details about the architecture can be found in
Appendix A between this template and the signal region data. The scores are
then aggregated by taking the mean score and used to classify the data in the
signal region. The top 0.1% most signal-like stars are selected as the
candidates for the next step. This threshold is selected to remove the most
background like stars from the signal region, which is a tunable parameter and
can be adjusted based on the desired purity and signal efficiency.

Figure 4: Distribution of $\mu_{\lambda}$ of background stars (in blue) and
GD-1 stream stars (in red) in the sidebands and signal region defined in bins
of $\mu_{\lambda}$. The signal region is defined as $\mu_{\lambda}$
$\in\left[-14,-8\right)$, and the sidebands are defined as $\mu_{\lambda}$
$\in\left[-18,-14\right)$ and $\mu_{\lambda}$ $\in\left[-8,-6\right)$. In this
SR, there are 184142 background like stars, and 688 PWB18 tagged GD-1 stream
stars, which is a signal fraction of only $0.0037$.

###  3.2 Line detection

The CurtainsF4F step gives us a set of stars which produce an overdensity in
the feature space. We still need to filter out the overdensities that are
particularly line like, as we are interested in stellar streams. We employ a
well known line finding algorithm to estimate the line parameters of the
stream via the Hough transform, as was done in (Shih et al., 2021, 2023).
Since we do not apply a fiducial cut to eliminate stars outside a $10^{\circ}$
radius, we can use the full set of stars that pass the CurtainsF4F cut in the
patch to estimate the line parameters.

For a given star located at ($\phi$ ′, $\lambda$ ′), the Hough transform is a
mapping from the $\phi$-$\lambda$ space to the parameter space $\rho-\theta$,
defined as:

$\rho=\phi^{\prime}\cos\theta-\lambda^{\prime}\sin\theta$ (3)

where $\rho$ is the perpendicular distance from the origin to the line and
$\theta$ is the angle between the line passing through ($\phi$ ′, $\lambda$ ′)
and the $\phi$-axis. The origin is defined as the point $(0,0)$ in the
$\phi$-$\lambda$ space. A line in the $\phi$-$\lambda$ space is represented as
a point in the $\rho-\theta$ space. If a set of points lie on a line in the
$\phi$-$\lambda$ space, they will map to a single point in the $\rho-\theta$
space. Therefore, identifying line candidates in the $\phi$-$\lambda$ space is
equivalent to identifying high density regions in the $\rho-\theta$ space. For
each patch, we bin the Hough space in a $100\times 100$ grid from
$-15^{\circ}\leq\rho\leq 15^{\circ}$ and $0\leq\theta\leq\pi$, such that each
bin $(i,j)\left(i\in[0,99],j\in[0,99]\right)$ is related to the $\rho$ and
$\theta$ values by:

$\displaystyle\rho_{i}$ $\displaystyle=$
$\displaystyle-15+i\times\frac{30}{99}$ (4) $\displaystyle\theta_{j}$
$\displaystyle=$ $\displaystyle\pi\times\frac{j}{99}$ (5)

Figure 5: Left: The significance map of the Hough space for candidate signal
stars in the patch with coordinates $\alpha$ = $146.9^{\circ}$, and $\delta$ =
$35.6^{\circ}$. Each pixel represents the number of stars whose Hough curve
passes through a box of width
$(\Delta\rho=1.5^{\circ},\Delta\theta=0.03\,\mathrm{rad})$ centred at that
pixel. The bright spots correspond to the peaks in the Hough space. The
highest significance pixel is marked with a black cross. Right:
$\phi$-$\lambda$ scatter plot of stars post CWoLa step in the same patch. The
red coloured stars correspond to the GD-1 stream as identified by PWB18. The
grey coloured stars are those selected by SkyCURTAINs as the most signal like
stars in this patch. The blue coloured stars are those selected after applying
a hough filter and fall within the acceptance region defined by the Hough fit
lines in purple.

Stellar streams have a finite width in space. To account for this, we
calculate the number of lines passing through a box of width ($\Delta
i=5,\Delta j=3$) centred at $(i,j)$. This is done by convolving the Hough
space with a uniform kernel of size $(\Delta i,\Delta j)$. These widths
correspond to $\Delta\rho=1.5^{\circ}$ and $\Delta\theta=0.09$ rad.
Thereafter, similar to (Shih et al., 2023), we scan this convolved Hough space
for high significance peaks. This allows us to calculate a range of $\rho$ and
$\theta$ values that correspond to the line-like structures in the
$\phi$-$\lambda$ space, and account for the finite width of the stream. Figure
5 (left) shows the Hough space for the GD-1 stream in one of the patches.

##  4 Results

The CurtainsF4F stage was trained on NVIDIA® RTX 3080 GPUs, and the Hough
stage was run on a single CPU core. The CurtainsF4F stage took $\sim 4$ hours
per patch, amounting to a total of $\sim 80$ GPU hours. The CurtainsF4F stage
consists of training a base and a top flow. Both base and top flow took $\sim
2$ hours to train. The line fitting took about a minute per patch, and was a
negligible fraction of the total computational cost.

Figure 6: Feature correlation plots in SB1 (orange), SR (red), and SB2
(green). The diagonal panels show the marginals of the features in the SB1
(orange), SR (red), and SB2 (green). The off-diagonal panels show the
correlation between the features. The left panel shows the feature correlation
plots of true sideband region data and true signal region data. The right
panel shows the feature correlation plots of the CurtainsF4F generated
template (SB1 $\rightarrow$ SR, SB2 $\rightarrow$ SR), and true signal region
data. The generated template has a much better agreement with the true signal
region data, and preserves the correlation between the features. The features
correspond to the patch centred at $\alpha$ = $146.9^{\circ}$, and $\delta$ =
$35.6^{\circ}$

.

The most crucial step in the SkyCURTAINs method is the generation of a
background enriched template in the signal region. In Figure 6, we show the
marginals and correlations of features in the sidebands and signal region in
the left panel. The features are strongly correlated with the proper motion,
which would bias the classifier in the CWoLa step to produce false positives
in the signal region even in the absence of a stream. In the right panel, we
show the marginals and correlations of the features in the generated template
by CurtainsF4F in the same patch. The generated template leverages the
correlation of the features with the proper motion to construct a background
enriched template in the signal region. This allows for a more representative
template of the background in the signal region, and reduces the false
positives in the search for stellar streams. With the generated template, we
can now train a classifier in the CWoLa step to tag the stars in the signal
region.

###  4.1 Metrics

We now demonstrate the performance of the SkyCURTAINs method on GDR2 data. To
quantify the discovery potential of SkyCURTAINs method, we measure the
Significance Improvement Characteristic (SIC) curve for the GD-1 stream. In
Figure 7, we show the SIC curve as a function of the signal efficiency for the
GD-1 stream in one of the 21 patches. This metric is defined as the ratio of
the signal efficiency to the square root of the background efficiency, and
essentially quantifies the improvement in the discovery significance of the
signal from the method. SkyCURTAINs achieves a maximum significance
improvement of $\sim 10$ at $\sim 50\%$ signal efficiency. Although a direct
comparison with Via Machinae is difficult on account of different SR being
used for the analysis, one can look at the maximum value of the SIC as a
heuristic measure, which are comparable for both methods.

Figure 7: Significance improvement characteristic curve as a function of
signal efficiency for the GD-1 stream in the patch with coordinates $\alpha$ =
$146.9^{\circ}$, and $\delta$ = $35.6^{\circ}$.

We track two other metrics to quantify the performance of SkyCURTAINs: purity
$p$: The fraction of candidate CurtainsF4F stars that overlap with the PWB18
identified GD-1 stream members; and signal efficiency, $\epsilon_{S}$ which is
the fraction of GD-1 stream members that have been flagged as candidates by
CurtainsF4F step. Figure 5 (right) shows the candidates from the CurtainsF4F
step in the $\phi$-$\lambda$ space that corresponds to the GD-1 stream with a
$p=75\%$ and $\epsilon_{S}=36.82\%$. We note that it also predicts a few stars
that do not form a line like structure in the $\phi$-$\lambda$ space. This is
expected, as this stage is designed to flag any overdensity in the feature
space as a potential signal candidate. To filter out the line like
overdensities we perform a Hough transform on the output of CurtainsF4F step.
After applying the Hough filter, the purity is improved to $91.79\%$, albeit
at the cost of a slightly reduced signal efficiency of $34.3\%$.

###  4.2 Full GD-1 stream scan

Figure 8: Top: Comparison of PWB18 and the 753 SkyCURTAINs identified stream
candidates. The first panel shows the PWB18 GD-1 stream members in the GD-1
stream stream aligned coordinate system $(\phi_{1},\phi_{2})$. The middle
panel shows the potential members of the GD-1 stream stream in the GD-1 stream
aligned coordinate system identified by the SkyCURTAINs method. The bottom
panel is a comparison of the number of candidate stream stars identified by
the SkyCURTAINs method (red) and the number of PWB18 (black) in the GD-1
stream in $\phi_{1}$ bins of width $1^{\circ}$. The shaded region to the left
of the stream correspond to the patches which are excluded from the analysis
due to their proximity to the galactic disk. The shaded region to the right of
the stream correspond to the patches containing stars with very low proper
motions, which reduces the signal sensitivity of the CurtainsF4F to potential
streams. Bottom: Coverage plots in the $G$ vs $G_{\text{BP}}-G_{\text{RP}}$
and $\mu_{\alpha}$ vs $\mu_{\alpha}\cos\delta$ space. The contours show the
68.2, 95.4 and 99.7 percentiles of the GD-1 stream members identified by
PWB18. The SkyCURTAINs identified candidates are shown in red.

Table 1 shows the performance of the SkyCURTAINs method in the 21 patches that
contain the GD-1 stream. We quote the purity $p$ after applying the Hough
filter for each patch. SkyCURTAINs is able to identify the GD-1 stream members
with a high purity in most of the patches, and significantly improves the
performance compared to standalone CWoLa. In Table 2 we report the total PWB18
identified GD-1 stream members and SkyCURTAINs candidates (after Hough filter)
in the patches. The combined result is shown in Figure 8.

Table 1: Performance of the SkyCURTAINs method in the 21 patches that contain the GD-1 stream. The patches are identified by the central $\alpha$ and $\delta$ of the patch. We quote the purity $p$ after applying the Hough filter for each patch, and compare the performance with standalone CWoLa.  Patch ($\alpha$, $\delta$) | $p$
---|---
SkyCURTAINs | CWoLa
$\left(128.4^{\circ},\,\,28.8^{\circ}\right)$ | $\mathbf{82.99}$ | $77.0$
$\left(132.6^{\circ},\,\,16.9^{\circ}\right)$ | $\mathbf{78.05}$ | $62.0$
$\left(136.5^{\circ},\,\,36.1^{\circ}\right)$ | $\mathbf{90.56}$ | $86.0$
$\left(138.8^{\circ},\,\,25.1^{\circ}\right)$ | $\mathbf{90.79}$ | $84.0$
$\left(142.7^{\circ},\,\,14.5^{\circ}\right)$ | $\mathbf{86.79}$ | $65.0$
$\left(146.9^{\circ},\,\,35.6^{\circ}\right)$ | $\mathbf{91.79}$ | $90.0$
$\left(148.6^{\circ},\,\,24.2^{\circ}\right)$ | $\mathbf{94.9}$ | $87.0$
$\left(148.6^{\circ},\,\,47.0^{\circ}\right)$ | $\mathbf{93.15}$ | $78.0$
$\left(156.2^{\circ},\,\,57.5^{\circ}\right)$ | $\mathbf{70.14}$ | $54.0$
$\left(156.9^{\circ},\,\,34.1^{\circ}\right)$ | $\mathbf{88.17}$ | $86.0$
$\left(160.5^{\circ},\,\,45.5^{\circ}\right)$ | $\mathbf{87.43}$ | $73.0$
$\left(171.4^{\circ},\,\,43.0^{\circ}\right)$ | $\mathbf{89.52}$ | $72.0$
$\left(171.8^{\circ},\,\,54.7^{\circ}\right)$ | $\mathbf{89.66}$ | $53.0$
$\left(174.3^{\circ},\,\,65.1^{\circ}\right)$ | $\mathbf{64.94}$ | $47.0$
$\left(185.4^{\circ},\,\,50.0^{\circ}\right)$ | $\mathbf{84.0}$ | $57.0$
$\left(192.0^{\circ},\,\,58.7^{\circ}\right)$ | $\mathbf{83.87}$ | $66.0$
$\left(138.1^{\circ},\,\,5.7^{\circ}\right)$ | $0.0$ | $0.0$
$\left(203.7^{\circ},\,\,49.1^{\circ}\right)$ | $\mathbf{0.13}$ | $0.0$
$\left(212.7^{\circ},\,\,55.2^{\circ}\right)$ | $0.0$ | $0.0$
$\left(224.7^{\circ},\,\,60.6^{\circ}\right)$ | $\mathbf{2.58}$ | $0.0$
$\left(202.4^{\circ},\,\,66.5^{\circ}\right)$ | $0.0$ | $\mathbf{50.0}$
Table 2: PWB18 identified GD-1 stream members and SkyCURTAINs candidates in the 21 patches that contain the GD-1 stream.  Patch ($\alpha$, $\delta$) | PWB18 |  SkyCURTAINs
---|---|---
$\left(128.4^{\circ},\,\,28.8^{\circ}\right)$ | $307$ | $147$
$\left(132.6^{\circ},\,\,16.9^{\circ}\right)$ | $321$ | $41$
$\left(136.5^{\circ},\,\,36.1^{\circ}\right)$ | $428$ | $180$
$\left(138.8^{\circ},\,\,25.1^{\circ}\right)$ | $421$ | $76$
$\left(142.7^{\circ},\,\,14.5^{\circ}\right)$ | $312$ | $106$
$\left(146.9^{\circ},\,\,35.6^{\circ}\right)$ | $564$ | $207$
$\left(148.6^{\circ},\,\,47.0^{\circ}\right)$ | $470$ | $73$
$\left(148.6^{\circ},\,\,24.2^{\circ}\right)$ | $415$ | $98$
$\left(156.2^{\circ},\,\,57.5^{\circ}\right)$ | $473$ | $144$
$\left(156.9^{\circ},\,\,34.1^{\circ}\right)$ | $541$ | $186$
$\left(160.5^{\circ},\,\,45.5^{\circ}\right)$ | $584$ | $167$
$\left(171.4^{\circ},\,\,43.0^{\circ}\right)$ | $551$ | $229$
$\left(171.8^{\circ},\,\,54.7^{\circ}\right)$ | $585$ | $116$
$\left(174.3^{\circ},\,\,65.1^{\circ}\right)$ | $453$ | $77$
$\left(185.4^{\circ},\,\,50.0^{\circ}\right)$ | $551$ | $50$
$\left(192.0^{\circ},\,\,58.7^{\circ}\right)$ | $583$ | $62$
$\left(138.1^{\circ},\,\,5.7^{\circ}\right)$ | $209$ | $3$
$\left(203.7^{\circ},\,\,49.1^{\circ}\right)$ | $380$ | $4$
$\left(212.7^{\circ},\,\,55.2^{\circ}\right)$ | $336$ | $0$
$\left(224.7^{\circ},\,\,60.6^{\circ}\right)$ | $244$ | $28$
$\left(202.4^{\circ},\,\,66.5^{\circ}\right)$ | $380$ | $0$
Figure 9: GD-1 stream purity obtained by the SkyCURTAINs method as a function
of PWB18 signal to background ratio (shown here in percentages).

SkyCURTAINs has a very low purity in 5 of the 21 patches of the sky where the
GD-1 stream is present. On closer inspection we find that in patches
($\alpha$, $\delta$) =
$\left[\left(203.7^{\circ},49.1^{\circ}\right),\left(212.7^{\circ},55.2^{\circ}\right),\left(224.7^{\circ},60.6^{\circ}\right),\left(202.4^{\circ},66.5^{\circ}\right)\right]$,
the GD-1 stream members peak at very low proper motion ($\mu_{\lambda}$). This
results in a SR that is dominated by distant stars, and the sensitivity of the
CurtainsF4F step to actual stream stars is reduced. These patches correspond
to $\phi_{1}\geq-10^{\circ}$ in the GD-1 stream aligned coordinates, which
explains the low yield of the SkyCURTAINs method to the right of the stream in
Figure 8. SkyCURTAINs also has a low purity in the patch centred at ($\alpha$,
$\delta$) = $\left(138.1^{\circ},5.7^{\circ}\right)$. The low GD-1 stream
purity in this patch is likely due to the extremely low signal to background
ratio in the corresponding SR. In Figure 9 we show the GD-1 stream purity as a
function of the PWB18 signal to background ratio. We find the purity has a
sharp drop to zero when the signal to background ratio is near 0.01%. Patch
($\alpha$, $\delta$) = $\left(138.1^{\circ},5.7^{\circ}\right)$ has a signal
to background ratio of 0.01%, which is the lowest in the 21 patches. This
patch corresponds to $-80^{\circ}\leq\phi_{1}\leq-60^{\circ}$ in the GD-1
stream aligned coordinates, and explains the low yield of the SkyCURTAINs
method to the left of the stream in Figure 8. These patches (marked in red)
are the patches where the GD-1 stream members peak at very low proper motion
($\mu_{\lambda}$), and the sensitivity of the CurtainsF4F step to actual
stream stars is reduced. This patch corresponds to
$-80^{\circ}\leq\phi_{1}\leq-60^{\circ}$ in the GD-1 stream aligned
coordinates, and explains the low yield of the SkyCURTAINs method to the left
of the stream in Figure 8.

It is crucial to note that SkyCURTAINs method assumes very little
astrophysical information about the stream, allowing it to be agnostic to the
stream’s properties. The only information used in the method is the proper
motion which is used to define the SR and SB regions. For stream
identification, fiducial cuts on $G_{\text{BP}}-G_{\text{RP}}$ and $G$ (there
no requirements on streams to lie on an isochrone) are applied. This is in
parity with the fiducial cuts applied in (Pettee et al., 2023; Shih et al.,
2021). SkyCURTAINs flags 753 unique stars as potential GD-1 stream members, of
which 568 are also identified by PWB18, thereby attaining an overall GD-1
stream purity of $75.4\%$. This surpasses the standalone CWoLa method which
has a purity of $56\%$, and Via Machinae 1.0 which has a purity of $49\%$.
SkyCURTAINs also outperforms Via Machinae 2.0, which has a purity of $65\%$,
despite the latter employing additional fiducial cuts and performs an
augmented scan over both proper motions. There are 1498 PWB18 identified GD-1
stream stars in our fiducial region, which gives us a global signal efficiency
of $37.9\%$. Furthermore, an important result of the SkyCURTAINs method is
that it produces no spurious streams in the 21 patches that were scanned. This
can be attributed to the very stringent selection criteria applied in the
CurtainsF4F stage of the method, designed to reduce false positives.

Figure 10: Coverage plot for the GD-1 stream in the $G$ vs
$G_{\text{BP}}-G_{\text{RP}}$ space. The contours show the 68.2, 95.4 and 99.7
percentiles of the GD-1 stream members identified by PWB18. The 185 additional
SkyCURTAINs identified candidates are shown in red.

Of the remaining 185 stars, some may potentially be new undiscovered members
of the GD-1 stream. Figure 10 shows the isochrone plot for the GD-1 stream
members identified by PWB18, with the additional SkyCURTAINs candidates
overlaid. There is a significant overlap between these 185 stars and the PWB18
labelled members, which suggests that the SkyCURTAINs method is able to
identify some members of the GD-1 stream that may have been missed by PWB18.
There are also a few stars that are not part of the GD-1 stream isochrone, and
are likely to be false positives.

Despite the lack of prior astrophysical information, the SkyCURTAINs method is
able to recover well known density perturbations in the GD-1 stream. In the
GD-1 stream stream aligned coordinates $(\phi_{1},\phi_{2})$ (Koposov et al.,
2010) shown in Figure 8, we see that SkyCURTAINs recovers the "gaps" at
$\phi_{1}\approx-40^{\circ}$ and $\phi_{1}\approx-20^{\circ}$, as well as the
"offshoot" or "spur" at $\phi_{1}\approx-35^{\circ}$, which are well known
features of the GD-1 stream. Furthermore, SkyCURTAINs predictions of the
overdensity regions at $\phi_{1}\approx-50^{\circ}$ and
$\phi_{1}\approx-10^{\circ}$ are in good agreement with the PWB18 members. The
low yield regions at $\phi_{1}\geq-10^{\circ}$ and
$-80^{\circ}\leq\phi_{1}\leq-60^{\circ}$ are due to the reasons discussed
above. The region $\phi_{1}\leq-80^{\circ}$ correspond to the patches that are
excluded from the analysis due to their proximity to the galactic disk.

##  5 Conclusion

In this work, we described the SkyCURTAINs method, a model-agnostic, template
based, data-driven approach to detecting stellar streams in the Milky Way
using the _Gaia_ DR2 data. Originally developed for anomaly detection in High
Energy Physics, SkyCURTAINs joins the ranks of Via Machinae and CWoLa, in the
search for stellar streams, which highlights the versatility of these tools in
their performance across different domains. Synergies between the High Energy
Physics, Astrophysics, and other communities should be further encouraged in
order to identify similar problems that can be solved using the same tools
developed in respective fields.

We demonstrated the performance of SkyCURTAINs on the GD-1 stream, and its
ability to identify the line-like overdensity in the $\phi$-$\lambda$ space
that corresponds to the GD-1 stream with a very high purity across most
patches. The main advantage of SkyCURTAINs is the minimal assumptions it makes
about the underlying signal, thereby making it a versatile tool for
identifying any localized overdensities (streams, globular clusters, dwarf
galaxies) in the feature space of the stars in a model agnostic manner. In
this work we chose GD-1 stream, as it provides a good test case for
SkyCURTAINs, and is a well known stream in the Milky Way, where we show that
SkyCURTAINs currently outperforms the other weakly supervised machine learning
methods like Via Machinae and CWoLa in terms of purity by over 10%. For a full
sky scan for streams, where the locations of the streams are unknown, the
method will need to scan a larger number of $\mu_{\lambda}$. In this case,
training two conditional generative models for each patch may not be feasible.
However, the modularity of the design of CurtainsF4F step allows for much
easier scaling. The base flow can be trained on the entire patch and frozen
and then individual top flows can be trained on respective regions of
interest, for efficient scaling. SkyCURTAINs builds the template by leveraging
the correlations of features with the proper motions, which results in a more
background representative template. This allows the method to use more
discriminatory features to be used in the downstream task regardless of their
correlation with the proper motions, which otherwise would lead to a biased
classifier and thus false positives.

The follow-up work will involve applying SkyCURTAINs in a full sky search for
stellar streams in the Milky Way, and comparing the performance of the method
with other existing methods. The latest data release GDR3 from _Gaia_ contains
over 1.8 billion sources with improved astrometric and photometric
measurements for a significantly larger number of sources compared to GDR2,
which will provide a more detailed view of the Milky Way. The improved
measurements of radial velocities of sources in GDR3 will also allow for a
more detailed study of the kinematics of the streams. The improvements in data
quality in GDR3 would likely further improve the performance of SkyCURTAINs
and as such is a natural next step for the method.

## Acknowledgements

We would like to acknowledge funding through the SNSF Sinergia grant
CRSII5_193716 “Robust Deep Density Models for High-Energy Particle Physics and
Solar Flare Analysis (RODEM)”. This work has made use of data from the
European Space Agency (ESA) mission _Gaia_ (https://www.cosmos.esa.int/gaia),
processed by the _Gaia_ Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. The computations were
performed at University of Geneva using Baobab HPC service. Special thanks to
Samuel Klein and Kinga Anna Wozniak for their inputs during the development of
the method and this manuscript.

## Data Availability

SkyCURTAINs uses the publicly available GDR2 data. A curated set of 21 patches
used in this work is available at https://zenodo.org/records/7897936. The GD-1
stream membership labels are taken from https://zenodo.org/records/1295543.

## References

  * Banik & Bovy (2019) Banik N., Bovy J., 2019, Monthly Notices of the Royal Astronomical Society, 484, 2009
  * Belokurov et al. (2006) Belokurov V., et al., 2006, The Astrophysical Journal, 642, L137–L140
  * Bonaca et al. (2019) Bonaca A., Hogg D. W., Price-Whelan A. M., Conroy C., 2019, The Astrophysical Journal, 880, 38
  * Bonaca et al. (2020) Bonaca A., et al., 2020, The Astrophysical Journal Letters, 892, L37
  * Borsato et al. (2019) Borsato N. W., Martell S. L., Simpson J. D., 2019, Monthly Notices of the Royal Astronomical Society, 492, 1370–1384
  * Carlberg (2017) Carlberg R. G., 2017, The Astrophysical Journal, 838, 39
  * Carlberg et al. (2012) Carlberg R. G., Grillmair C. J., Hetherington N., 2012, The Astrophysical Journal, 760, 75
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Golling et al. (2023) Golling T., Klein S., Mastandrea R., Nachman B., Raine J. A., 2023, Phys. Rev. D, 108, 096018
  * Grillmair & Dionatos (2006) Grillmair C. J., Dionatos O., 2006, The Astrophysical Journal, 643, L17
  * Helmi & White (1999) Helmi A., White S. D. M., 1999, Monthly Notices of the Royal Astronomical Society, 307, 495
  * Hough (1962) Hough P. V., 1962, Method and means for recognizing complex patterns
  * Ibata et al. (2001) Ibata R., Lewis G. F., Irwin M., Totten E., Quinn T., 2001, The Astrophysical Journal, 551, 294–311
  * Ibata et al. (2021) Ibata R., et al., 2021, The Astrophysical Journal, 914, 123
  * Johnston (1998) Johnston K. V., 1998, The Astrophysical Journal, 495, 297–308
  * Johnston et al. (1999) Johnston K. V., Zhao H., Spergel D. N., Hernquist L., 1999, The Astrophysical Journal, 512, L109–L112
  * Koposov et al. (2010) Koposov S. E., Rix H.-W., Hogg D. W., 2010, The Astrophysical Journal, 712, 260–273
  * Malhan & Ibata (2018) Malhan K., Ibata R. A., 2018, Monthly Notices of the Royal Astronomical Society, 477, 4063–4076
  * Malhan et al. (2018) Malhan K., Ibata R. A., Goldman B., Martin N. F., Magnier E., Chambers K., 2018, Monthly Notices of the Royal Astronomical Society, 478, 3862–3870
  * Meingast, Stefan & Alves, João (2019) Meingast, Stefan Alves, João 2019, A&A, 621, L3
  * Meingast, Stefan et al. (2019) Meingast, Stefan Alves, João Fürnkranz, Verena 2019, A&A, 622, L13
  * Metodiev et al. (2017) Metodiev E. M., Nachman B., Thaler J., 2017, Journal of High Energy Physics, 2017
  * Nachman & Shih (2020) Nachman B., Shih D., 2020, Physical Review D, 101
  * Necib et al. (2019) Necib L., Lisanti M., Garrison-Kimmel S., Wetzel A., Sanderson R., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2019, The Astrophysical Journal, 883, 27
  * Papamakarios et al. (2021) Papamakarios G., Nalisnick E., Rezende D. J., Mohamed S., Lakshminarayanan B., 2021, Journal of Machine Learning Research, 22, 1
  * Pettee et al. (2023) Pettee M., Thanvantri S., Nachman B., Shih D., Buckley M. R., Collins J. H., 2023, Weakly-Supervised Anomaly Detection in the Milky Way (arXiv:2305.03761)
  * Price-Whelan & Bonaca (2018) Price-Whelan A. M., Bonaca A., 2018, The Astrophysical Journal Letters, 863, L20
  * Purcell et al. (2012) Purcell C. W., Zentner A. R., Wang M.-Y., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 027–027
  * Raine et al. (2023) Raine J. A., Klein S., Sengupta D., Golling T., 2023, Frontiers in Big Data, 6
  * Sanders & Binney (2013) Sanders J. L., Binney J., 2013, Monthly Notices of the Royal Astronomical Society, 433, 1813
  * Sanders et al. (2016) Sanders J. L., Bovy J., Erkal D., 2016, Monthly Notices of the Royal Astronomical Society, 457, 3817–3835
  * Sengupta et al. (2023) Sengupta D., Klein S., Raine J. A., Golling T., 2023, CURTAINs Flows For Flows: Constructing Unobserved Regions with Maximum Likelihood Estimation (arXiv:2305.04646)
  * Shih et al. (2021) Shih D., Buckley M. R., Necib L., Tamanas J., 2021, Monthly Notices of the Royal Astronomical Society, 509, 5992
  * Shih et al. (2023) Shih D., Buckley M. R., Necib L., 2023, Via Machinae 2.0: Full-Sky, Model-Agnostic Search for Stellar Streams in Gaia DR2 (arXiv:2303.01529)
  * Varghese et al. (2011) Varghese A., Ibata R., Lewis G. F., 2011, Monthly Notices of the Royal Astronomical Society, 417, 198–215
  * Vera-Casanova et al. (2022) Vera-Casanova A., et al., 2022, Monthly Notices of the Royal Astronomical Society, 514, 4898
  * Yuan et al. (2018) Yuan Z., Chang J., Banerjee P., Han J., Kang X., Smith M. C., 2018, The Astrophysical Journal, 863, 26

##  Appendix A CurtainsF4F Training and Hyperparameter Tuning Details

###  A.1 CurtainsF4F features preprocessing

The first step in training the CurtainsF4F model is to define the SR and SB
regions in $\mu_{\lambda}$. In this work, we chose the SR in a given patch to
ensure that the GD-1 stream is fully contained. The SB region, defined as the
region adjacent to the SR is chosen to be $\sim 6$ mas/yr wide. This ensures
sufficient training statistics for the base and the top flow model. The
features used for training the CurtainsF4F model are: [$\phi$, $\lambda$, $G$,
$G_{\text{BP}}-G_{\text{RP}}$, $\mu_{\phi}^{*}$ ] and the conditional feature
is $\mu_{\lambda}$. As these features have different dynamic ranges, we opt to
further scale them to ensure a stable model training. All features are first
scaled to be in the range $[0,1]$. The $G$ feature has a sharp cutoff at 20.2
which proves to be a difficult feature for generative models to learn. To
mitigate this, we apply a logit transformation to the $G$ feature. The logit
transformation is defined as:

$\text{logit}(x)=\log\left(\frac{x}{1-x}\right),$ (6)

Finally, all features are scaled to be in the range $[-3,3]$.

The data for the base and top flow training is divided into training and
validation sets in a $80:20$ ratio.

###  A.2 Hyperparameter tuning

The three neural network components in the SkyCURTAINs method are the base
flow, the top flow, and the classifier in the CWoLa step. The base flow
comprises of a stack of autoregressive transformations parametrised by a
Rational Quadratic Spline (RQS) function. The top flow is a stack of coupling
transformations, also parametrised by a RQS function. The hyperparameters for
the base and the top flow are the number of stacked transformations, the
number of bins in the RQS function, and the number of hidden units and layers
in the multi-layer perceptron (MLP) used to estimate the parameters of the RQS
function. Other hyperparameters include the learning rate, the batch size, and
the number of epochs for training the base and top flow. The classifier in the
CWoLa step is an MLP with 3 hidden layers of 32 units, and are trained using
the Adam optimizer with a learning rate of $10^{-3}$ with k-Fold cross
validation (with k = 5). We found this architecture to be robust across
different patches, and did not perform any hyperparameter scans. For the
CurtainsF4F training, we want to ensure that the generated template is in
accordance with data. Since we do not know a priori if there is a signal in
the SR, we test the performance of the CurtainsF4F model in the sidebands. We
select the hyperparameters that minimise the AUC score for the CWoLa
classifier on SB1 and SB2 vs template classification. For a well-trained
CurtainsF4F model the generated template should have an AUC score close to
0.5.

The hyperparameters for the base and top flow are listed in Table 3, where the
hyperparameters for the base flows were found to give robust performance
regardless of the patch, and so held constant. For the top flows, there could
be significant variation in performance related to the hyperparameter
selection depending on the patch, and so hyperparameter tuning was performed
to find the values that performed well regardless of the patch. Both base and
top flow were capped at a maximum number of 150, and 100 epochs respectively.
While the base flow seemed to improve with higher number of epochs, the top
flow converged much more quickly at $\sim 30-40$ epochs of training.

Table 3: Hyperparameters for CurtainsF4F Training Model | Hyperparameter Name | Value
---|---|---
Base Flow | Number of Stacked Transformations | 4
Number of Bins in RQS Function | 12
Number of Hidden Units in MLP | 64
Number of Layers in MLP | 2
Learning Rate | 0.0001
Batch Size | 512
Maximum number of Epochs | 150
Top Flow | Number of Stacked Transformations | 6
Number of Bins in RQS Function | 10
Number of Hidden Units in MLP | 32
Number of Layers in MLP | 2
Learning Rate | 0.001
Batch Size | 512
Maximum number of Epochs | 100


# SISE-PC: Semi-supervised Image Subsampling for Explainable Pathology
Classification

Sohini Roychowdhury  Director, Curriculum,
Fourthbrain.ai<EMAIL_ADDRESS>Kwok Sun Tang  University of Illinois
<EMAIL_ADDRESS>Mohith Ashok  AggDirect
<EMAIL_ADDRESS>Anoop Sanka  Fourthbrain.ai
<EMAIL_ADDRESS>
###### Abstract

Although automated pathology classification using deep learning (DL) has
proved to be predictively efficient, DL methods are found to be data and
compute cost intensive. In this work, we aim to reduce DL training costs by
pre-training a Resnet feature extractor using SimCLR contrastive loss for
latent encoding of OCT images. We propose a novel active learning framework
that identifies a minimal sub-sampled dataset containing the most uncertain
OCT image samples using label propagation on the SimCLR latent encodings. The
pre-trained Resnet model is then fine-tuned with the labelled minimal sub-
sampled data and the underlying pathological sites are visually explained. Our
framework identifies upto 2% of OCT images to be most uncertain that need
prioritized specialist attention and that can fine-tune a Resnet model to
achieve upto 97% classification accuracy. The proposed method can be extended
to other medical images to minimize prediction costs.

###### Index Terms:

SimCLR, contrastive loss, Resnet, semi-supervised, label spreading

##  I Introduction

Automated pathology classification has shown to significantly improve patient
prioritization and resourcefulness of treatment procedures and patient care
[1]. Although deep learning algorithms such as Resnet and InceptionV3 have
been established as state-of-the-art [2] for several pathology classification
tasks, training these models from scratch can be expensive from the labelled
data acquisition and compute resource perspectives. In this work, we present a
semi-supervised image sub-sampling method that identifies a minimal sub-
sampled data set that represents the most sample uncertainty in a latent
feature space that is encoded using a self-supervised contrastive model [3].
We demonstrate the classification performance of a pre-trained Resnet model
that is fine tuned using only the minimal sub-sampled data for multi-class
pathology classification.

Optical Coherence Tomography (OCT) is a commonly performed diagnostic test
designed to assist doctors in identifying retinal diseases, such as choroidal
neovascularization (CNV), Diabetic macular edema (DME), and Drusen, that are
the most common pathologies resulting in acquired blindness [1]. With
approximately 30 million invasive OCT scans being performed each year
worldwide, there is a need to identify patients with OCT images that need
prioritized attention. The existing work in [2] shows 92-99% classification
accuracy for OCT image classification using large annotated OCT image data
sets to train a CapsuleNet model. Contrarily, in this work we propose a novel
active learning framework [4] [5] to significantly reduce the training
complexity by reducing the size of annotated training data set. The proposed
system and steps are explained in the Fig. 1.

Figure 1: Overview of the proposed workflow. a) A unsupervised model SimCLR is
trained on the constrastive learning objective to learn useful representation
of each image. b) Dimension reduction is applied followed by unsupervised
clustering and identification of cluster peripheral samples. c) Label
Propagation is applied using a pre-labelled annotated image set. Images with
the most uncertain labels are selected as the optimal training sample subset.
d) The Resnet pretrained with SimCLR is now fine-tuned with the annotated
images selected in c) for classification.

This paper makes three key contributions. First, we present a novel semi-
supervised framework to identify the most uncertain set of images that are not
straightforward to classify or annotate. The proposed method can identify upto
2% of training images samples that are encoded using a self-supervised method
[3] and that are found to lie along classification decision boundaries. The
proposed active learning algorithm enables prioritization of patients with
such uncertain images to be seen early for treatment. Second, we present a
two-step classifer training method where feature extractor layers are first
pre-trained using image augmentations without labels and SimCLR [3] framework.
Next, the classifier is fine-tuning using the labelled uncertain image set for
training. This process achieves 96.45% classification test accuracy for OCT-
based pathologies by training on atmost 2000 images only. Third, we utilize
the Gradcam library to analyze the regions of interest (ROIs) [1] that explain
underlying pathology. We apply random image occlusion followed by feature
analysis to identify primary and secondary ROIs for pathology to lie between
the external limiting membrane and choroidal layers.

##  II Data and Methods

The data set and methods used in this work are described below.

###  II-A Data: OCT Image Dataset

In this work we apply the OCT data set in [1], where, training data set
contains images annotated for 4 classes {CNV, DME, Drusen, Normal} with
{26318, 8616, 11350, 37205}, images per class, respectively. The test set has
968 images with the same 4 class labels.

###  II-B Methods and Mathematical Frameworks

The first step to medical image classification is feature extraction while
preserving the structural, contextual and textural features. The methods used
to encode the OCT images to a latent vector space, and the proposed semi-
supervised image sampling methods are described below 111Code is available at
https://github.com/anoopsanka/retinal_oct.

####  II-B1 Image Encoding by Self-supervised Learning

SimCLR [3] is a recent model proposed for unsupervised visual representation
learning. A minibatch of $n$ images is sampled and applied with pairs of
augmentation to produce $2n$ images. These augmented pairs are fed into a
feature encoder (Resnet model) followed by a projection layer (Multilinear
Perceptron) to obtain an latent vector ${\bf z}$. The training objective is to
maximize the agreement of the latent vectors from the same image with
augmented views (positive pair), while repulsing from other different images
($2n-1$ negatives). In SimCLR, the similarity between two views is defined
using cosine similarity, i.e. $\rm{sim({\bf u,v})}={\bf u}^{T}{\bf v}/||{\bf
u}||~{}||{\bf v}||$. Thus, the loss function for each positive pair of sample
$(a,b)$ is defined to be

$\mathcal{L}_{a,b}=-\rm{log}\frac{exp(\rm{sim({\bf
z_{a},z_{b}})}/\tau)}{\sum_{\begin{subarray}{c}c=1\\\ c\neq
a\end{subarray}}^{2n}exp(\rm{sim({\bf z_{a},z_{c}})/\tau})},$ (1)

where $\tau$ refers to the temperature parameter.

The SimCLR model is trained with batch size of 128 for 120 epochs. A batch
size of 128 provides 255 negative samples per positive pair from two augmented
views. We apply the LAMB optimizer since training with a large batch size may
lead to instability using SGD with momentum. Similar to the original SimCLR
implementation in [3], we use a linear warmup for the first 10 epochs and
decay the learning rate with the cosine decay schedule. Here, the default
augmentation strategy from [3] is adopted. Sample augmentations are shown in
Figure 2. The model training is stopped at epoch 75 with a contrastive
accuracy of 99.77%.

Figure 2: Row 1 shows the original image with different label class. Rows 2, 3
show the same images in Row 1 with SimCLR default augmentations.

Each OCT image $I$ is resized to [224x224] followed by SimCLR self-supervised
latent vector encoding process. The output of the Resnet encoder for each
image of size [7x7x512] is subjected to global averaging to extract $J$, which
is a [1x512] dimensional encoded vector representation for each image. We
further reduce the dimensions of the latent vector representation by
subjecting $J$ to PCA-based decomposition, thereby retaining 8 top principal
components per image that are found to be most representative. Thus, the data
set used for subsequent selective sub-sampling is represented by
{$X\in\mathcal{R}^{[n\times d]},Y\in\mathcal{R}^{[n\times 1]}$}, where $X$
depicts sample features with $d=8$ dimensions for $n=83,484$ samples, and $Y$
represents the pathological class labels. Further, cluster-specific sub-
sampling methods described below.

####  II-B2 Semi-supervised Sub-sampling

With our intention to train a classification model with minimal training
labels, we use k-means to cluster $X$ into four categories/clusters with means
and standard deviations, respectively, as $[\mu_{k},\Sigma_{k}]\forall
k=[1:4]$. Next, we identify uncertain samples as the ones that lie farther
away from the cluster centers and are therefore more likely to lie along the
classifier decision boundaries. In (2), we compute the probability for a
sample $x$ to belong to cluster $k$.

$\displaystyle
p_{k}(x)=\frac{exp(-\frac{1}{2}(x-\mu_{k})\Sigma_{k}^{-1}(x-\mu_{k})^{T})}{\sqrt{(2\pi)^{d}|\Sigma|}}$
(2)

Samples that have normalized probabilities of belonging to any cluster
distribution in the uncertain range of [0.4-0.6] can be considered to lie at
cluster peripheries. These samples are collected as $(X_{S}\subset
X),X_{S}\in\mathcal{R}^{[mxd]}$, where $m<n$. Other methods such as relative
z-scores and Silhouette Scores per sample were also evaluated to identify
uncertain samples, but the Gaussian distribution-based sample isolation proved
most effective in identifying the corner cases per cluster. Next, we use a
small set of labelled data ($L=\\{X_{l},Y_{l}\\}$) with atmost $n_{l}=80$
labelled samples as seed to apply label propagation [6] to all unlabelled
samples in $X_{S}$. The goal here is identification of samples that
demonstrate high variances for transductive labelling. Thus, uncertain samples
that acquire different class labels at different sample runs can be indicative
of the uncertain space around decision boundaries. Using the semi-supervised
label propagation method shown in Algorithm 1, we identify the most uncertain
minimum sub-sampled data set $(X_{T}\subset
X_{S}),X_{T}\in\mathcal{R}^{[r\times d]}$, $r<<n$, that needs to be labelled
to adequately train a deep learning classifier [7].

Output: Most uncertain sample set $X_{T}$

Input: $X_{S}$, $L=\\{X_{l},Y_{l}\\}$

for  _$e=$ 0 to $10$_ do

Randomly select 5 samples per class from $L$;

for  _$x_{s}\in X_{S}$_ do

Label-spreading with $\alpha=0.01$, kernel$=rbf$;

end for

return _labels $Y_{S}$ for $X_{S}$_;

U[:,e] $\longleftarrow Y_{S}$;

end for

$Given$: $U\in\mathcal{R}^{[m\times 10]},X_{S}\in\mathcal{R}^{[m\times
d]},X_{T}=[.]$;

for  _$j=$ 0 to m_ do

$l$=unique(U[j,:]);

if  _length $(l)>=5$_ then

append $X_{S}(j,:)$ to $X_{T}$;

end if

end for

Algorithm 1 Semi-supervised Sub-sampling

In Algorithm 1, label-spreading using rbf-kernel is applied with $\alpha=0.01$
implying that once a label is acquired by a sample in $X_{S}$ it is less
likely to be modified again. Thus, we identify the samples that acquire highly
variable labels across 10 label spread runs. The minimal sub-sampled data set
$X_{T}$ is returned for labelling and classifier training thereafter.

####  II-B3 Pathology Classification and Explainability

Once the labels for the minimal sub-sampled data set are collected, we use
this minimal training data to further fine-tune the SimCLR pre-trained Resnet
classifier. Next we apply the Gradcam library and tf_explain libraries in
Python [8] to identify primary and secondary features corresponding to ROIs
that explain the underlying pathology as shown in Fig. 3. Here, we observe
that the feature extraction layers of Resnet focus on the regions between the
inner segment layer and retinal pigment epithelium layer in the OCT images to
classify pathological vs. normal images. In case this ROI is occluded, the
secondary regions that are analyzed for each OCT image include the choroidal
regions for macular OCT images.

Figure 3: Examples of primary and secondary ROIs identified on test data.

##  III Experiments and Results

To enable minimal data sub-sampling followed by pathology classification, we
perform three major experiments. First, we analyze the repeatability of the
proposed semi-supervised sub-sampling method by analyzing across multiple runs
and assessing the numbers of images detected per class. Second, we compare the
classification performance of the proposed sub-sampling with random stratified
sampling for Resnet classification. Third, we analyze the robustness of
features explained by the sub-sampled uncertain images.

###  III-A Analysis of Sub-sampling Repeatability

We repeat Algorithm 1 for 20 runs and record the numbers of images
corresponding to each class category detected per run. The average and
standard deviations in the number of samples per class are shown in Table I.

TABLE I: Analysis of sub-sampled images across 20 runs.

Metric | Normal | Drusen | DME | CNV | Total
---|---|---|---|---|---
Average | 711 | 274 | 261 | 466 | 1707
Std. dev. | 78 | 50 | 26 | 141 | 251

Here, we observe significant down sampling and relative consistency across the
number of images sub-sampled across multiple runs when compared to the
original sample size of 83,484 samples. Thus, the proposed method selects less
than 2000 training samples, that represents about 2% of the total sample size,
to be annotated for subsequent classification.

###  III-B Classification of sub-sampled Images

To analyze the importance of the novel semi-supervised sub-sampling method
proposed in this work, we further fine-tune the Resnet model that was
previously trained using the contrastive loss for SimCLR. For this fine-
tuning, we use several sets of sub-sampled labelled data to asses the impact
of transfer learning and also to identify the minimum number of samples
required to train an explainable pathology classifier.

For this experiment, a set of labelled data ($L$) with $n_{l}$ samples is
subjected to stratified random 80/20 train/validation split. This data set is
then used to train the Resnet model layers with categorical crossentropy loss,
balanced class weights and Adam Optimizer with a learning rate $10^{-5}$ .
Early stopping is applied when the validation loss has not improved for more
than 3 epochs. To analyze the importance of data sub-sampling over using the
complete training data set, we apply stratified random sampling to isolate 1,
5, 10% of samples per class to create labelled data $L$ that is then used to
fine-tune the Resnet model. The classification performances using stratified
random sub-sampling in terms of average classification accuracy, precision,
recall and F1 score [1] are shown in Table II.

TABLE II: Averaged classification performances from Stratified and Proposed
Semi-supervised Sampling across several runs.

Stratified Sample Size (%) | Accuracy | Precision | Recall | F1
---|---|---|---|---
1% | 0.9101 | 0.9221 | 0.9101 | 0.9083
5% | 0.9834 | 0.9844 | 0.9834 | 0.9834
10% | 0.9880 | 0.9883 | 0.9880 | 0.9880
100% | 0.9989 | 0.9989 | 0.9989 | 0.9989
Proposed Semi-supervised Sampling
2% | 0.9645 | 0.9675 | 0.9625 | 0.9625

Here, we observe that training by atmost 5% random stratified samples, the
Resnet model is significantly well trained to visually explain pathology
classification. Improvements in classification performances thereafter are
incremental and often not generalizable. Also, in Fig. 4 we observe the UMAP
representations of classified samples with 1% and 98% of the data samples. We
observe significant cluster separability as the training data set increases.

(a) UMAP with 1% data

(b) UMAP with 98% data

Figure 4: UMAP representations of classified samples trained with $L$
corresponding to 1% and 98% stratified random samples, respectively.
Separability across clusters is significantly higher with 98% of training
samples.

Our goal is to identify the minimal training data set to achieve significant
cluster separability. Thus, we analyze the performance of the proposed sub-
sampling method to train the classifier in Table II. We observe that the
proposed sub-sampled data set that isolates 2% training images around decision
boundaries achieves about 96.45% classification accuracy. This performance is
comparable to InceptionV3 based model in [2] that achieved 96.1% accuracy with
complete training data set. Thus, the proposed sub-sampling method intuitively
lies within the 5% sampling size as learned from the random stratified
sampling. Additionally, the Resnet model fine-tuned on the proposed minimal
sub-sampled data achieves consistent precision, and recall performances, which
indicates consistencies in the numbers of classified false positives and false
negatives.

The confusion matrix for Resnet classifier after being trained with 1%
stratified random samples and with the 2% proposed sub-sampling method are
shown in Fig. 5.

(a) 1% Stratified Random sampling.

(b) 2% Proposed Sub-sampling.

Figure 5: Confusion Matrix for Classification with varying training samples.

Here, we observe that using the proposed sub-sampled data, some DME images get
mis-classified as CNV where small cystic regions develop close to the inner
sub-retinal layers. Also some images with thin drusen membranes are mistaken
for normal images and some images with CNV and poor contrast are classified as
DME. Most of these classifications can be improved by pre-processed zooming in
and contrast corrections applied on test images in future works.

###  III-C Explainability of Pathology

Finally, we incorporate the Gradcam library to analyze the ROI that
contributes to pathology classification. In Fig. 6 (a), we observe some sub-
sampled images using Algorithm 1 that lie along the border of two distinct
clusters, thereby enhancing classification performance for samples of those
clusters. Further, Fig. 6 (b) shows the spatial robustness of the fine-tuned
Resnet for pathology explainability.

(a) Examples of proposed sub-sampled images along cluster peripheries.

(b) ROI explanations for SimCLR augmentations shows the same ROI is
highlighted per augmentation.

Figure 6: Examples of pathology explanation on (a) sub-sampled training
images, (b) test images, respectively.

##  IV Conclusion

In this work we present a novel semi-supervised algorithm that uses encoded
latent features per image to identify a minimal sub-sampled data set that can
be useful to train and visually explain pathology. We incorporate a Resnet
model for pathology classification that is first trained using contrastive
loss to distinguish an original and augmented sample from other images. Next,
we perform label spreading on the encoded latent space to identify a training
set of upto 2% of samples that represent the most uncertainty in feature
space. Once trained on this labelled minimal saub-sampled data, the final
classifier can achieve upto 97% classification accuracy and is capable of
visually explaining pathological sites. Future work can be directed towards
extending the proposed image sub-sampling method and Resnet training to other
pathologies.

## Acknowledgement

We gratefully acknowledge the direction setting and help by the AI team at
Samsung SDSA.

## References

  * [1] D. S. Kermany, M. Goldbaum, W. Cai, C. C. Valentim, H. Liang, S. L. Baxter, A. McKeown, G. Yang, X. Wu, F. Yan _et al._ , “Identifying medical diagnoses and treatable diseases by image-based deep learning,” _Cell_ , vol. 172, no. 5, pp. 1122–1131, 2018.
  * [2] T. Tsuji, Y. Hirose, K. Fujimori, T. Hirose, A. Oyama, Y. Saikawa, T. Mimura, K. Shiraishi, T. Kobayashi, A. Mizota _et al._ , “Classification of optical coherence tomography images using a capsule network,” _BMC ophthalmology_ , vol. 20, no. 1, pp. 1–9, 2020.
  * [3] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” 2020.
  * [4] S. Dasgupta and D. Hsu, “Hierarchical sampling for active learning,” in _Proceedings of the 25th international conference on Machine learning_ , 2008, pp. 208–215.
  * [5] H. Hao, S. Didari, J. O. Woo, H. Moon, and P. Bangert, “Highly efficient representation and active learning framework for imbalanced data and its application to covid-19 x-ray classification,” _https://arxiv.org/abs/2103.05109_ , 2021.
  * [6] K. Sun, Z. Min, and J. Wang, “Pp-pll: Probability propagation for partial label learning,” in _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_. Springer, 2019, pp. 123–137.
  * [7] C. Szegedy, S. Ioffe, V. Vanhoucke, and A. Alemi, “Inception-v4, inception-resnet and the impact of residual connections on learning,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 31, no. 1, 2017.
  * [8] R. Meudec. (Accessed, Jan, 2021) tf_explain. [Online]. Available: https://github.com/sicara/tf-explain


# Quantifying the role of antiferromagnetic fluctuations in the
superconductivity of the doped Hubbard model

Xinyang Dong  Department of Physics, University of Michigan  Emanuel Gull
Department of Physics, University of Michigan  A. J. Millis  Center for
Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York,
NY 10010  Department of Physics, Columbia University, 538 West 120th Street,
New York, New York 10027

Superconductivity arises from the pairing of charge $e$ electrons into charge
$2e$ bosons (“Cooper pairs”) and their condensation into a coherent quantum
state. In conventional superconductors such as lead, a comparison of the
frequency dependence of the superconducting gap function to the frequency
spectrum of the phonons (quantized lattice vibrations) [1, 2] establishes that
the electron-phonon interaction provides the “pairing glue” that binds
electrons into Cooper pairs. Many “unconventional” superconductors are now
known [3, 4, 5, 6] in which the pairing glue is believed not to be provided by
phonons. Substantial indirect evidence indicates that in many cases the
relevant interaction is the exchange of spin fluctuations [7, 8, 9, 10, 11]
but direct evidence has been lacking and many other mechanisms have been
proposed [12, 13, 14, 15, 16, 17]. This paper quantifies the strength of the
electron-spin fluctuation coupling in the two dimensional Hubbard model, a
minimal theoretical model of materials with strong electron-electron
interactions, by analysing the frequency dependence of the computationally
determined electron self energy, superconducting gap function and spin
fluctuation spectrum. The analysis shows that about half of the
superconductivity is attributable to RPA-like weak coupling spin fluctuations,
with the other half coming from higher energy processes.

Unconventional superconductivity occurs in quantum materials exhibiting
unusual electronic properties believed to arise from strong electron-electron
interactions. Theoretical analysis requires a model that captures the
essentials of the correlated electron physics and can be studied non-
perturbatively. The Hubbard model has been proposed as the minimal theoretical
model of quantum materials such as the copper-oxide based high-Tc
superconductors [18]. This model describes electrons hopping among sites of a
lattice (here we consider the two dimensional square lattice case with
nearest-neighbor hopping of amplitude $t$) and subject to a site-local
repulsive interaction $U$ (here we focus on the intermediate-strong coupling
value $U=6t$). Cluster dynamical mean field theory (DCA) [19] provides non-
perturbative access to both the static phase diagram and dynamical properties.
Importantly, the electron propagator and spin fluctuation spectrum are
computed within the same formalism and at the same level of approximation,
enabling a quantitative analysis of the electron-spin fluctuation interaction.
The resulting solution [20, 21] produces a good qualitative description of the
physics of the high-Tc copper oxide superconductors, including a high-doping
Fermi liquid regime, a Mott insulator, a low doping pseudogap and an
intermediate-doping dome of d-wave superconductivity. The extent to which a
stripe magnetic phase preempts the superconducting phase found in the DCA is
currently under debate [22, 23] but we emphasize that the superconductivity
found in DCA is well defined and locally stable, with properties that we now
study.

The DCA [19] proceeds by tiling the Brillouin zone into $N$ equal-area non-
overlapping tiles $a$ centered at momentum points $K_{a}$ and approximating
the electron self energy as
$\Sigma(k,\omega)\rightarrow\Sigma_{K_{a}}(\omega)$ for $k$ in tile $a$, so
that the momentum dependence is approximated as a sum of piecewise constant
functions and the full frequency dependence is retained. The
$\Sigma_{K_{a}}(\omega)$ are obtained from the solution of a $N$-site quantum
impurity model with the same interaction $U$ as in the original model and
single particle parameters obtained by a self-consistency condition. We have
chosen $N=8$ which provides sufficient momentum resolution while allowing for
calculation of the detailed dynamical information needed here.

We used continuous-time quantum Monte Carlo methods [24, 25] to solve the
impurity model for a range of temperatures and carrier concentrations. We
present here results obtained for doping $x=0.1$ (carrier concentration
$n=1-x$ per site) and temperatures as low as $T=t/50$. For this carrier
concentration at $U=6t$ the normal state is a momentum-space differentiated
Fermi liquid outside the pseudogap regime corresponding to the over-doped side
of the cuprates and the superconducting state, which we explicitly construct,
appears below a transition temperature $T_{c}\approx t/40$. We calculate the
normal (N) and anomalous (A) components of the electron self energy. Using
recent algorithmic developments [26] we also calculated the impurity model
spin susceptibilities $\chi$ in both normal and superconducting states.

Spin fluctuation theories yield the spin fluctuation (SF) contribution to the
normal (N) and anomalous (A) self energies in terms of the spin susceptibility
and normal and anomalous components of the Green function $G_{K_{a}}^{N/A}$ as
[27, 28, 8, 7] (see Fig. 1)

$\Sigma_{K_{a}}^{SF;N/A}(\omega)=g^{2}\frac{1}{\beta
N}\sum_{Q,\Omega}\chi_{Q}(\Omega)G_{K_{a}-Q}^{N/A}(\omega-\Omega).$ (1)

We assess the relevance of spin fluctuations by using our calculated $G$ and
$\chi$, along with an estimated coupling constant $g$ to compare $\Sigma^{SF}$
(Eq. 1) to our numerically calculated self energies.

Figure 1:  Spin-fluctuation diagrams for normal and anomalous self energy.
Solid lines: Normal or anomalous Green’s function; Wavy lines: Magnetic
susceptibility.  Figure 2: Upper panel: Imaginary part of normal component of
the Matsubara self energy for the $M$-point at $U=6$, $\beta=35$, and
$\mu=-1.0$ compared to spin fluctuation self energy computed with $g^{2}=3.8$.
Inset: $\text{Im}\Sigma^{N}(i\omega_{n})/\omega_{n}$. Lower panel: Negative of
the analytically continued real-axis M-point self-energy and spin-fluctuation
contribution computed with $g^{2}=3.8$. Inset: Self energy over a wide
frequency range.

To estimate $g^{2}$ we partition the imaginary part of the real frequency self
energy into a low frequency part that we suppose arises mainly from spin
fluctuations and a higher frequency part that represents all of the other
processes contributing to the imaginary part of the self energy:
$\text{Im}\Sigma^{N}_{K}(\omega)=\text{Im}\Sigma^{SF;N}_{K}(\omega)+\text{Im}\Sigma^{\text{high;N}}_{K}(\omega)$.
We take $\Sigma^{SF;N}$ to have the functional form of Eq. 1 and determine
$g^{2}$ by requiring consistency with our numerically computed self energies.
We have computed the self energies in all momentum tiles but focus here on the
self energies corresponding to the tiles centered on the points
$M=(\pi,0)/(0,\pi)$ where the superconducting gap is maximal and the normal
component of the self energy is largest. We consider consistency both directly
on the Matsubara axis (avoiding the ambiguities associated with analytical
continuation) and on the real axis. For the imaginary axis analysis we note
that the quantity
$Z^{N}=\frac{\partial[\text{Re}\Sigma^{N}(\omega)]}{\partial\omega}$ related
to the normal state mass enhancement may be estimated from Matsubara axis
results as
$Z^{N}=\frac{\text{Im}(\Sigma^{N}(i\omega_{1})-\Sigma^{N}(i\omega_{0}))}{\omega_{1}-\omega_{0}}$
(see inset of Fig. 2) and cannot be larger than the contribution from the spin
fluctuation sector. We therefore write $Z^{N}\approx
Z^{\text{high};N}+Z^{SF;N}$. $Z^{\text{high};N}$ is computed from a minimal
two-parameter fit of
$\text{Im}\Sigma^{\text{high};N}(i\omega_{n})=-\frac{A}{\pi}\frac{\omega_{n}}{\omega_{n}^{2}+x_{0}^{2}}$,
where $x_{0}$ and $A$ are determined by the maximum of
$|\text{Im}(\Sigma^{N}-\Sigma^{SF;N})|$. We performed the same analysis using
the analytically continued numerical data and the real axis version of Eq. 1,
bounding $g^{2}$ from the argument that at low frequencies the imaginary part
of the self energy should not be larger than spin fluctuation term.

The upper panel of Fig. 2 shows the Matsubara analysis of the M-point self
energy and the lower panel shows the real axis fits. Both cases are consistent
with a value of $g^{2}=3.80$ implying that about $2/3$ of $Z^{N}$ comes from
the electron-spin fluctuation interaction.

With the spin fluctuation spectrum and the electron-spin fluctuation coupling
constant in hand we next determine the extent to which superconductivity
arises from spin fluctuations by solving the anomalous component of Eq. 1 and
comparing the result to the CTQMC solution. We begin with the equation for the
transition temperature $T_{c}$, obtained by linearizing Eq. 1 in the anomalous
component of the self energy. The resulting equation is a linear eigenvalue
equation for eigenvector $\Sigma^{A}(i\omega_{n})$; the largest eigenvalue
increases as temperature decreases and $T_{c}$ is the temperature at which the
leading eigenvalue equals unity. Using our estimated $g^{2}=3.8$ we find that
at temperature $T=t/40$ the leading eigenvalue is about $0.5$, so that
increasing the net pairing strength by a factor of about two would be needed
to bring the leading eigenvalue up to $1$ (in fact a larger increase would be
required because the coupling constant of the normal state self energy means
the transition temperature does not vary linearly with the coupling).

Figure. 3 compares the QMC anomalous self energy to the spin fluctuation self
energy $\Sigma^{SF;A}$. We note that the spin fluctuation interaction has two
components, one from fluctuations near the antiferromagnetic wavevector
$(\pi,\pi)$ and one from fluctuations at small momenta near $Q=(0,0)$. The
small momentum fluctuations make a negative contribution to $\Sigma^{A}$. At
the lowest Matsubara frequency the $\Sigma^{SF;A}$ produced by the spin
fluctuation theory is approximately half of the QMC self energy, again
indicating that spin fluctuation theory [1] alone cannot account for the
superconductivity.

Figure 3:  Total measured anomalous self-energy $\Sigma^{A}$ and estimated
spin-fluctuation contribution $\Sigma^{SF;A}$ at $U=6$, $\beta=50$ and
$\mu=-1.0$ ($n=0.90$). Also shown are the individual contributions to
$\Sigma^{SF;A}$ from transferred momenta $Q=(\pi,\pi)$ and $Q=(0,0)$.

We now examine in Fig. 4 the frequency dependence of the gap function
$\Delta(\omega)$, a complex function of real frequency defined in terms of the
normal and anomalous self energies as [29, 30]

$\displaystyle\Delta(\omega)=\Sigma^{A}(\omega)\Big{/}\left[1-\frac{\Sigma^{N}(\omega)-\Sigma^{N}(-\omega)}{2\omega}\right]\,.$
(2)

Following Ref. [1] we compare the frequency dependence of the spin fluctuation
spectrum, the imaginary part of the DCA-computed gap function and the
estimated gap function computed by solving Eq. 1 using the CT-QMC-computed
$\Sigma^{A}$ and $\chi$. The real frequency quantities are obtained from
maximum entropy analytical continuation of imaginary frequency data obtained
at $T=t/50$, well below the superconducting transition temperature. As noted
in Ref. [1] the presence of a gap in the electron Greens function means that a
peak in $\chi$ at a frequency $\omega_{peak}$ implies a peak in $\Delta$ at
$\omega+\omega_{peak}$ so we shift $\chi$ by the zero frequency gap function
in the comparison.

We emphasize that the uncertainties in the analytical continuation are not
small; while areas are reliably estimated, peak heights and widths are subject
to some uncertainty. We see from Fig. 4 that while the peaks in the gap
function and shifted $\chi$ roughly coincide, the spin fluctuation
contribution to the imaginary part of the gap function is concentrated at low
frequencies, decaying much more rapidly than the DCA-computed gap function,
further demonstrating the importance of a high-frequency non-spin-fluctuation
contribution to the electron self energy.

Figure 4:  Comparison of the true gap function, gap function from spin
fluctuation, the AFM susceptibility
$\text{Im}\chi^{\omega+\Delta_{0}}_{m,(\pi,\pi)}$ and the FM susceptibility
$\text{Im}\chi^{\omega+\Delta_{0}}_{m,(0,0)}$ shifted by the
$\Delta_{0}=\text{Re}\Delta(\omega=0)=0.057$ at $U=6$, $\beta=50$ and
$\mu=-1.0$ ($n=0.90$). Upper panel: Imaginary part. Inset: Integral of
$\text{Im}\Delta(\omega)/\omega$ starting from $\omega=0$. Lower panel: Real
part.

Spin fluctuation theories are widely considered to be promising candidates for
theories of superconductivity. Here we have performed a quantitative study, in
a well defined, numerically controlled theory, of the extent to which this is
actually the case. We take advantage of recent developments [26] that permit
us to obtain numerically exact results for the spin susceptibility of the 8
site DCA approximation to the Hubbard model in both normal and superconducting
states. Access to the spin fluctuation spectrum enables us to compare the spin
fluctuation theory calculation of the normal state self energies to
numerically exact results for the same quantities, thereby allowing an
estimate of the electron-spin fluctuation coupling constant. Knowledge of the
coupling constant then enables a quantitative analysis of the contribution of
spin fluctuations to the superconducting transition temperature and to the
magnitude and form of the superconducting gap function. In qualitative
consistency with previous results [10] we find that low-frequency spin
fluctuations of the RPA type contribute to the superconductivity, but we find
that quantitatively only about half of the pairing can be attributed to these
fluctuations. The other half of the pairing arises from higher frequency
fluctuations, whose nature and precise physical origin remains to be
determined.

The cluster dynamical mean field theory used in this work has limitations. It
does not adequately capture for example the stripe physics that may preempt
superconductivity in some parameter ranges, and the cluster sizes available,
while large enough to provide results that compare well to experiment and more
exact calculations, cannot capture many of the interesting specifics of
superconducting phenomenology. However, it is important to emphasize that the
method provides a single internally consistent computational scheme that
produces a well defined locally stable superconducting phase whose properties
can be studied, and that provides, at the same level of approximation, normal
and anomalous self energies and spin fluctuation spectra, enabling a
theoretically meaningful comparison.

Our finding that the RPA-type spin fluctuations examined are not the dominant
form of superconductivity suggests more generally that spin fluctuation
theories of this type may miss important aspects of correlated electron
superconductivity. Our finding also suggests that if the nature of the higher
frequency contributions to the pairing could be elucidated, tuning these
degrees of freedom might be an effective strategy for raising the transition
temperature.

Acknowledgements: EG and XD are supported by NSF DMR 2001465. The Flatiron
Institute is a division of the Simons Foundation.

## References

  * Scalapino  _et al._ [1966] D. J. Scalapino, J. R. Schrieffer, and J. W. Wilkins, Phys. Rev. 148, 263 (1966).
  * McMillan [1968] W. L. McMillan, Phys. Rev. 167, 331 (1968).
  * Steglich  _et al._ [1979] F. Steglich, J. Aarts, C. D. Bredl, W. Lieke, D. Meschede, W. Franz, and H. Schäfer, Phys. Rev. Lett. 43, 1892 (1979).
  * Bednorz and Müller [1986] J. G. Bednorz and K. A. Müller, Zeitschrift für Physik B Condensed Matter 64, 189 (1986).
  * Maeno and Hashimoto [1994] Y. Maeno and H. e. a. Hashimoto, Nature 372, 532–534 (1994).
  * Kamihara  _et al._ [2006] Y. Kamihara, H. Hiramatsu, M. Hirano, R. Kawamura, H. Yanagi, T. Kamiya, and H. Hosono, Journal of the American Chemical Society 128, 10012 (2006).
  * Miyake  _et al._ [1986] K. Miyake, S. Schmitt-Rink, and C. M. Varma, Phys. Rev. B 34, 6554 (1986).
  * Scalapino  _et al._ [1986] D. J. Scalapino, E. Loh, and J. E. Hirsch, Phys. Rev. B 34, 8190 (1986).
  * Scalapino [1999] D. J. Scalapino, Journal of Low Temperature Physics 117, 179 (1999).
  * Maier  _et al._ [2008] T. A. Maier, D. Poilblanc, and D. J. Scalapino, Phys. Rev. Lett. 100, 237001 (2008).
  * Maier and Scalapino [2019] T. A. Maier and D. J. Scalapino, npj Quantum Materials 4, 30 (2019).
  * Castellani  _et al._ [1996] C. Castellani, C. Di Castro, and M. Grilli, Zeitschrift für Physik B Condensed Matter 103, 137 (1996).
  * Varma [1997] C. M. Varma, Phys. Rev. B 55, 14554 (1997).
  * Capone  _et al._ [2004] M. Capone, M. Fabrizio, C. Castellani, and E. Tosatti, Phys. Rev. Lett. 93, 047001 (2004).
  * Anderson [2007] P. W. Anderson, Science 316, 1705 (2007).
  * Stanescu  _et al._ [2008] T. D. Stanescu, V. Galitski, and S. Das Sarma, Phys. Rev. B 78, 195114 (2008).
  * Saito  _et al._ [2015] T. Saito, Y. Yamakawa, S. Onari, and H. Kontani, Phys. Rev. B 92, 134522 (2015).
  * Anderson [1987] P. W. Anderson, Science 235, 1196 (1987).
  * Maier  _et al._ [2005] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev. Mod. Phys. 77, 1027 (2005).
  * Gull  _et al._ [2010] E. Gull, M. Ferrero, O. Parcollet, A. Georges, and A. J. Millis, Phys. Rev. B 82, 155101 (2010).
  * Gull  _et al._ [2013] E. Gull, O. Parcollet, and A. J. Millis, Phys. Rev. Lett. 110, 216405 (2013).
  * Zheng  _et al._ [2017] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Science 358, 1155 (2017).
  * Qin  _et al._ [2020] M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 10, 031016 (2020).
  * Gull  _et al._ [2008] E. Gull, P. Werner, O. Parcollet, and M. Troyer, EPL (Europhysics Letters) 82, 57003 (2008).
  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
  * Chen  _et al._ [2015] X. Chen, J. P. F. LeBlanc, and E. Gull, Phys. Rev. Lett. 115, 116402 (2015).
  * Layzer and Fay [1971] A. Layzer and D. Fay, Int. J. Magn. 1, 135 (1971).
  * Anderson and Brinkman [1973] P. W. Anderson and W. F. Brinkman, Phys. Rev. Lett. 30, 1108 (1973).
  * Poilblanc and Scalapino [2002] D. Poilblanc and D. J. Scalapino, Phys. Rev. B 66, 052513 (2002).
  * Gull and Millis [2014] E. Gull and A. J. Millis, Phys. Rev. B 90, 041110 (2014).