# Efficient Batch Dynamic Graphlet Counting

Hriday G<EMAIL_ADDRESS>BITS Pilani, Hyderabad
CampusHyderabadIndia ,  Pranav Saikiran Sista<EMAIL_ADDRESS>pilani.ac.in BITS Pilani, Hyderabad CampusHyderabadIndia and  Apurba Das
<EMAIL_ADDRESS>BITS Pilani, Hyderabad CampusHyderabadIndia

###### Abstract.

Graphlet counting is an important problem as it has numerous applications in
several fields, including social network analysis, biological network
analysis, transaction network analysis, etc. Most of the practical networks
are dynamic. A graphlet is a subgraph with a fixed number of vertices and can
be induced or non-induced. There are several works for counting graphlets in a
static network where graph topology never changes. Surprisingly, there have
been no scalable and practical algorithms for maintaining all fixed-sized
graphlets in a dynamic network where the graph topology changes over time. We
are the first to propose an efficient algorithm for maintaining graphlets in a
fully dynamic network. Our algorithm is efficient because (1) we consider only
the region of changes in the graph for updating the graphlet count, and (2) we
use an efficient algorithm for counting graphlets in the region of change. We
show by experimental evaluation that our technique is more than 10x faster
than the baseline approach.

##  1\. Introduction

Graphlets are induced subgraphs with significant frequency counts that can be
found in any network. The graphlets which are considered to be important are
those of sizes 3, 4, and 5 and can be identified in both directed and
undirected graphs. Given an undirected graph $G$, with $n$ vertices and $m$
edges, counting graphlets of size $k$ is defined as counting the number of
subgraphs isomorphic to subgraphs of a fixed size $k$. It is a computationally
intensive problem since there are $O(|V|^{k})$ possible subgraphs of size $k$.
Despite that, graphlet counting is an important problem due and it has served
as a building block in solving problems in several areas.

In computational biology, they have been used for detecting cancer cells
(Milenković et al., 2010), and for analyzing protein-protein interactions
(Milenković and Pržulj, 2008; Pržulj et al., 2004). Other areas where they
have been used include bioinformatics (Vacic et al., 2010; Wang et al., 2012;
Saha and Hasan, 2015; Wang et al., 2017), computer vision (Harchaoui and Bach,
2007; Zhang et al., 2012, 2013), social network analysis (Seshadhri et al.,
2013; Rahman et al., 2014; Wang et al., 2014; Janssen et al., 2012) and so on.
They have also been used to compare graphs using graph kernels (Shervashidze
et al., 2009).

Although most research involving graphlets has been done on static graphs,
progress has been made on dynamic networks. This is pretty significant since
this will open up ways to analyze networks that change over time. Solving such
problems will give graphlets more real-world applications since most problems
are dynamic in nature. This is especially true in the case of biological
networks. Knowing the counts of graphlets in a biological network can provide
information about the network’s properties. However, these networks can
frequently change over time. For example, when a stem cell differentiate into
other cell types during its development process (Mukherjee et al., 2018), or
when chromosomes’ chromatin structures change due to various events (Mukherjee
et al., 2018). Due to the ever-changing nature of these networks, maintaining
the counts of graphlets statically will not be sufficient as they will be
outdated. It is necessary to update the counts whenever the network changes.
Similar problems are present in other areas as well. For instance, social
networks are constantly changing with new connections being added and existing
connections are becoming obsolete, or new transactions are getting added and
old transactions are getting removed from the historical data store in
transaction databases.

In order to solve such problems, we have proposed an algorithm that maintains
the counts of 4-node graphlets in a network which is changing over time due to
the addition and deletion of edges.

##  2\. Related Works

Most of the previous works have focused on counting only a special type of
graphlets, such as cliques and cycles (Danisch et al., 2018; Bera and
Seshadhri, 2020; McGregor and Vorotnikova, 2020; Shin et al., 2020). However,
counting all graphlets of a fixed size has been used extensively in recent
times in graph mining and machine learning applications. We discuss the state-
of-the-art algorithms for graphlet counting in static and dynamic graphs.

Graphlet counting in static network: Counting graphlets exactly (Ahmed et al.,
2015a; Pinar et al., 2017) and approximately (Rossi et al., 2018; Rahman et
al., 2014; Wang et al., 2017; Jha et al., 2015) in a static network has been
studied extensively. In (Ahmed et al., 2015a), the authors develop an
algorithm for counting graphlets up to size $4$. In their work, they iterate
over edges to compute smaller-size graphlets and use combinatorial arguments
to extend them to large graphlets in constant time. Their algorithm is
scalable and easy to parallelize. In (Pinar et al., 2017), the authors have
developed an algorithm for exactly counting all graphlets of size $5$. Their
framework is based upon cutting a big subgraph pattern into smaller ones and
using the smaller patterns to get the counts of larger patterns. These two
algorithms are state-of-the-art for exact counting up to size $5$. Approximate
graphlet counting algorithms are developed based on techniques such as path
sampling (Jha et al., 2015), edge and local neighborhood sampling (Rossi et
al., 2018; Rahman et al., 2014), subgraph sampling (Wang et al., 2017).

Graphlet counting in dynamic network: Over the last few years, a specific type
of subgraph such as cycle (Chen et al., 2022; Shin et al., 2020; Stefani et
al., 2017), clique (Dhulipala et al., 2021; Das et al., 2020) has been studied
extensively. Hanauer, Henzinger, & Hua (Hanauer et al., 2022) first proposed
an algorithm for maintaining the number of all graphlets of size $4$ and
provided rigorous theoretical analysis on the bounds on update and query
complexity, update time, query time, and space complexity. Our work differs
from theirs in several aspects: (1) Our algorithm is practical - it is based
on one of the most efficient algorithms for exact counting $4$ node graphlets
in practice such as PGD. (2) Space complexity of our algorithm is typically
low as opposed to the work of Hanauer et al. where the space complexity is
$O(n^{2})$ as it requires maintaining several data structures along with the
dynamic graph in their work. In contrast, we only need to maintain the dynamic
graph. (3) Our algorithm works in a batch dynamic mode, meaning we process all
the edges in a batch at a time instead of updating one edge at a time.

##  3\. Problem Definition

We consider simple undirected graphs and model the dynamic graph as a sequence
of edge additions and deletions.

We consider monitoring a graph that changes over time. Assume that, for any
$t\geq 0$, $G^{t}=(V^{t},E^{t})$ be the graph up to time $t$ where $V^{t}$ is
the set of vertices, and $E^{t}$ is the set of edges. For any $t\geq 0$, at
time $t+1$, we add a tuple $<e,o>$ from the stream, where $o\in\\{+,-\\}$
represents an update process, add edge or delete edge, $e=(u,v)$ a pair of
vertices. A graph $G^{t+1}=(V^{t+1},E^{t+1})$ is generated by adding a new
edge or deleting an existing edge as follows:

$E^{t+1}=\begin{cases}E^{t}\cup\\{u\\}&\textit{if }o=+\\\
E^{t}\setminus\\{u\\}&\textit{if }o=-\end{cases}$

The addition and deletion of vertices can be handled similarly. Moreover, we
assume that by adding an edge $e=(u,v)$, the vertices $u$ and $v$ are added to
the graph if they are not present at time $t$. Furthermore, when a vertex is
deleted, we assume all the incident edges are deleted before deleting the
vertex. For simplicity, the rest of the paper deals with edge addition and
deletion - vertex operations can be handled easily by iteratively deleting all
the edges adjacent to a vertex.

The neighborhood of a vertex $u\in V^{t}$ at time $t$ is defined as
$N^{t}(u)=\\{v|(u,v)\in E^{t}\\}$. Similarly, the $k$-hop neighborhood of $u$
at time $t$ is denoted as $N^{t}_{k}(u)$. When the context is clear, we will
use $N(u)$ and $N_{k}(u)$ to mean neighborhood and $k$-hop neighborhood.

Given a simple undirected graph $G=(V,E)$, a graphlet of size $k$ nodes is
defined as a subgraph $g_{k}\subset G$ where $g_{k}$ has $k$ vertices. A
graphlet $g_{k}$ is an induced graphlet if $g_{k}$ is an induced subgraph of
$G$ and it is a connected graphlet if $g_{k}$ is a connected subgraph. Given a
subset $S\subset V$, $G[S]$ denotes the subgraph of $G$ induced by $S$.

Next, we define the problem of maintaining the frequency of connected induced
graphlets of size-$k$.

###### Problem 1.

Given an evolving graph $G^{t}=(V,E^{t})$, an integer $k$, and a batch of
edges $\mathcal{B}=\mathcal{A}\cup\mathcal{D}$, maintain the frequency of
connected induced size-$k$ graphlets in
$G^{t+1}=(V,E^{t}\cup\\{\mathcal{A}\\}\setminus\\{\mathcal{D}\\})$.

We focus on a fully dynamic setting when the addition and deletion of edges
are allowed in an arbitrary order in a batch. We aim to maintain the frequency
of all connected induced size-$4$ graphlets (shown in Figure 1) without
recomputing from scratch.

Figure 1. A list of different kinds of graphlets from the perspective of the
bottom edge. From top to bottom, left to right, the graphlets are as follows:
3-paths ($g_{1}^{4}$), a 3-star ($g_{2}^{4}$), a 4-cycle ($g_{3}^{4}$), tailed
triangles ($g_{4}^{4}$), diamonds ($g_{5}^{4}$), and a 4-clique ($g_{6}^{4}$).

In this work, we develop algorithms for the maintenance of size-$4$ graphlets.
In doing so, we identify the local subgraph and count graphlets in the local
subgraph using a state-of-the-art deterministic graphlet counting algorithm
such as PGD.

PGD Algorithm: The algorithm counts graphlets of size $3$ and $4$ in a static
simple undirected graph. The algorithm leverages combinatorial arguments for
different graphlets.

The algorithm iterates over the edges of the graph. Each iteration counts only
a few graphlets and uses combinatorial arguments to derive the exact count of
the rest of the graphlets in constant time. We leverage this algorithm to
maintain the graphlet count in a dynamic network by appropriately choosing the
subgraph and applying PGD on that subgraph with a guarantee that we can
maintain the global graphlet count by counting graphlets in that subgraph
only.

##  4\. Algorithms

In this section, we design efficient algorithms for maintaining size-$4$
graphlets in a dynamic graph. First, we will describe the algorithm for an
incremental batch where the goal is to update the graphlet count once a batch
of new edges is added to the graph. Next, we will describe the algorithm for a
fully dynamic batch consisting of edges for the addition and edges for the
deletion. The algorithm’s objective is to maintain the frequency of all
induced graphlets in the dynamic graph.

###  4.1. Baseline Algorithm

A straightforward algorithm for maintaining graphlet counts is to count all
the graphlets in the updated graph using a state-of-the-art graphlet counting
algorithm such as PGD (Ahmed et al., 2015b). We name this baseline as PGDN.
This counting-from-scratch strategy is inefficient when the batch of edges
touches a tiny portion of the graph. At the same time, the static algorithm
for graphlet counting has to work on the much larger region of the graph for
graphlet counting, which is required for maintaining the global count. An
efficient technique would be to update the graphlet based on the count from
the graph region where the change is located.

###  4.2. Incremental Stream

Now we will present an efficient algorithm for updating the graphlet count
once a batch of edges $\mathcal{B}$ is added to the original graph $G$ through
systematically exploring the subgraphs local to the changes in the graph. Once
the subgraphs are explored, we count the graphlets in the local subgraphs.
Finally, we aggregate these local counts to maintain the global graphlet
counts in the updated graph.

Suppose a batch of new edges $\mathcal{B}=\\{e_{1},e_{2},...,e_{l}\\}$ is
added to the graph $G=(V,E)$ and the updated graph is
$G^{\prime}=G+\mathcal{B}=(V,E\cup\mathcal{B})$. This algorithm maintains two
sub-graphs $\mathcal{G}_{a}$ and $\mathcal{G}_{b}$ of $\mathcal{G}$. It is
easy to see that the diameter of size-$4$ graphlet is $3$, which occurs in the
case of a $3$-path (Figure 1). Hence, the construction of $\mathcal{G}_{a}$ is
as follows: For each edge $e_{i}=(u_{i},v_{i})$, we construct a subgraph
$g_{i}$ induced by the vertex set $W_{i}=N_{3}(u_{i})\cup N_{3}(v_{i})$.
Finally, $\mathcal{G}_{a}=g_{1}\cup g_{2}\cup...\cup g_{l}$. We construct
$\mathcal{G}_{b}$ by removing all the edges in $\mathcal{B}$ from
$\mathcal{G}_{a}$, formally,
$\mathcal{G}_{b}=\mathcal{G}_{a}\setminus\mathcal{B}$.

The pseudocode is presented in Algorithm 1. The algorithm first updates the
original graph $\mathcal{G}$ by adding the batch $\mathcal{B}$ to it. Next, it
creates $\mathcal{G}_{a}$ and $\mathcal{G}_{b}$ and counts graphlets in
$\mathcal{G}_{a}$ and $\mathcal{G}_{b}$. The following lemma shows that we can
maintain the exact count of the graphlets using the graphlet counts in
$\mathcal{G}_{a}$ and $\mathcal{G}_{b}$.

1

2Initialization:

3 Global Graphlet Counts $\mathcal{F}_{4}$ = [0, 0, 0, 0, 0, 0]

4 $\mathcal{G}_{a}$ Graphlet Counts $\mathcal{C}_{a}$ = [0, 0, 0, 0, 0, 0]

5 $\mathcal{G}_{b}$ Graphlet Counts $\mathcal{C}_{b}$ = [0, 0, 0, 0, 0, 0]

6

7foreach  _Batch $\mathcal{B}\in Inputs$_ do

8

9 $\mathcal{G}\leftarrow\mathcal{G}+\mathcal{B}$

10 Generate graph $\mathcal{G}_{a}$ = ($\mathcal{V}_{a}$, $\mathcal{E}_{a}$)

11 Generate graph $\mathcal{G}_{b}$ = ($\mathcal{V}_{b}$, $\mathcal{E}_{b}$)

12

13 $\mathcal{C}_{a}\leftarrow PGD(\mathcal{G}_{a})$

14 $\mathcal{C}_{b}\leftarrow PGD(\mathcal{G}_{b})$

15 $\mathcal{C}^{\prime}\leftarrow\mathcal{C}_{a}-\mathcal{C}_{b}$

16 $\mathcal{F}_{4}\leftarrow\mathcal{F}_{4}+\mathcal{C}^{\prime}$

17

18 end foreach

19return $\mathcal{F}_{4}$

Algorithm 1 IGC: Incremental Graphlet Counting

###### Lemma 0.

Suppose a batch of new edges $\mathcal{B}$ is added to the graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ and the graph is updated to
$\mathcal{G}^{\prime}$. Then,
$\mathcal{F}^{\mathcal{G}^{\prime}}_{4}=\mathcal{F}^{\mathcal{G}}_{4}+\mathcal{C}_{a}-\mathcal{C}_{b}$.

###### Proof.

There are three cases to consider:

case-1: all graphlets containing at least one new edge are added to the count.
Observe that $\mathcal{C}_{a}$ captures all the graphlets that contain at
least one new edge. Further, $\mathcal{C}_{b}$ will never capture any of these
graphlets since they are missing some edges in $\mathcal{G}_{b}$. So, the
count of every graphlet containing at least a new edge is added exactly once
to the global count to update the answer.

case-2: all graphlets extended to a new graphlet by addition of edges are
removed from the count. Observe that $\mathcal{C}_{b}$ captures all the
graphlets as we remove all the new edges from $\mathcal{G}_{a}$ to get
$\mathcal{G}_{b}$ and $\mathcal{C}_{a}$ will never capture the count of any
such graphlets as they become non-induced once a batch is added. So, the count
of these graphlets is subtracted exactly once.

case-3: All graphlets that do not change after the addition of new edges do
not contribute to the changes in the global graphlet counts. Since the
graphlets were present before and after the addition of edges, they must be
captured by both $\mathcal{C}_{a}$ and $\mathcal{C}_{b}$. So the count of
every such graphlet is added exactly once and subtracted exactly once. So, the
global count does not change. ∎

###  4.3. Fully Dynamic Stream

In this section, we describe our algorithm for maintaining global graphlet
count in a fully dynamic setting where both insertion and deletion of edges
are possible.

Suppose $\mathcal{B}=\mathcal{A}\cup\mathcal{D}$ is a batch of edges where
$\mathcal{A}$ is a batch of edges to add and $\mathcal{D}$ is a batch of edges
to delete. and Let,
$H=\cup_{e=(u,v)\in\mathcal{A}\cup\mathcal{D}}G[\mathcal{N}_{3}(u)\cup\mathcal{N}_{3}(v)]$.
Next we define $\mathcal{G}_{a}=(H\setminus\mathcal{D})\cup\mathcal{A}$ and
$\mathcal{G}_{b}=(H\setminus\mathcal{A})\cup\mathcal{D}$ and then count the
number of graphlets in $\mathcal{G}_{a}$ and $\mathcal{G}_{b}$. We use these
counts in a systematic way to update the global graphlet count after
processing the batch $\mathcal{B}$. The pseudocode is presented in Algorithm 2

1

2Initialization:

3 Global Graphlet Counts $\mathcal{F}_{4}$ = [0, 0, 0, 0, 0, 0]

4 $\mathcal{G}_{a}$ Graphlet Counts $\mathcal{C}_{a}$ = [0, 0, 0, 0, 0, 0]

5 $\mathcal{G}_{b}$ Graphlet Counts $\mathcal{C}_{b}$ = [0, 0, 0, 0, 0, 0]

6

7foreach  _Batch $\mathcal{B}=\mathcal{A}\cup\mathcal{D}\in Inputs$_ do

8

9 $H=\cup_{e\in\mathcal{A}\cup\mathcal{D}}G[\mathcal{N}_{3}(e)]$

10 $\mathcal{G}_{a}=(H\setminus\mathcal{D})\cup\mathcal{A}$

11 $\mathcal{G}_{b}=(H\setminus\mathcal{A})\cup\mathcal{D}$

12

13 $\mathcal{C}_{a}\leftarrow PGD(\mathcal{G}_{a})$

14 $\mathcal{C}_{b}\leftarrow PGD(\mathcal{G}_{b})$

15 $\mathcal{C}^{\prime}\leftarrow\mathcal{C}_{a}-\mathcal{C}_{b}$

16 $\mathcal{F}_{4}\leftarrow\mathcal{F}_{4}+\mathcal{C}^{\prime}$

17 $\mathcal{G}\leftarrow\mathcal{G}\setminus\mathcal{D}\cup\\{\mathcal{A}\\}$

18

19 end foreach

20return $\mathcal{F}_{4}$

Algorithm 2 FDGC: Fully Dynamic Graphlet Counting

###### Lemma 0.

Algorithm 2 correctly updates the global graphlet frequency.

###### Proof.

There are four types of graphlets. (1) type-1 graphlets containing edges in
$\mathcal{D}$ but no edges in $\mathcal{A}$; (2) type-2 graphlets containing
edges in both $\mathcal{A}$ and $\mathcal{D}$ (3) type-3 graphlets containing
edges in $\mathcal{A}$ but no edge in $\mathcal{D}$ (4) type-4 graphlets
containing edges neither from $\mathcal{A}$ nor from $\mathcal{D}$. All type-4
graphlets are captured by both $\mathcal{C}_{a}$ and $\mathcal{C}_{b}$, So
$\mathcal{C}_{a}-\mathcal{C}_{b}=0$. Counts due to type-3 graphlets are
captured in $\mathcal{C}_{a}$, added to the global count, but they are not
captured by $\mathcal{C}_{b}$. type-2 graphlets are captured once in
$\mathcal{C}_{a}$ and once in $\mathcal{C}_{b}$. So
$\mathcal{C}_{a}-\mathcal{C}_{b}=0$. type-1 graphlets are removed, and
therefore, corresponding counts are deleted, captured in $\mathcal{C}_{b}$,
but not captured in $\mathcal{C}_{a}$. Transitioning of the graphlets as in
case-2 of Lemma 4.1 are handled similarly as in Lemma 4.1. This completes the
proof.

∎

##  5\. Experiments

We evaluated our algorithm on large datasets to show the efficiency over the
baseline PGDN. Every time we add a batch of edges and perform recomputation
from scratch using the PGD. We show the description of the dataset in Table 1
taken from the Stanford SNAP (Leskovec and Krevl, 2014). We implemented all
the algorithms in C++ and evaluated them in a system equipped with Intel(R)
Xeon(R) W-2223 CPU $@$ $3.60$GHz with $4$ physical cores and $32$G RAM. We
execute PGD with OpenMP to run in parallel.

Table 1.  Characteristics of the Data sets

Dataset | $|V|$ | $|E|$
---|---|---
WikiTalk (WT) | $2394385$ | $4659565$
WikiVote (WV) | $7115$ | $100762$
Soc-Pokec (SP) | $1632803$ | $22301964$
soc-LiveJournal (LJ) | $4033137$ | $27933062$

Batch creation: For the incremental batch computation, we group edges in a
batch and execute the algorithms by inserting one batch at a time. For the
fully dynamic batch, we create the batches as follows: First, we give $+$ive
and $-$ive labels to the edges by sampling each edge with probability $p$ for
selecting it to be a positive edge and putting the sampled edges in a queue.
Later, we extract edges from the queue and sample them with probability $1-p$
to assign a $-$ive sign. This way, we create both positive and negative edges
where the positive edge refers to the addition and the negative refers to the
deletion. Next, we group them as before for creating batches. For our
experiment, we choose $p=0.7$.

Computation on incremental batch: We empirically evaluated our algorithm IGC
to show that IGC is much faster than PGDN on large graphs such as Soc-Pokec
and soc-LiveJournal. This is because IGC counts graphlets on a relatively
small size subgraph, which is local to the changes in the graph due to the
batch update, while PGDN access the entire network built so far to count the
graphlets. Also, the performance of IGC and PGDN is almost similar on WikiVote
(Figure 2(a)) and WikiTalk (Figure 2(b)) because these networks are relatively
small and sparse. Also, IGC needs to compute $3$-hop neighborhood for
constructing the local subgraph, which is not substantially smaller than the
original network where PGD work in the baseline. The results on batches of
size $10$ and $100$ are shown in Figure 2 and Figure 3, respectively. In both
cases, we observe that as the size of the graph increases with the increase in
the edge counts, IGC is substantially faster than the baseline PGDN in two
large networks Soc-Pokec (Figure 2(c) and Figure 3(b)) and soc-LiveJournal
(Figure 2(d) and Figure 3(a)). This is as expected because PGDN accesses a
larger graph each time compared to the size of the graph accessed by IGC.

Computation on fully dynamic batch: In a fully dynamic batch, we first removed
the edges, which are both added and deleted, then executed Algorithm FDGC on
the rest of the edges. Similar to the incremental stream, we observed similar
behavior on two large networks Soc-Pokec and soc-LiveJournal as shown in
Figure 4. Similar to the earlier observation, as the network size increases,
the runtime of FDGC is significantly faster compared to PGDN than the initial
state of the graph.

(a) WikiVote

(b) WikiTalk

(c) Soc-Pokec

(d) soc-LiveJournal

Figure 2.  Cumulative time to add insert only batches (of size $10$). Time at
$y$-axis shows the overall time taken to add all the batches so far
(represented by the batch number on the $x$-axis)

(a) soc-LiveJournal

(b) Soc-Pokec

Figure 3. Cumulative time to add insert only batches (of size $100$). Time at
$y$-axis shows the overall time taken to add all the batches so far
(represented by the batch number on the $x$-axis)

(a) soc-LiveJournal

(b) Soc-Pokec

Figure 4. Cumulative time to add fully-dynamic batches (of size $100$). Time
at $y$-axis shows the overall time taken to add all the batches so far
(represented by the batch number on the $x$-axis)

##  6\. Conclusion

In this work, we have developed an efficient algorithm for updating the counts
of size-$4$ induced graphlets. Our algorithm is based on the exploration of
subgraphs local to change in the graph by the fully dynamic stream of edges.
We experimentally evaluated our algorithm to show the efficiency over the
baseline.

## References

  * (1)
  * Ahmed et al. (2015a) Nesreen K Ahmed, Jennifer Neville, Ryan A Rossi, and Nick Duffield. 2015a.  Efficient graphlet counting for large networks. In _2015 IEEE international conference on data mining_. IEEE, 1–10.
  * Ahmed et al. (2015b) Nesreen K Ahmed, Jennifer Neville, Ryan A Rossi, and Nick Duffield. 2015b.  Efficient graphlet counting for large networks. In _2015 IEEE International Conference on Data Mining_. IEEE, 1–10.
  * Bera and Seshadhri (2020) Suman K Bera and C Seshadhri. 2020.  How the degeneracy helps for triangle counting in graph streams. In _Proceedings of the 39th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems_. 457–467.
  * Chen et al. (2022) Justin Y Chen, Talya Eden, Piotr Indyk, Honghao Lin, Shyam Narayanan, Ronitt Rubinfeld, Sandeep Silwal, Tal Wagner, David P Woodruff, and Michael Zhang. 2022\.  Triangle and Four Cycle Counting with Predictions in Graph Streams. In _International Conference on Learning Representations_.
  * Danisch et al. (2018) Maximilien Danisch, Oana Balalau, and Mauro Sozio. 2018\.  Listing k-cliques in sparse real-world graphs. In _Proceedings of the 2018 World Wide Web Conference_. 589–598.
  * Das et al. (2020) Apurba Das, Seyed-Vahid Sanei-Mehri, and Srikanta Tirthapura. 2020\.  Shared-memory parallel maximal clique enumeration from static and dynamic graphs.  _ACM Transactions on Parallel Computing (TOPC)_ 7, 1 (2020), 1–28.
  * Dhulipala et al. (2021) Laxman Dhulipala, Quanquan C Liu, Julian Shun, and Shangdi Yu. 2021.  Parallel batch-dynamic k-clique counting. In _Symposium on Algorithmic Principles of Computer Systems (APOCS)_. SIAM, 129–143.
  * Hanauer et al. (2022) Kathrin Hanauer, Monika Henzinger, and Qi Cheng Hua. 2022\.  Fully Dynamic Four-Vertex Subgraph Counting. In _1st Symposium on Algorithmic Foundations of Dynamic Networks_.
  * Harchaoui and Bach (2007) Zaid Harchaoui and Francis Bach. 2007.  Image classification with segmentation graph kernels. In _2007 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE, 1–8.
  * Janssen et al. (2012) Jeannette Janssen, Matt Hurshman, and Nauzer Kalyaniwalla. 2012\.  Model selection for social networks using graphlets.  _Internet Mathematics_ 8, 4 (2012), 338–363.
  * Jha et al. (2015) Madhav Jha, C Seshadhri, and Ali Pinar. 2015.  Path sampling: A fast and provable method for estimating 4-vertex subgraph counts. In _Proceedings of the 24th international conference on world wide web_. 495–505.
  * Leskovec and Krevl (2014) Jure Leskovec and Andrej Krevl. 2014.  SNAP Datasets: Stanford Large Network Dataset Collection.  http://snap.stanford.edu/data.
  * McGregor and Vorotnikova (2020) Andrew McGregor and Sofya Vorotnikova. 2020.  Triangle and four cycle counting in the data stream model. In _Proceedings of the 39th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems_. 445–456.
  * Milenković et al. (2010) Tijana Milenković, Vesna Memišević, Anand K Ganesan, and Nataša Pržulj. 2010\.  Systems-level cancer gene identification from protein interaction network topology applied to melanogenesis-related functional genomics data.  _Journal of the Royal Society Interface_ 7, 44 (2010), 423–437.
  * Milenković and Pržulj (2008) Tijana Milenković and Nataša Pržulj. 2008\.  Uncovering biological network function via graphlet degree signatures.  _Cancer informatics_ 6 (2008), CIN–S680.
  * Mukherjee et al. (2018) Kingshuk Mukherjee, Md Mahmudul Hasan, Christina Boucher, and Tamer Kahveci. 2018.  Counting motifs in dynamic networks.  _BMC systems biology_ 12, 1 (2018), 1–12.
  * Pinar et al. (2017) Ali Pinar, C Seshadhri, and Vaidyanathan Vishal. 2017\.  Escape: Efficiently counting all 5-vertex subgraphs. In _Proceedings of the 26th international conference on world wide web_. 1431–1440.
  * Pržulj et al. (2004) Natasa Pržulj, Derek G Corneil, and Igor Jurisica. 2004\.  Modeling interactome: scale-free or geometric?  _Bioinformatics_ 20, 18 (2004), 3508–3515.
  * Rahman et al. (2014) Mahmudur Rahman, Mansurul Alam Bhuiyan, and Mohammad Al Hasan. 2014.  Graft: An efficient graphlet counting method for large graph analysis.  _IEEE Transactions on Knowledge and Data Engineering_ 26, 10 (2014), 2466–2478.
  * Rossi et al. (2018) Ryan A Rossi, Rong Zhou, and Nesreen K Ahmed. 2018.  Estimation of graphlet counts in massive networks.  _IEEE transactions on neural networks and learning systems_ 30, 1 (2018), 44–57.
  * Saha and Hasan (2015) Tanay Kumar Saha and Mohammad Al Hasan. 2015.  Finding network motifs using MCMC sampling.  (2015), 13–24.
  * Seshadhri et al. (2013) Comandur Seshadhri, Ali Pinar, and Tamara G Kolda. 2013\.  Triadic measures on graphs: The power of wedge sampling. In _Proceedings of the 2013 SIAM international conference on data mining_. SIAM, 10–18.
  * Shervashidze et al. (2009) Nino Shervashidze, SVN Vishwanathan, Tobias Petri, Kurt Mehlhorn, and Karsten Borgwardt. 2009\.  Efficient graphlet kernels for large graph comparison. In _Artificial intelligence and statistics_. PMLR, 488–495.
  * Shin et al. (2020) Kijung Shin, Sejoon Oh, Jisu Kim, Bryan Hooi, and Christos Faloutsos. 2020.  Fast, accurate and provable triangle counting in fully dynamic graph streams.  _ACM Transactions on Knowledge Discovery from Data (TKDD)_ 14, 2 (2020), 1–39.
  * Stefani et al. (2017) Lorenzo De Stefani, Alessandro Epasto, Matteo Riondato, and Eli Upfal. 2017.  Triest: Counting local and global triangles in fully dynamic streams with fixed memory size.  _ACM Transactions on Knowledge Discovery from Data (TKDD)_ 11, 4 (2017), 1–50.
  * Vacic et al. (2010) Vladimir Vacic, Lilia M Iakoucheva, Stefano Lonardi, and Predrag Radivojac. 2010.  Graphlet kernels for prediction of functional residues in protein structures.  _Journal of Computational Biology_ 17, 1 (2010), 55–72.
  * Wang et al. (2012) Jianxin Wang, Yuannan Huang, Fang-Xiang Wu, and Yi Pan. 2012\.  Symmetry compression method for discovering network motifs.  _IEEE/ACM Transactions on Computational Biology and Bioinformatics_ 9, 6 (2012), 1776–1789.
  * Wang et al. (2014) Pinghui Wang, John CS Lui, Bruno Ribeiro, Don Towsley, Junzhou Zhao, and Xiaohong Guan. 2014\.  Efficiently estimating motif statistics of large networks.  _ACM Transactions on Knowledge Discovery from Data (TKDD)_ 9, 2 (2014), 1–27.
  * Wang et al. (2017) Pinghui Wang, Junzhou Zhao, Xiangliang Zhang, Zhenguo Li, Jiefeng Cheng, John CS Lui, Don Towsley, Jing Tao, and Xiaohong Guan. 2017.  MOSS-5: A fast method of approximating counts of 5-node graphlets in large graphs.  _IEEE Transactions on Knowledge and Data Engineering_ 30, 1 (2017), 73–86.
  * Zhang et al. (2013) Luming Zhang, Mingli Song, Zicheng Liu, Xiao Liu, Jiajun Bu, and Chun Chen. 2013\.  Probabilistic graphlet cut: Exploiting spatial structure cue for weakly supervised image segmentation. In _Proceedings of the IEEE conference on computer vision and pattern recognition_. 1908–1915.
  * Zhang et al. (2012) Luming Zhang, Mingli Song, Qi Zhao, Xiao Liu, Jiajun Bu, and Chun Chen. 2012\.  Probabilistic graphlet transfer for photo cropping.  _IEEE Transactions on Image Processing_ 22, 2 (2012), 802–815.


# Searching for Afterglow: Light Dark Matter Boosted by Supernova Neutrinos

Yen-Hsun Lin<EMAIL_ADDRESS>Physics Division, National Center for
Theoretical Sciences, Taipei 106, Taiwan  Institute of Physics, Academia
Sinica, Taipei 115, Taiwan  Wen-Hua Wu  Department of Physics, National Taiwan
University, Taipei 106, Taiwan  Institute of Physics, Academia Sinica, Taipei
115, Taiwan  Meng-Ru Wu<EMAIL_ADDRESS>Institute of Physics, Academia
Sinica, Taipei 115, Taiwan  Institute of Astronomy and Astrophysics, Academia
Sinica, Taipei 106, Taiwan  Physics Division, National Center for Theoretical
Sciences, Taipei 106, Taiwan  Henry Tsz-King Wong<EMAIL_ADDRESS>Institute of Physics, Academia Sinica, Taipei 115, Taiwan

###### Abstract

A novel analysis is performed, incorporating _time-of-flight_ (TOF)
information to study the interactions of dark matter (DM) with standard model
particles. After supernova (SN) explosions, DM with mass
$m_{\chi}\lesssim\mathcal{O}({\rm MeV})$ in the halo can be boosted by SN
neutrinos (SN$\nu$) to relativistic speed. The SN$\nu$ boosted DM (BDM)
arrives on Earth with TOF which depends only on $m_{\chi}$ and is independent
of the cross section. These BDMs can interact with detector targets in low-
background experiments and manifest as _afterglow_ events after the arrival of
SN$\nu$. The characteristic TOF spectra of the BDM events can lead to large
background suppression and unique determination of $m_{\chi}$. New cross
section constraints on $\sqrt{\sigma_{\chi e}\sigma_{\chi\nu}}$ are derived
from SN1987a in the Large Magellanic Cloud with data from the Kamiokande and
Super-Kamiokande experiments. Potential sensitivities for the next galactic SN
with Hyper-Kamiokande are projected. This analysis extends the existing bounds
on $\sqrt{\sigma_{\chi e}\sigma_{\chi\nu}}$ over a broad range of
$r_{\chi}=\sigma_{\chi\nu}/\sigma_{\chi e}$. In particular, the improvement is
by 1–3 orders of magnitude for $m_{\chi}<\mathcal{O}(100\,{\rm keV})$ for
$\sigma_{\chi e}\sim\sigma_{\chi\nu}$. Prospects of exploiting TOF information
in other astrophysical systems to probe exotic physics with other DM
candidates are discussed.

Introduction.—Although there is compelling evidence on the existence of dark
matter (DM) as an additional gravity source, its properties and interactions
remain unknown Battaglieri  _et al._ (2017); Workman  _et al._ (2022).
Experimental searches of DM are intensely pursued worldwide Aguilar  _et al._
(2015); Akerib  _et al._ (2017a); Ackermann  _et al._ (2017); Akerib  _et al._
(2017b); Aprile  _et al._ (2018); Ambrosi  _et al._ (2017); Aprile  _et al._
(2019a, b); Abramoff  _et al._ (2019); Agnese  _et al._ (2018). Direct
detection (DD) experiments focus on the weakly interacting massive particle
(WIMP) scenario of DM mass $m_{\chi}\gtrsim\mathcal{O}({\rm GeV})$, with
sensitivities approaching the neutrino floor Workman  _et al._ (2022). The
search for lighter WIMPs is an active area of research. One scenario with
rapidly expanding interest is where light DM is upscattered or boosted by
known cosmic particles including baryons, electrons, and neutrinos Bringmann
and Pospelov (2019); Ema  _et al._ (2019); Cappiello and Beacom (2019); Dent
_et al._ (2020); Wang  _et al._ (2020); Zhang (2022); Guo  _et al._ (2020a);
Ge  _et al._ (2021); Cao  _et al._ (2021); Jho  _et al._ (2020); Cho  _et al._
(2020); Lei  _et al._ (2020); Guo  _et al._ (2020b); Xia  _et al._ (2021);
Dent  _et al._ (2021); Ema  _et al._ (2021); Flambaum  _et al._ (2020); Jho
_et al._ (2021); Das and Sen (2021); Bell  _et al._ (2021); Chao  _et al._
(2021); Ghosh  _et al._ (2022); Feng  _et al._ (2022); Wang  _et al._ (2021);
Xia  _et al._ (2022); Wang  _et al._ (2022); Cui  _et al._ (2022); Xu  _et
al._ (2022); Granelli  _et al._ (2022). The boosted DM (BDM) then carries
kinetic energy $T_{\chi}$ much larger than when it is nonrelativistic (with
velocity $v_{\chi}\sim 10^{-3}$) according to the Halo model. Nuclear and
electron recoil events from BDM interaction with the detector targets will
therefore have increased energy deposition, making DM with
$m_{\chi}\lesssim\mathcal{O}({\rm GeV})$ experimentally accessible.

_Time-of-flight_ (TOF) techniques are matured laboratory tools for
differentiation or measurement of particle masses. This technique, however,
has not been well exploited to probe exotic physics in astrophysical systems.
One notable exception is the neutrino mass constraints derived from the timing
distributions of supernova neutrinos (SN$\nu$s) from SN1987a Arnett and Rosner
(1987); Kolb  _et al._ (1987). We explore in this Letter a novel scenario of
BDM with kinetic energy injected by SN$\nu$ interactions, and in particular
where the prompt SN$\nu$ burst is also detected, providing a time-zero
definition in terrestrial experiments. The prompt SN$\nu$ events will be
followed by time-evolving BDM afterglow events where energy and time can be
measured. The delay time between BDM and SN$\nu$ is a distinctive “smoking-
gun” signature and provides unique information to infer $m_{\chi}$,
independent of the interaction cross section. Specifically, a delay time of
$\Delta t\simeq 10$ days$\times[R/(8~{}{\rm kpc})][m_{\chi}/(10~{}{\rm
keV})]^{2}[T_{\chi}/(10~{}{\rm MeV})]^{-2}$ for SN$\nu$BDM traveling an
astronomical distance $R$ before reaching the Earth highlights that although
BDM has $v_{\chi}\sim c$, the delay can be substantial but measurable in a
duration post the arrival of SN$\nu$. In contrast, most proposed BDM scenarios
rely on steady sources, e.g., cosmic rays Bringmann and Pospelov (2019); Ema
_et al._ (2019); Cappiello and Beacom (2019); Dent  _et al._ (2020); Wang  _et
al._ (2020); Guo  _et al._ (2020a); Ge  _et al._ (2021); Cao  _et al._ (2021);
Jho  _et al._ (2020); Cho  _et al._ (2020); Lei  _et al._ (2020); Guo  _et
al._ (2020b); Xia  _et al._ (2021); Dent  _et al._ (2021); Ema  _et al._
(2021); Flambaum  _et al._ (2020); Bell  _et al._ (2021); Feng  _et al._
(2022); Wang  _et al._ (2021); Xia  _et al._ (2022), stellar $\nu$ Zhang
(2022); Jho  _et al._ (2021), diffuse SN$\nu$ Das and Sen (2021); Ghosh  _et
al._ (2022), etc, for which the BDM flux is constant with time and lacks any
time-dependent feature.

We explore the signatures of SN$\nu$BDM with SN1987a in the Large Magellanic
Cloud (LMC) and a future supernova (SN) in the Galactic Center (GC) to derive
the fluxes and the associated electron-recoil event rates via $\sigma_{\chi
e}$ in multikiloton water Cherenkov detectors, including Kamiokande, Super-
Kamiokande (Super-K), and Hyper-Kamiokande. The scenario of SN$\nu$BDM depends
on finite DM cross section with $\nu$ ($\sigma_{\chi\nu}$), which may
originate from an effective Lagrangian
$\bar{\chi}\Gamma{\chi}\bar{\ell_{i}}\Gamma\ell_{i}/\Lambda$ where $\chi$ and
$\ell_{i}=(\nu_{i},i)$ are the DM and SM fields with $i=e,\mu,\tau$. The
vertex $\Gamma$ denotes the interaction type and $\Lambda$ indicates certain
cutoff scale. Possible interactions between $\chi$ with $\nu$ is a subject of
intense recent interest Holdom (1986); Argüelles  _et al._ (2017); Chang  _et
al._ (2018); Murase and Shoemaker (2019); Lin and Lin (2021); Croon  _et al._
(2021); Escudero  _et al._ (2019); Foldenauer (2019). They can naturally arise
in many particle physics models such as the extensively studied $B-L$ and
$L_{\mu}-L_{\tau}$, where the new gauge bosons can kinematically mix with the
standard model photon. Further constraints will be provided by this work.

DM boosted by SN$\nu$.—Assuming a SN explodes near the center of a galaxy
(location $O$ in Fig. 1), it emits a large amount of $\mathcal{O}(10)$ MeV
neutrinos within $\tau\approx 10$ s carrying total luminosity $L_{\nu,{\rm
tot}}\approx 3\times 10^{52}$ erg s-1. We approximate these SN$\nu$ by an
expanding thin spherical shell with a radius $r$ away from $O$ and a thickness
$d\approx c\tau$ (see Fig. 1). The radially propagating SN$\nu$ within the
shell has a number density of

$\frac{dn_{\nu}}{dE_{\nu}}=\sum_{i}\frac{L_{\nu_{i}}}{4\pi r^{2}\langle
E_{\nu_{i}}\rangle}E_{\nu}^{2}f_{\nu_{i}}(E_{\nu}),$ (1)

where $L_{\nu_{i}}=L_{\nu,{\rm tot}}/6$ is the luminosity of each flavor
($\nu_{e}$, $\nu_{\mu}$, $\nu_{\tau}$ and their antineutrinos). We take the
average energy $\langle E_{\nu_{e}}\rangle$, $\langle
E_{\bar{\nu}_{e}}\rangle$, and $\langle E_{\nu_{x}}\rangle$
($\nu_{x}\in\\{\nu_{\mu},\nu_{\tau},\bar{\nu}_{\mu},\bar{\nu}_{\tau}\\}$) to
be 11, 16, 25 MeV, respectively Duan  _et al._ (2006). The energy distribution
follows a Fermi-Dirac distribution $f_{\nu_{i}}$ with a pinch parameter
$\eta_{\nu_{i}}\equiv\mu_{\nu_{i}}/T_{\nu_{i}}=3$, such that
$T_{\nu_{i}}\approx\langle E_{\nu_{i}}\rangle/3.99$.

Figure 1:  Schematic plot of DM boosted by SN$\nu$ within an expanding
spherical shell with width $d$ at radius $r$. The SN occurs at $O$. BDM from
$A$ arrives $B$ with an scattering angle $\alpha$.

With a nonvanishing DM-$\nu$ interaction, these neutrinos can upscatter DM in
the halo [with number density $n_{\chi}(r)$] when they propagate outward. The
BDM from location $A$ can reach the Earth at $B$ (with a distance $R$ away
from the center) with a scattering angle $\alpha$ after traveling a length
$\ell$. At neutrino energy $E_{\nu}$ much larger than the typical DM kinetic
energy in the halo, DM can be approximated as at rest, and the BDM kinetic
energy is given by

$T_{\chi}=\frac{E_{\nu}^{2}}{E_{\nu}+m_{\chi}/2}\left(\frac{1+\cos\theta_{c}}{2}\right),$
(2)

where $\theta_{c}\in[0,\pi]$ is the scattering angle in the center-of-mass
(c.m.) frame. One can relate $\theta_{c}$ to the lab frame scattering angle
$\alpha\in[0,\pi/2]$ by $\theta_{c}=2\tan^{-1}(\gamma\tan\alpha)$ and
$\gamma=(E_{\nu}+m_{\chi})/\sqrt{m_{\chi}(2E_{\nu}+m_{\chi})}$. Assuming
$\sigma_{\chi\nu}$ is independent of $\theta_{c}$ in the c.m. frame, the
normalized BDM angular distribution in the lab frame is given by

$f_{\chi}(\alpha,E_{\nu})=\frac{\gamma^{2}\sec^{3}\alpha}{\pi(1+\gamma^{2}\tan^{2}\alpha)^{2}},$
(3)

such that $\int d\Omega_{\alpha}f_{\chi}(\alpha,E_{\nu})=1$ for any given
$E_{\nu}$, where $d\Omega_{\alpha}=2\pi\sin\alpha d\alpha$. In Fig. 2, we plot
$2\pi\sin\alpha f_{\chi}(\alpha)$ for a fixed $T_{\chi}=10\,{\rm MeV}$
(corresponding to different $E_{\nu}$) with different $m_{\chi}$. It shows
that for BDM with $m_{\chi}/T_{\chi}\ll 1$, they are confined within a small
scattering angle relative to the direction of SN$\nu$.

Figure 2:  The BDM angular distribution $f_{\chi}(\alpha)$ times
$2\pi\sin\alpha$ for fixed $T_{\chi}=10$ MeV and different $m_{\chi}$.

The BDM emissivity $j_{\chi}$ at location $A$ can be written as

$j_{\chi}(r,T_{\chi},\alpha)=c\sigma_{\chi\nu}n_{\chi}\left(\frac{dn_{\nu}}{dE_{\nu}}\right)\left(\frac{dE_{\nu}}{dT_{\chi}}\frac{v_{\chi}}{c}\right)f_{\chi},$
(4)

where the BDM velocity
$v_{\chi}/c=\sqrt{T_{\chi}(2m_{\chi}+T_{\chi})}/(m_{\chi}+T_{\chi})$, and can
be evaluated using Eqs. (1) to (3).

Time-dependent BDM flux at Earth.—To obtain the BDM flux (number of BDM per
unit time per unit energy per solid angle) at Earth
$d\Phi_{\chi}/(dT_{\chi}d\Omega)$ (location $B$ in Fig. 1), we shall integrate
all $j_{\chi}$ along the line of sight $\ell$,

$\displaystyle\frac{d\Phi_{\chi}(T_{\chi},\theta,t^{\prime})}{dT_{\chi}d\Omega}$
$\displaystyle=\int d\ell
j_{\chi}(r,T_{\chi},\alpha)H\left(t^{\prime}-\frac{r}{c}-\frac{\ell}{v_{\chi}}\right)$
$\displaystyle\quad\times
H\left(\frac{r}{c}+\frac{\ell}{v_{\chi}}+\tau-t^{\prime}\right),$ (5)

where $d\Omega=2\pi\sin\theta d\theta$ is viewed from $B$. The Heaviside
functions limit $j_{\chi}$ to being nonzero only within the spherical shell of
width $d$ where SN$\nu$ are present. The arrival time of BDM, $t^{\prime}$,
relative to the time of SN explosion, includes the propagation time of SN$\nu$
from $O$ to $A$ ($r/c$) and the traveling time of BDM from $A$ to $B$
($\ell/v_{\chi}$).

Integrating Eq. (5) over $d\Omega$ and approximating
$H(x-x_{0})H(x_{0}+\epsilon-x)\sim\epsilon\delta(x_{0})$ for $\epsilon\ll
x_{0}$, we obtain

$\displaystyle\frac{d\Phi_{\chi}(T_{\chi},t^{\prime})}{dT_{\chi}}$
$\displaystyle=2\pi\tau\int d\cos\theta d\ell
j_{\chi}(r,T_{\chi},\alpha)\delta\left(t^{\prime}-\frac{r}{c}-\frac{\ell}{v_{\chi}}\right)$
$\displaystyle=\left.2\pi\tau\int_{0}^{1}d\cos\theta\mathcal{J}j_{\chi}(r,T_{\chi},\alpha)\right|_{t^{\prime}=\frac{r}{c}+\frac{\ell}{v_{\chi}}},$
(6)

where

$\mathcal{J}=\left(\frac{\ell-R\cos\theta}{rc}+\frac{1}{v_{\chi}}\right)^{-1}$
(7)

appears due to the change of variable $d\ell=\mathcal{J}dt^{\prime}$. Note
that for a given $(t^{\prime},\theta)$, one can find a unique solution of
$(r,\ell,\alpha)$ and compute the integration.

Figure 3: The BDM flux at Earth vs. $t$ with different $m_{\chi}$ for
$T_{\chi}=10\,{\rm MeV}$ and $\sigma_{\chi\nu}=10^{-35}\,{\rm cm^{2}}$. Fluxes
resulting from a SN in the GC and from SN1987a in LMC are shown with solid and
dash-dotted lines. The black dashed line indicates the maximum exposure time
$t_{\rm exp}=35$ years (see text for details).

BDM flux from SN in the GC and LMC.—We now compute the BDM fluxes at the Earth
from SN1987a in LMC and from a SN in the GC. We characterize $n_{\chi}$ in the
Milky Way (MW) and LMC by Navarro-Frenk-White (NFW) and Hernquist profiles
respectively. Both share the same expression

$n_{\chi}(r)=\frac{\rho_{s}}{m_{\chi}}\frac{1}{\frac{r}{r_{s}}(1+\frac{r}{r_{s}})^{n}},$
(8)

with $(n,\rho_{s},r_{s})=(2,184\,{\rm MeV\,cm^{-3}},24.4\,{\rm kpc})$ for MW
Bertone  _et al._ (2005) and $(n,\rho_{s},r_{s})=(3,68\,{\rm
MeV\,cm^{-3}},31.9\,{\rm kpc})$ for LMC Erkal  _et al._ (2019). The distances
$R$ for the two are $(R_{\rm GC},R_{\rm LMC})=(8.5,50)\,{\rm kpc}$. We neglect
the contribution from $r<10^{-5}\,{\rm kpc}$ since the profile in the inner
region is highly uncertain and the adopted profile diverges when $r\rightarrow
0$.

Fig. 3 shows $d\Phi_{\chi}/dT_{\chi}$ versus $t>\tau$ for $T_{\chi}=10$ MeV
with different $m_{\chi}$ for SN in the GC (solid) and in LMC (dot-dashed),
assuming $\sigma_{\chi\nu}=10^{-35}$ cm2. Note that we define a shifted time
coordinate $t=t^{\prime}-R/c$ as the delayed arrival time for BDM relative to
SN$\nu$. For $m_{\chi}=1$ keV and $1$ MeV, the most prominent feature is that
the BDM fluxes contain a rising part and peak at $t_{p}\approx
R(1/v_{\chi}-1/c)$. This is mainly due to the increase of $n_{\chi}\propto
r^{-1}$ toward the halo center. The postpeak tails are due to BDM
contributions with larger scattering angles. For $m_{\chi}=1\,{\rm eV}$,
$t_{p}\approx 0.004\,{\rm s}$ is too short and overlaps with the 10 s duration
of SN$\nu$ to be shown in Fig. 3. Comparing BDM fluxes coming from the GC to
LMC, the LMC cases have smaller fluxes and larger $t_{p}$ due to larger $R$
and smaller halo density.

Fig. 3 also shows another important feature—the BDM flux for a given
$T_{\chi}$ and $m_{\chi}$ vanishes after some time post $t_{p}$, which is
related to the sharp cutoff of $f_{\chi}$ shown in Fig. 2. This allows us to
consider a reduced duration for BDM searches after the arrival of SN$\nu$.
Practically, a detector that can probe BDM has a threshold energy $T_{\rm
th}$, below which the detector is insensitive to BDM. Thus, for a given
$m_{\chi}$, one can define the latest possible arrival time of BDM with
$T_{\chi}=T_{\rm th}$ as the vanishing time $t_{\rm van}$ to analyze the data.
We stress that all these time-dependent features only depend on $m_{\chi}$ but
not $\sigma_{\chi\nu}$. Consequently, if such BDM is detected, analyzing the
time profile of the signal will allow direct measurements of $m_{\chi}$.

Events in Kamiokande and Super-K.—For BDM that also interact with electrons
with a cross section $\sigma_{\chi e}$, they can produce signals in neutrino
or DM experiments. The total event number $N_{\chi}$ induced by BDM with
$T_{\rm th}\leq T_{\chi}\leq T_{\rm max}$ within an exposure time $t_{0}\leq
t\leq t_{\rm exp}$ is given by

$N_{\chi}=N_{e}\sigma_{\chi e}\int^{T_{\rm max}}_{T_{\rm
th}}dT_{\chi}\int_{t_{0}}^{t_{\rm
exp}}dt~{}\epsilon\frac{d\Phi_{\chi}}{dT_{\chi}},$ (9)

with $N_{e}$ the total target number of electrons and $\epsilon$ the signal
efficiency. We consider the water Cherenkov experiments, Kamiokande and
Super-K, to calculate $N_{\chi}$ for BDM from LMC (by SN1987a) and from a SN
in GC. They have $N_{e}=(M_{T}/m_{\rm H_{2}O})N_{A}n_{e}$ with $M_{T}$ the
fiducial detector mass, $m_{\rm H_{2}O}$ the water molar mass, $N_{A}$ the
Avogadro constant and $n_{e}$ the electron number per water molecule. We take
$M_{T}=2.2$ and $22.5$ kton for Kamiokande and Super-K, respectively Hirata
_et al._ (1987); Abe  _et al._ (2016a), and set $(T_{\rm th},T_{\rm
max})=(5,100)\,{\rm MeV}$ for both. We make a conservative choice of taking
$\epsilon=50\%$, lower than the energy-dependent efficiency roughly ranging
from 50% to 75% reported for solar $\nu$ detection in Super-K Abe  _et al._
(2016a). For signal duration, we consider $t_{0}=10\,{\rm s}$ to approximately
exclude events produced by SN$\nu$, and let $t_{\rm exp}={\rm min}(t_{\rm
van},t_{\rm cut}=35~{}\rm{yrs})$ depending on $m_{\chi}$. For the LMC case,
the considered duration thus includes the running time of Kamiokande from 1987
to 1996 and Super-K after 1996 for heavier $m_{\chi}$. For the GC case, we
consider Super-K only. The main background for both comes from the solar and
atmospheric neutrinos for $T_{\chi}\lesssim 20\,{\rm MeV}$ and
$T_{\chi}\gtrsim 20\,{\rm MeV}$, respectively. We adopt values in Table XIV of
Ref. Abe  _et al._ (2016a) and the FLUKA simulation result in Ref. Battistoni
_et al._ (2005) to estimate the background.

Figure 4: BDM events in water Cherenkov detectors $N_{\chi}$ as a function of
$m_{\chi}$ for both the GC and LMC (SN1987a) cases. For LMC, events in
Kamiokande (hollow-green squares) and Super-K (hollow-orange triangles) are
shown separately. For the GC, only Super-K (red dots) is considered.
Background counts (dashed lines) are also shown for both cases.

Fig. 4 shows $N_{\chi}$ vs. $m_{\chi}$ resulting from the GC and LMC given
$\sigma_{\chi\nu}=\sigma_{\chi e}=10^{-35}$ cm2. We first discuss the GC case
where only Super-K is considered. The red-solid dots show that
$N_{\chi}\propto m_{\chi}^{-1}$ perfectly for $m_{\chi}\leq 25$ keV, which
corresponds to having $t_{\rm exp}=t_{\rm van}\leq t_{\rm cut}$. This is
because for smaller $m_{\chi}$, all BDM arrive at the detector before $t_{\rm
cut}$ so that $N_{\chi}$ is proportional to the amount of DM in the halo. For
heavier $m_{\chi}$, however, a larger part of BDM flux only arrives after
$t_{\rm cut}$ (see Fig. 3), leading to a faster decrease of $N_{\chi}$ with
increasing $m_{\chi}$. For the same reason, the background counts (red-dashed
curve) $N_{b}\simeq 526M_{T}t_{\rm exp}$ for $m_{\chi}\leq 25$ keV due to a
constant background rate of $\sim 526$ events per kton per year Abe  _et al._
(2016a). For $m_{\chi}>25$ keV where $t_{\rm exp}=t_{\rm cut}$ is applied,
$N_{b}$ stays constant.

For BDM associated with SN1987a in LMC, we plot $N_{\chi}$ in Kamiokande
(1987–1996) and Super-K (after 1996) by hollow green squares and orange
triangles separately. The behavior of $N_{\chi}(m_{\chi})$ in Kamiokande is
similar to that of the GC case, but falls off faster for large $m_{\chi}$ due
to the maximal exposure time of 9 yr only. The difference at small $m_{\chi}$
is mainly due to different detector fiducial mass $M_{T}$, geometric dilution
factor $1/R^{2}$, and the characteristic density $\rho_{s}$ of DM profiles. A
simple estimate gives $N_{\chi}^{\rm GC}/N_{\chi}^{\rm LMC}\sim(\rho_{s}^{\rm
NFW}/\rho_{s}^{\rm Hernquist})(R_{\rm LMC}^{2}/R_{\rm GC}^{2})(M_{T}^{\rm
SK}/M_{T}^{\rm Kamioka})\sim\mathcal{O}(10^{3})$ consistent with Fig. 4.
Super-K here only starts to accumulate events for $m_{\chi}\gtrsim 1.1$ keV
whose $t_{\rm van}>9~{}{\rm yr}$, and eventually dominates the contribution to
$N_{\chi}$ more than that from Kamiokande for larger $m_{\chi}$. For
comparison, we also plot the combined background numbers $N_{b}$ from both
detectors.

Sensitivity and constraint.—We use $N_{\chi}$ and $N_{b}$ derived above to
estimate the constraint and sensitivity on light DM, taking for simplicity

$n_{\sigma}=\frac{N_{\chi}}{\sqrt{N_{\chi}+N_{b}}}.$ (10)

The dependencies of sensitivity ($s$) versus $m_{\chi}$ are displayed in Fig.
5, where $s=\sqrt{\sigma_{\chi\nu}\sigma_{\chi e}}$ for the experimental limit
at $n_{\sigma}=1.64$ [90% confidence level (CL)] for SN1987a in LMC, and the
projected sensitivity at $n_{\sigma}=2.0$ for a SN in the GC. In order to
compare $s$ with existing constraints based exclusively on $\sigma_{\chi e}$
Cappiello and Beacom (2019); An  _et al._ (2018); Aprile  _et al._ (2019b, a);
Abramoff  _et al._ (2019); Agnese  _et al._ (2018), a model-dependent choice
relating $\sigma_{\chi\nu}$ and $\sigma_{\chi e}$ has to be made111Our bounds
are considerably better than those reported in Ref. Jho  _et al._ (2021); Das
and Sen (2021), which are not shown in Fig. 5. Also noted is that such a
comparison implicitly assumes a cross section that is independent of the
center-of-mass energy, which was adopted similarly in pioneer works for
cosmic-ray upscattered dark matter scenario Bringmann and Pospelov (2019); Ema
_et al._ (2019); Cappiello and Beacom (2019).. Under a generic description of
$r_{\chi}=\sigma_{\chi\nu}/\sigma_{\chi e}$, the specific case of $r_{\chi}=1$
was selected as illustration, with which the resulting bounds are superimposed
in Fig. 5. The Super-K constraints are derived from the average background
rates Abe  _et al._ (2016a) and statistical uncertainties. Time stability can
be inferred from the absence of anomalous time variations in the solar $\nu$
annual modulation analysis Yoo  _et al._ (2003); Cravens  _et al._ (2008).
Limits derived with BDM from SN1987a in LMC leads to orders of magnitude
improvement over existing bounds for $m_{\chi}<2\,{\rm keV}$ over a large
range of $r_{\chi}$. For instance, more stringent limits are derived at
$m_{\chi}\sim 10^{-6}$ MeV for $r_{\chi}>10^{-6}$. Moreover, a future SN in
the GC can improve the sensitivity by a factor of $\sim 30$ with Super-K,
since $s^{\rm GC}_{\chi}/s^{\rm LMC}_{\chi}\sim\sqrt{N^{\rm LMC}_{\chi}/N^{\rm
GC}_{\chi}}\sim\mathcal{O}(0.03$), allowing one to probe $s\lesssim 10^{-36}$
cm2 for $m_{\chi}\lesssim 10$ keV. For $m_{\chi}\lesssim 100\,{\rm keV}$, the
sensitivity curves follow $s\propto m_{\chi}^{1/2}$ simply due to
$N_{\chi}\propto m_{\chi}^{-1}$ (see Fig. 4).

Figure 5:  Sensitivity of BDM searches in $(m_{\chi}–s)$ plane.
$s=\sqrt{\sigma_{\chi\nu}\sigma_{\chi e}}$ for this work on the experimental
limits at $n_{\sigma}=1.64$ (90% CL) for LMC (green) and projected $s$ for the
GC at $n_{\sigma}=2$ (red). Current leading limits from cosmic-ray BDM
Cappiello and Beacom (2019), solar reflection An  _et al._ (2018), and DD
Aprile  _et al._ (2019b, a); Abramoff  _et al._ (2019); Agnese  _et al._
(2018) at $s=\sigma_{\chi e}$ are superimposed.

On the other hand, the sensitivities for $m_{\chi}\gtrsim 100\,{\rm keV}$
weaken considerably due to the reduced BDM that can arrive at the detector
within 35 years. Finally, we include an additional projection with Hyper-K for
the GC case (red-dashed curve). The analysis is similar to Super-K, with
fiducial mass and background rate scaled up by a factor of 10, which then
leads to another improvement of $\sim$ 2–3 over the Super-K result.

Summary and prospects.—We have examined the scenario of halo DM being boosted
by prompt SN$\nu$, and extracted a wealth of information from its TOF
measurements. The BDM events on Earth are characterized by unique timing
distributions, which vanish beyond $m_{\chi}$-dependent end points and are
independent of the interaction cross sections, while their peak positions
provide information on the SN locations and $m_{\chi}$.

A new constraint was derived on $s=\sqrt{\sigma_{\chi\nu}\sigma_{\chi e}}$
using Kamiokande and Super-K data on the SN1987a in LMC. Our results probe and
exclude new parameter space over a large range of $r_{\chi}$ and in particular
improve over the existing cosmic-ray BDM bounds for $m_{\chi}<100$ keV by 1–3
orders of magnitude at $\sigma_{\chi e}\sim\sigma_{\chi\nu}$. A future SN in
the GC can provide improved sensitivity by another factor of 30–100 with
Super-K or Hyper-K. The improvement over other probes Cappiello and Beacom
(2019); An  _et al._ (2018); Aprile  _et al._ (2019b, a); Abramoff  _et al._
(2019); Agnese  _et al._ (2018); Jho  _et al._ (2021); Das and Sen (2021) in
the sub-MeV mass range originates from the transient BDM flux arriving in a
short duration that can be calibrated by the detection of SN$\nu$, thereby
minimizing the background counts. The constraint and sensitivity of this work
were derived by a conservative analysis which stands on the BDM rates not
being larger than those of background. A detailed analysis that optimally
exploits the $m_{\chi}$-dependent TOF temporal profile or combines multiple
detectors is beyond the scope of this work but will further enhance the
sensitivities. Furthermore, most BDM arriving on Earth are within a small
solid angle relative to the SN direction for $m_{\chi}\lesssim\mathcal{O}({\rm
MeV})$ (see Fig. 2). Coupled with the good pointing capability for galactic SN
Beacom and Vogel (1999); Tomas  _et al._ (2003); Abe  _et al._ (2016b); Linzer
and Scholberg (2019), the angular information can be exploited to greatly
reduce the background.

Other effects such as the distortion of SN$\nu$ spectra, the recently proposed
SN$\nu$ echo Carpio  _et al._ (2022), and the impact of $\chi-\nu$ interaction
on SN$\nu$ emission have been neglected here. Estimations suggest that SN$\nu$
spectra be minimally affected for the parameter space examined. These effects
may be combined with the TOF profiles of SN$\nu$BDM to provide severe
constraints on specific phenomenological models relating $\sigma_{\chi e}$ and
$\sigma_{\chi\nu}$. With all the rich information, the next galactic SN will
offer new insights to the nature of DM. Furthermore, TOF analysis following SN
or other transient astrophysical events can be applied in a similar vein to
studies of other exotic physics interactions. A broad range of interesting
scenario will be explored in our future research.

Code availability.—We provide a Python package snorer sno , which can fully
reproduce the results in this Letter. The package is available on both GitHub
and PyPI. It offers numerous new features, such as including DM spikes, user-
specified SN locations in arbitrary distant galaxies, and an implementation of
any particle physics model. See its official page for further details.

###### Acknowledgements.

We thank Gang Guo, Kin-Wang Ng, Gianluca Petrillo, Yun-Tse Tsai for useful
discussions. Y.-H. L. acknowledges the Postdoctoral Scholar Program of the
Academia Sinica, Taiwan, during which the major part of this work was done, as
well as supports from National Science and Technology Council, the Ministry of
Education under Project No. NTU-112L104022, and the National Center for
Theoretical Sciences of Taiwan. M.-R. W. acknowledges supports from the
National Science and Technology Council, Taiwan under Grant No.
110-2112-M-001-050 and 111-2628-M-001-003-MY4, the Academia Sinica under
Project No. AS-CDA-109-M11, and Physics Division, National Center for
Theoretical Sciences of Taiwan. H. T.-K. W. acknowledges supports from the
National Science and Technology Council, Taiwan under Grant No.
106-2923-M-001-006-MY5. We also acknowledge the computing resources provided
by the Academia Sinica Grid-computing Center.

## References

  * Battaglieri  _et al._ (2017) M. Battaglieri _et al._ , in  _U.S. Cosmic Visions: New Ideas in Dark Matter_ (2017) arXiv:1707.04591 [hep-ph] .
  * Workman  _et al._ (2022) R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Aguilar  _et al._ (2015) M. Aguilar _et al._ (AMS), Phys. Rev. Lett. 115, 211101 (2015).
  * Akerib  _et al._ (2017a) D. S. Akerib _et al._ (LUX), Phys. Rev. Lett. 118, 021303 (2017a), arXiv:1608.07648 [astro-ph.CO] .
  * Ackermann  _et al._ (2017) M. Ackermann _et al._ (Fermi-LAT), Astrophys. J. 840, 43 (2017), arXiv:1704.03910 [astro-ph.HE] .
  * Akerib  _et al._ (2017b) D. S. Akerib _et al._ (LUX), Phys. Rev. Lett. 118, 251302 (2017b), arXiv:1705.03380 [astro-ph.CO] .
  * Aprile  _et al._ (2018) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Ambrosi  _et al._ (2017) G. Ambrosi _et al._ (DAMPE), Nature 552, 63 (2017), arXiv:1711.10981 [astro-ph.HE] .
  * Aprile  _et al._ (2019a) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 251801 (2019a), arXiv:1907.11485 [hep-ex] .
  * Aprile  _et al._ (2019b) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 241803 (2019b), arXiv:1907.12771 [hep-ex] .
  * Abramoff  _et al._ (2019) O. Abramoff _et al._ (SENSEI), Phys. Rev. Lett. 122, 161801 (2019), arXiv:1901.10478 [hep-ex] .
  * Agnese  _et al._ (2018) R. Agnese _et al._ (SuperCDMS), Phys. Rev. Lett. 121, 051301 (2018), [Erratum: Phys.Rev.Lett. 122, 069901 (2019)], arXiv:1804.10697 [hep-ex] .
  * Bringmann and Pospelov (2019) T. Bringmann and M. Pospelov, Phys. Rev. Lett. 122, 171801 (2019), arXiv:1810.10543 [hep-ph] .
  * Ema  _et al._ (2019) Y. Ema, F. Sala, and R. Sato, Phys. Rev. Lett. 122, 181802 (2019), arXiv:1811.00520 [hep-ph] .
  * Cappiello and Beacom (2019) C. V. Cappiello and J. F. Beacom, Phys. Rev. D 100, 103011 (2019), [Erratum: Phys.Rev.D 104, 069901 (2021)], arXiv:1906.11283 [hep-ph] .
  * Dent  _et al._ (2020) J. B. Dent, B. Dutta, J. L. Newstead, and I. M. Shoemaker, Phys. Rev. D 101, 116007 (2020), arXiv:1907.03782 [hep-ph] .
  * Wang  _et al._ (2020) W. Wang, L. Wu, J. M. Yang, H. Zhou, and B. Zhu, JHEP 12, 072 (2020), [Erratum: JHEP 02, 052 (2021)], arXiv:1912.09904 [hep-ph] .
  * Zhang (2022) Y. Zhang, PTEP 2022, 013B05 (2022), arXiv:2001.00948 [hep-ph] .
  * Guo  _et al._ (2020a) G. Guo, Y.-L. S. Tsai, and M.-R. Wu, JCAP 10, 049 (2020a), arXiv:2004.03161 [astro-ph.HE] .
  * Ge  _et al._ (2021) S.-F. Ge, J. Liu, Q. Yuan, and N. Zhou, Phys. Rev. Lett. 126, 091804 (2021), arXiv:2005.09480 [hep-ph] .
  * Cao  _et al._ (2021) Q.-H. Cao, R. Ding, and Q.-F. Xiang, Chin. Phys. C 45, 045002 (2021), arXiv:2006.12767 [hep-ph] .
  * Jho  _et al._ (2020) Y. Jho, J.-C. Park, S. C. Park, and P.-Y. Tseng, Phys. Lett. B 811, 135863 (2020), arXiv:2006.13910 [hep-ph] .
  * Cho  _et al._ (2020) W. Cho, K.-Y. Choi, and S. M. Yoo, Phys. Rev. D 102, 095010 (2020), arXiv:2007.04555 [hep-ph] .
  * Lei  _et al._ (2020) Z.-H. Lei, J. Tang, and B.-L. Zhang, (2020), arXiv:2008.07116 [hep-ph] .
  * Guo  _et al._ (2020b) G. Guo, Y.-L. S. Tsai, M.-R. Wu, and Q. Yuan, Phys. Rev. D 102, 103004 (2020b), arXiv:2008.12137 [astro-ph.HE] .
  * Xia  _et al._ (2021) C. Xia, Y.-H. Xu, and Y.-F. Zhou, Nucl. Phys. B 969, 115470 (2021), arXiv:2009.00353 [hep-ph] .
  * Dent  _et al._ (2021) J. B. Dent, B. Dutta, J. L. Newstead, I. M. Shoemaker, and N. T. Arellano, Phys. Rev. D 103, 095015 (2021), arXiv:2010.09749 [hep-ph] .
  * Ema  _et al._ (2021) Y. Ema, F. Sala, and R. Sato, SciPost Phys. 10, 072 (2021), arXiv:2011.01939 [hep-ph] .
  * Flambaum  _et al._ (2020) V. V. Flambaum, L. Su, L. Wu, and B. Zhu, (2020), arXiv:2012.09751 [hep-ph] .
  * Jho  _et al._ (2021) Y. Jho, J.-C. Park, S. C. Park, and P.-Y. Tseng, (2021), arXiv:2101.11262 [hep-ph] .
  * Das and Sen (2021) A. Das and M. Sen, Phys. Rev. D 104, 075029 (2021), arXiv:2104.00027 [hep-ph] .
  * Bell  _et al._ (2021) N. F. Bell, J. B. Dent, B. Dutta, S. Ghosh, J. Kumar, J. L. Newstead, and I. M. Shoemaker, Phys. Rev. D 104, 076020 (2021), arXiv:2108.00583 [hep-ph] .
  * Chao  _et al._ (2021) W. Chao, T. Li, and J. Liao, (2021), arXiv:2108.05608 [hep-ph] .
  * Ghosh  _et al._ (2022) D. Ghosh, A. Guha, and D. Sachdeva, Phys. Rev. D 105, 103029 (2022), arXiv:2110.00025 [hep-ph] .
  * Feng  _et al._ (2022) J.-C. Feng, X.-W. Kang, C.-T. Lu, Y.-L. S. Tsai, and F.-S. Zhang, JHEP 04, 080 (2022), arXiv:2110.08863 [hep-ph] .
  * Wang  _et al._ (2021) W. Wang, L. Wu, W.-N. Yang, and B. Zhu, (2021), arXiv:2111.04000 [hep-ph] .
  * Xia  _et al._ (2022) C. Xia, Y.-H. Xu, and Y.-F. Zhou, JCAP 02, 028 (2022), arXiv:2111.05559 [hep-ph] .
  * Wang  _et al._ (2022) J.-W. Wang, A. Granelli, and P. Ullio, Phys. Rev. Lett. 128, 221104 (2022), arXiv:2111.13644 [astro-ph.HE] .
  * Cui  _et al._ (2022) X. Cui _et al._ (PandaX-II), Phys. Rev. Lett. 128, 171801 (2022), arXiv:2112.08957 [hep-ex] .
  * Xu  _et al._ (2022) R. Xu _et al._ (CDEX), (2022), arXiv:2201.01704 [hep-ex] .
  * Granelli  _et al._ (2022) A. Granelli, P. Ullio, and J.-W. Wang, (2022), arXiv:2202.07598 [astro-ph.HE] .
  * Arnett and Rosner (1987) W. D. Arnett and J. L. Rosner, Phys. Rev. Lett. 58, 1906 (1987).
  * Kolb  _et al._ (1987) E. W. Kolb, A. J. Stebbins, and M. S. Turner, Phys. Rev. D 35, 3598 (1987).
  * Holdom (1986) B. Holdom, Phys. Lett. B 166, 196 (1986).
  * Argüelles  _et al._ (2017) C. A. Argüelles, A. Kheirandish, and A. C. Vincent, Phys. Rev. Lett. 119, 201801 (2017), arXiv:1703.00451 [hep-ph] .
  * Chang  _et al._ (2018) J. H. Chang, R. Essig, and S. D. McDermott, JHEP 09, 051 (2018), arXiv:1803.00993 [hep-ph] .
  * Murase and Shoemaker (2019) K. Murase and I. M. Shoemaker, Phys. Rev. Lett. 123, 241102 (2019), arXiv:1903.08607 [hep-ph] .
  * Lin and Lin (2021) G.-L. Lin and Y.-H. Lin, Phys. Rev. D 104, 063021 (2021), arXiv:2102.11151 [hep-ph] .
  * Croon  _et al._ (2021) D. Croon, G. Elor, R. K. Leane, and S. D. McDermott, JHEP 01, 107 (2021), arXiv:2006.13942 [hep-ph] .
  * Escudero  _et al._ (2019) M. Escudero, D. Hooper, G. Krnjaic, and M. Pierre, JHEP 03, 071 (2019), arXiv:1901.02010 [hep-ph] .
  * Foldenauer (2019) P. Foldenauer, Phys. Rev. D 99, 035007 (2019), arXiv:1808.03647 [hep-ph] .
  * Duan  _et al._ (2006) H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys. Rev. D 74, 105014 (2006), arXiv:astro-ph/0606616 .
  * Bertone  _et al._ (2005) G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005), arXiv:hep-ph/0404175 .
  * Erkal  _et al._ (2019) D. Erkal _et al._ , Mon. Not. Roy. Astro. 487, 2685 (2019).
  * Hirata  _et al._ (1987) K. Hirata _et al._ , Phys. Rev. Lett. 58, 1490 (1987).
  * Abe  _et al._ (2016a) K. Abe _et al._ (Super-Kamiokande), Phys. Rev. D 94, 052010 (2016a), arXiv:1606.07538 [hep-ex] .
  * Battistoni  _et al._ (2005) G. Battistoni, A. Ferrari, T. Montaruli, and P. R. Sala, Astropart. Phys. 23, 526 (2005).
  * An  _et al._ (2018) H. An, M. Pospelov, J. Pradler, and A. Ritz, Phys. Rev. Lett. 120, 141801 (2018), [Erratum: Phys.Rev.Lett. 121, 259903 (2018)], arXiv:1708.03642 [hep-ph] .
  * Yoo  _et al._ (2003) J. Yoo _et al._ (Super-Kamiokande), Phys. Rev. D 68, 092002 (2003), arXiv:hep-ex/0307070 .
  * Cravens  _et al._ (2008) J. P. Cravens _et al._ (Super-Kamiokande), Phys. Rev. D 78, 032002 (2008), arXiv:0803.4312 [hep-ex] .
  * Beacom and Vogel (1999) J. F. Beacom and P. Vogel, Phys. Rev. D 60, 033007 (1999), arXiv:astro-ph/9811350 .
  * Tomas  _et al._ (2003) R. Tomas, D. Semikoz, G. G. Raffelt, M. Kachelriess, and A. S. Dighe, Phys. Rev. D 68, 093013 (2003), arXiv:hep-ph/0307050 .
  * Abe  _et al._ (2016b) K. Abe _et al._ (Super-Kamiokande), Astropart. Phys. 81, 39 (2016b), arXiv:1601.04778 [astro-ph.HE] .
  * Linzer and Scholberg (2019) N. B. Linzer and K. Scholberg, Phys. Rev. D 100, 103005 (2019), arXiv:1909.03151 [astro-ph.IM] .
  * Carpio  _et al._ (2022) J. A. Carpio, A. Kheirandish, and K. Murase, (2022), arXiv:2204.09650 [hep-ph] .
  * (66) GitHub: https://github.com/yenhsunlin/snorer and PyPI: https://pypi.org/project/snorer/.