# GCNs-Net: A Graph Convolutional Neural Network Approach for Decoding Time-
resolved EEG Motor Imagery Signals

Xiangmin Lun, Shuyue Jia, , Yimin Hou, Yan Shi, Yang Li, , Hanrui Yang, Shu
Zhang, and Jinglei Lv  This work was supported by the National Natural Science
Foundation of China under Grant 31772059. Xiangmin Lun, Yimin Hou, Yan Shi,
and Hanrui Yang are with the School of Automation Engineering, Northeast
Electric Power University, Jilin, China. Xiangmin Lun is also with the College
of Mechanical and Electric Engineering, Changchun University of Science and
Technology, Changchun, China. Shuyue Jia, the corresponding author, is with
the School of Computer Science, Northeast Electric Power University, Jilin,
China. (E-mail: shuyuej@ieee.org). Yang Li is with the School of Electrical
Engineering, Northeast Electric Power University, Jilin, China. Shu Zhang is
with the School of Computer Science, Northwestern Polytechnical University,
Xi’an, China. Jinglei Lv is with the Sydney Imaging & School of Biomedical
Engineering, The University of Sydney, Sydney, Australia.

###### Abstract

Towards developing effective and efficient brain-computer interface (BCI)
systems, precise decoding of brain activity measured by electroencephalogram
(EEG), is highly demanded. Traditional works classify EEG signals without
considering the topological relationship among electrodes. However,
neuroscience research has increasingly emphasized network patterns of brain
dynamics. Thus, the Euclidean structure of electrodes might not adequately
reflect the interaction between signals. To fill the gap, a novel deep
learning framework based on the graph convolutional neural networks (GCNs) was
presented to enhance the decoding performance of raw EEG signals during
different types of motor imagery (MI) tasks while cooperating with the
functional topological relationship of electrodes. Based on the absolute
Pearson’s matrix of overall signals, the graph Laplacian of EEG electrodes was
built up. The GCNs-Net constructed by graph convolutional layers learns the
generalized features. The followed pooling layers reduce dimensionality, and
the fully-connected softmax layer derives the final prediction. The introduced
approach has been shown to converge for both personalized and group-wise
predictions. It has achieved the highest averaged accuracy, 93.056% and 88.57%
(PhysioNet Dataset), 96.24% and 80.89% (High Gamma Dataset), at the subject
and group level, respectively, compared with existing studies, which suggests
adaptability and robustness to individual variability. Moreover, the
performance was stably reproducible among repetitive experiments for cross-
validation. To conclude, the GCNs-Net filters EEG signals based on the
functional topological relationship, which manages to decode relevant features
for brain motor imagery.

###### Index Terms:

Electroencephalography (EEG), Motor Imagery (MI), Deep Learning (DL), Graph
Convolutional Neural Networks (GCNs), Brain-Computer Interface (BCI)

##  I Introduction

Recently, the brain-computer interface (BCI) has become one of the hottest
research topics for broad applications in the field of therapeutic and medical
engineering[1]. It refers to the establishment of an innovative technology
that exchanges information directly between the brain and the surroundings,
which does not rely on traditional methods such as human muscle tissue or
peripheral nerves. BCI systems decode brain activity patterns to manipulate
assistant devices, such as wheelchairs and artificial limbs[2].
Electroencephalogram (EEG) is extensively applied because of its high temporal
resolution, noninvasiveness, and portability. The principle of EEG is to
record spontaneous, event-related, and stimulus-evoked electrical signals of
the brain on time scales, which reveals variations for different brain
activities[3]. EEG decodes discriminable brain patterns while carrying out
different types of actual movement or imagery[4, 5]. Motor imagery (MI) based
EEG mentally simulates multiple motor motions, such as imagining hand or foot
movements. Controlling machines via only the MI without physical movements of
the body is one of the elemental jobs in the area of BCI[6]. To realize such
BCI systems, accurate classification of MI brain activity is of great essence.
Although previous studies have shown promising performance, there is still
space to improve the classification accuracy towards building effective and
efficient BCI applications. For instance, the adaptability and robustness to
individual variability remain among the challenges of setting up an EEG MI-
based wheelchair. Traditional approaches did not consider the topological
relationship of electrodes while decoding EEG signals. However, a growing
number of neuroscience research has emphasized the brain network dynamics[7,
8, 9]. Thus, the interaction between signals might not be adequately reflected
and represented via the Euclidean structure of EEG electrodes. To address the
concern, the Graph Convolutional Neural Networks (GCNs) were introduced to
decode EEG signals, promoting the classification performance by cooperating
with the functional topological relationship of EEG electrodes and
implementing the Convolutional Neural Networks (CNNs) on graphs.

###  I-A Literature Survey

Traditional works manually designed features from EEG signals, e.g., via the
analytic intrinsic mode functions, or wavelet transform, and then employed
machine learning-based approaches to classify features[10, 11, 12, 13, 14].
Recently, deep learning (DL) has accomplished superhuman performance across
multiple domains[15]. The DL-based methods learn the underlying features from
signals, which alleviates the need for hands-on feature engineering. The CNNs
have been broadly employed to classify the Euclidean-structured signals on
account of its ability to learn informative features from local receptive
fields. References [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]
implemented the CNN-based approaches to address the challenge of EEG task
classification. Reference [18] introduced an innovative approach by combining
the Scout EEG Source Imaging (ESI) and CNNs to decode EEG tasks, which
achieved competitive results, 94.5% maximum accuracy for 10 subjects, and
92.5% for 14 subjects, on the PhysioNet Dataset[28]. Reference [29] presented
the CNNs-LSTM method (cascade convolutional recurrent neural network), and
obtained 98.31% averaged accuracy on the PhysioNet Dataset. Reference [20]
applied one-dimensional convolutional filters to learn temporal and spatial
features, and it reached 80.38%, 69.82%, and 58.58% accuracies on the
PhysioNet Dataset with two, three, and four MI tasks. References [22, 19, 24,
25] utilized variants of CNNs to decode EEG signals from the BCI Competition
IV-2a Dataset[30], and achieved 73.70% (CNNs), 75.70% (MCNNs), 79.90% (DFFN),
and 83.00% (CNNs-LSTMs) accuracy, respectively. References [26, 27, 19]
obtained 92.50% accuracy (CNNs), 93.70% (CP-MixedNet), and 95.4% (MCNNs) at
the subject level on the High Gamma Dataset[31]. Furthermore, some recent
approaches have focused on the subject-independent EEG MI decoding, which is
considered our future work[32, 33, 34]. Reference [32] has accomplished 59.10%
averaged accuracy via a convolutional recurrent attention model (CRAM), an
attention-based CNNs approach, on the BCI Competition IV-2a Dataset. The
following two references used the graph structure to represent EEG topologies
and then applied CNN-based methods to decode EEG tasks. References [33]
applied the Graph-based Convolutional Recurrent Attention Model (G-CRAM),
where the graph representation was regarded as the input to the CNNs and
achieved 74.71% and 60.11% accuracy on the PhysioNet Dataset and BCI
Competition IV-2a Dataset, respectively. Reference [34] introduced a Graph-
based Hierarchical Attention M odel (G-HAM) and obtained 76.36% accuracy on
the PhysioNet Dataset. Although the performances of the above CNN-based models
were encouraging, there was still space to promote the classification accuracy
to build a robust and reliable BCI system. The reasons why we applied the GCNs
were as follows. The typical CNNs cannot process non-Euclidean structured data
directly because the discrete convolution cannot keep translation invariance
on the non-Euclidean signals. However, the GCNs can directly extract features
from the non-Euclidean data and process graph-structured signals since the
GCNs consider the relationship properties (e.g., correlations) between
nodes[35]. Through the end-to-end GCNs-Net approach on two EEG MI benchmarks,
it has achieved dominant performances, 98.72% accuracy on the PhysioNet
Dataset, and 96.24% on the High Gamma Dataset, for the EEG MI decoding, which
were far ahead than the results produced by the CNN-based approaches.
Moreover, more and more neuroscience research has suggested that the
topological information can promote the analysis of brain network dynamics[7,
8, 9]. Although the Euclidean distance is one of the similarity measurements,
it might be superior to decode EEG signals from the non-Euclidean perspective
by taking consideration of the functional topological relationship of
electrodes (e.g., the correlation and degree properties between electrodes) to
enhance the decoding performance of EEG tasks.

Considering the topological relationship of EEG electrodes, the graph in the
non-Euclidean space was put forward[36, 37, 38, 39]. Researchers have
investigated CNNs with the graph theory, intending to apply the convolutional
operation on graphs. Two strategies were presented to define the convolutions,
i.e., either from the spatial domain or from a spectral domain. At first,
spatial GCNs were proposed[40, 41, 42, 43]. However, references [44, 35]
indicated that it faced the challenge of matching the local neighborhoods, and
there was no convincing mathematical definition supporting such operations. On
the other hand, the spectral method provided a well-defined localized operator
on graphs[45]. Thus, references [44, 45] proposed an innovative approach of
the GCNs by combining the CNNs with the spectral graph theory. Noticeably,
reference [44] first performed the GCNs from the spectral perspective. A few
works have applied the above model to decode EEG tasks, primarily in EEG-based
emotion recognition[36, 37, 38, 39]. In detail, reference [36] combined the
GCNs with the broad learning system and put forward the Graph Convolutional
Broad Network, which achieved 94.24% accuracy on the SJTU emotion EEG Dataset.
Reference [37] presented phase-locking value-based GCNs. Reference [38]
proposed dynamical GCNs that can dynamically learn the topological
relationship of EEG electrodes during training. Reference [39] improved the
above method via a broad learning system. Nevertheless, to date, the GCNs have
not been studied in the field of EEG MI. In this work, a novel structure of
the GCNs was introduced to decode EEG MI signals. First of all, based on the
absolute Pearson’s matrix of overall signals, the graph Laplacian was built up
to represent the topological relationship of EEG electrodes. Besides, the
GCNs-Net built on graph convolutional layers learned the generalized features.
The followed pooling layers reduced dimensionality. And the fully-connected
(FC) softmax layer derived the final prediction. Furthermore, the Chebyshev
polynomial was applied to approximate graph convolutional filters, which
significantly promoted the computation efficiency. Last but not least, the
GCNs-Net decoded time-resolved EEG MI signals, which paved the road towards
effective and efficient BCI applications.

Figure 1: The system framework composed of i) 64-channel raw EEG signals’
acquisition, ii) correlation analysis for graph weights and degrees by
presenting absolute PCC matrix, adjacency matrix, and graph Laplacian, iii)
graph representation of EEG electrodes, and iv) novel deep learning framework
of the GCNs (the GCNs-Net).

###  I-B Contributions

The main contributions of this work were summarized as follows.
(i) A novel structure of the GCNs was first introduced to detect four-class MI
intentions while cooperating with the functional topological relationship of
EEG electrodes.
(ii) The individual and group-wise performances of the GCNs-Net on two
benchmark datasets of EEG MI outperformed the existing studies, which
validated that the method can decode relevant features for brain motor
imagination.
(iii) The introduced GCNs-Net framework can be easily transferred and
implemented for other MI related tasks, and potentially for other EEG BCI
tasks.

##  II Materials and Methods

###  II-A Overview

The framework of this work was shown schematically in Figure 1.
(i) 64-channel raw EEG signals were acquired as one of the inputs for the
GCNs-Net.
(ii) The PCC matrix, absolute PCC matrix, adjacency matrix, and graph
Laplacian were introduced to represent the correlations between electrodes.
(iii) The graph representation, another input for the GCNs-Net, was
represented by the graph Laplacian.
(iv) The GCNs-Net was applied to decode EEG MI signals. Where, $\rm N$ denotes
the number of electrodes, and $l$ denotes the $l^{th}$ graph convolutional
layer.

###  II-B Dataset Description

In this work, two benchmark datasets were used to evaluate the effectiveness
and robustness of our presented method. The PhysioNet Dataset included over
1,500 EEG records from 109 subjects[28]. There were 64 electrodes based on the
international 10-10 system. Each subject performed 84 trials (3 runs $\times$
7 trials $\times$ 4 tasks). 160 Hz sampling rate and 4 seconds’ signals, i.e.,
640 time points per trail, were utilized considering the experiment
duration[18]. Four MI tasks were termed as L (Image Left Fist), R (Image Right
Fist), B (Image Both Fists), and F (Image Both Feet), respectively.

Collected from 14 subjects, the Public High Gamma Dataset performed four EEG
tasks, left-hand movement, right-hand movement, both feet movement, and rest.
44 electrodes’ data and 0-125 Hz frequency were applied in our experiments,
and the dataset was resampled to 250 Hz[22]. Each subject performed
approximately 880 trials for training and 160 trials for testing.

Previous studies applied segments (windows) of time points as samples[46, 47].
However, neuroscience research indicates that the brain is one of the most
complicated systems, and its state at every moment is changing. Time-resolved
signals can represent the condition of the brain at the instant moment, which
reflects the network patterns of brain dynamics. Therefore, every time point
was recognized as a sample to map the brain states with the corresponding MI
task in this work. Compared with applying time-window singles as samples, the
method was time-resolved, which was superior to pave the road towards
developing real-time and efficient EEG MI applications. In this paper,
multiple subjects (S1$\sim$S20, S1$\sim$S50, S1$\sim$S100 from the PhysioNet
Dataset, and 14 subjects from the High Gamma Dataset) were picked up to train
and evaluate the GCNs-Net. For all experiments, 90% of the data was randomly
chosen as the training set, and the left 10% served as the test set[18, 29].

###  II-C Graph Preliminary

####  II-C1 Graph Representation

An undirected and weighted graph is represented by
$\mathcal{G}=\\{{\mathcal{V},\mathcal{E},\mathcal{A}\\}}$, in which
$\mathcal{V}$ denotes a set of nodes with the number of
$|\mathcal{V}|=\emph{N}$, $\mathcal{E}$ denotes a set of edges connecting
nodes, and
$\mathcal{A}\in\mathbb{R}\textsuperscript{\emph{N}$\times$\emph{N}}$ is a
weighted adjacency matrix representing correlations between two nodes.

To present the degree matrix of a graph, the scale of the graph weights is
analysed regardless of the polar relevance, i.e., either the correlations are
positive or negative. So the Absolute Pearson Correlation Coefficient (PCC)
Matrix $|P|\in[0,1]$ was introduced, which was the absolute of the PCC Matrix
$P$ mapping the linear correlations between signals. $\mathcal{A}$ was
represented as $\mathcal{A}=|P|-\emph{I}$, where I was an identity matrix. The
PCC matrix, absolute PCC matrix, adjacency matrix, and graph Laplacian for 20
and 100 subjects are given in Figure 2.

(a) PCC Matrix for 20 Subjects

(b) Absolute PCC Matrix for 20 Subjects

(c) Adjacency Matrix for 20 Subjects

(d) Graph Laplacian for 20 Subjects

(e) PCC Matrix for 100 Subjects

(f) Absolute PCC Matrix for 100 Subjects

(g) Adjacency Matrix for 100 Subjects

(h) Graph Laplacian for 100 Subjects

Figure 2: The PCC, absolute PCC, adjacency matrices, and Graph Laplacian for
20 and 100 subjects, respectively, from the PhysioNet Dataset.

Figure 2 and Figure 2 demonstrate that there are no apparent distinctions on
correlations of various sizes of dataset. The reason for this phenomenon is
that all the signals of participants are utilized to compute the PCC Matrix.
Furthermore, the degree matrix $D$ is obtained, which is a diagonal matrix,
and the $\it{i}$-th diagonal element can be calculated by:
$D_{ii}=\sum\limits_{j=1}^{N}\mathcal{A}_{ij}$. Finally, combinatorial
Laplacian $L\in\mathbb{R}\textsuperscript{\emph{N}$\times$\emph{N}}$ is
represented as $L=D-\mathcal{A}$. The graph Laplacians for 20 and 100 subjects
are represented in Figure 2 and Figure 2, respectively. And the normalized
graph Laplacian is as follows: $L=I_{N}-D^{-1/2}\mathcal{A}D^{-1/2}$ to
represent the correlations between nodes.

####  II-C2 Spectral Graph Filtering

As $L$ is a real symmetric positive semidefinite matrix, its set of
eigenvectors, called the graph Fourier modes,
$\left\\{u_{l}\right\\}_{l=0}^{N-1}\in\mathbb{R}^{N}$ are complete and
orthonormal. And the associated eigenvalues, as known as the graph Fourier
Frequencies, $\left\\{\lambda_{l}\right\\}_{l=0}^{N-1}$ are ordered and real
nonnegative. The Fourier basis
$U=\left[u_{0},\ldots,u_{N-1}\right]\in\mathbb{R}^{N\times N}$ decomposes the
graph Laplacian, i.e., $L=U\Lambda U^{T}$ where
$\Lambda=\operatorname{diag}\left(\left[\lambda_{0},\ldots,\lambda_{N-1}\right]\right)\in\mathbb{R}^{N\times
N}$. The signal $x\in\mathbb{R}^{N}$ transformed by graph Fourier is
represented as $\hat{x}=\hat{U}^{T}x\in\mathbb{R}^{N}$, and its inverse graph
Fourier transform as $x=U\hat{x}$[45]. It has projected the input signals to
an orthonormal space where the basis is formed by the eigenvectors of
normalized graph Laplacian[48].

The convolution on graphs $*_{\mathcal{G}}$ is defined below:

$x*_{\mathcal{G}}g=U\left(\left(U^{T}x\right)\odot\left(U^{T}g\right)\right)$
(1)

In which $\odot$ represents element-wise Hadamard product. Followed by a non-
parametric filter $g_{\theta}$, where the vector of Fourier coefficients
$\theta\in\mathbb{R}^{N}$, and
$g_{\theta}(\Lambda)=\operatorname{diag}(\theta)$. The convolution operation
implemented for the GCNs is shown in the following.

$y=g_{\theta}(L)x=g_{\theta}\left(U\Lambda
U^{T}\right)x=Ug_{\theta}(\Lambda)U^{T}x$ (2)

TABLE I: Implementation details of the GCNs-Net for the PhysioNet Dataset. Layer | Type | Maps | Size | Edges |  | Polynomial
---
Orders
| Pooling
---
Size
Activation | Weights | Bias
Softmax | Fully-connected | - | O | - | - | - | Softmax |  $\frac{\rm N}{64}\times\frac{\rm N}{64}\times$F${}_{6}\times$O | O
Flatten | Flatten | - |  $\frac{\rm N}{64}\times\frac{\rm N}{64}\times$F6 | - | - | - | - | - | -
P6 | Max-pooling | F6 | $\frac{\rm N}{32}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{32}-1}i$ | - | 2 | - | - | -
C6 | Convolution | F6 | $\frac{\rm N}{32}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{32}-1}i$ | K | - | Softplus | F${}_{5}\times$F${}_{6}\times$K |  $\frac{\rm N}{32}\times$F6
P5 | Max-pooling | F5 | $\frac{\rm N}{16}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{16}-1}i$ | - | 2 | - | - | -
C5 | Convolution | F5 | $\frac{\rm N}{16}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{16}-1}i$ | K | - | Softplus | F${}_{4}\times$F${}_{5}\times$K |  $\frac{\rm N}{16}\times$F5
P4 | Max-pooling | F4 | $\frac{\rm N}{8}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{8}-1}i$ | - | 2 | - | - | -
C4 | Convolution | F4 | $\frac{\rm N}{8}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{8}-1}i$ | K | - | Softplus | F${}_{3}\times$F${}_{4}\times$K |  $\frac{\rm N}{8}\times$F4
P3 | Max-pooling | F3 | $\frac{\rm N}{4}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{4}-1}i$ | - | 2 | - | - | -
C3 | Convolution | F3 | $\frac{\rm N}{4}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{4}-1}i$ | K | - | Softplus | F${}_{2}\times$F${}_{3}\times$K |  $\frac{\rm N}{4}\times$F3
P2 | Max-pooling | F2 | $\frac{\rm N}{2}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{2}-1}i$ | - | 2 | - | - | -
C2 | Convolution | F2 | $\frac{\rm N}{2}$ | $\sum\nolimits_{i=1}^{\frac{\rm N}{2}-1}i$ | K | - | Softplus | F${}_{1}\times$F${}_{2}\times$K |  $\frac{\rm N}{2}\times$F2
P1 | Max-pooling | F1 | $\rm N$ | $\sum\nolimits_{i=1}^{N-1}i$ | - | 2 | - | - | -
C1 | Convolution | F1 | $\rm N$ | $\sum\nolimits_{i=1}^{N-1}i$ | K | - | Softplus | 1$\times$F${}_{1}\times$K | N$\times$F1
Input | Input | 1 | $\rm N$ | $\sum\nolimits_{i=1}^{N-1}i$ | - | - | - | - | -

The difference of spectral graph convolution lies in the choice of filter
$g_{\theta}$. The non-parametric filters are not localized in space, and their
computation complexity is too expensive. However, a polynomial filter can
address the problem[35]. One such method conventionally utilized to
approximate filters is the Chebyshev polynomials[49]. As the Chebyshev
polynomials approximation of the diagonal matrix of the eigenvalues, it can be
parametrized as truncated expansion of $k-1$ order:

$g_{\theta}(\Lambda)=\sum_{k=0}^{K-1}\theta_{k}T_{k}(\tilde{\Lambda})$ (3)

In which parameter $\theta\in\mathbb{R}^{K}$ is a set of Chebyshev
coefficients, $T_{k}(\tilde{\Lambda})\in\mathbb{R}^{K}$ is the $k^{th}$ order
Chebyshev polynomial evaluated at
$\tilde{\Lambda}=2\Lambda/\lambda_{max}-I_{n}$, and $I_{n}$ is a diagonal
matrix of the scaled eigenvalues.

Then, the signal $x$ is convolutioned by the defined filter $g_{\theta}$ as
follows:

$y=g_{\theta}(L)x=U\sum_{k=0}^{K-1}\theta_{k}T_{k}(\tilde{\Lambda})U^{T}x=\sum_{k=0}^{K-1}\theta_{k}T_{k}(\tilde{L})x$
(4)

Where $T_{k}(\tilde{L})\in\mathbb{R}^{n\times n}$ is the Chebyshev polynomial
of order $k$ evaluated at the scaled Laplacian
$\tilde{L}=2L/\lambda_{max}-I_{n}$. Let
$\bar{x}_{k}=T_{k}(\tilde{L})x\in\mathbb{R}^{n}$, the recursive relation is
utilized to compute $\bar{x}_{k}=2\tilde{L}\bar{x}_{k-1}-\bar{x}_{k-2}$ with
$\bar{x}_{0}=x$, and $\bar{x}_{1}=\tilde{L}x$. Another reason why Chebyshev
polynomial is applied to approximate convolutional filters is that it
implicitly avoids computations for graph Fourier basis, thus reducing the
computation complexity from $O(N^{2})$ to $O(KN)$.

####  II-C3 Graph Coarsening and Fast Pooling

Compared with the pooling operations of the regular CNNs, on graphs, it
involves nodes clustering and one-dimensional pooling. To carry out pooling to
diminish the dimensionality, the Graclus multilevel clustering algorithm is
performed[50]. A greedy algorithm was employed to measure the consecutive
coarser of the graph, and minimize the objective of spectral clustering[51].
(i) At each level, multiple numbers of coarser graphs are given. (ii) It picks
an unmarked node $i$, and matches with its unmarked neighborhood $j$, which
needs to maximize the local normalized cut
$W_{ij}\left(1/d_{i}+1/d_{j}\right)$. (iii) It will mark the two matched
nodes, and the sum of their weights will be the coarsened weight. (iv) All the
nodes will undergo the same procedure. At the coarsest level, the nodes will
be arbitrarily ordered. Next, the ordered nodes will be propagated to the most
finest level. Then, the graph signal is pooled in a one-dimensional pooling
manner[35]. This algorithm cuts the number of nodes by two between two levels.

###  II-D Model Initialization

A novel structure of the GCNs was introduced to classify the MI tasks. Based
on the absolute Pearson’s matrix of overall signals, the graph Laplacian was
built to represent the topological relationship of EEG electrodes. The graph
convolutional layers learned the generalized features. Built on a maximum of
$\rm\log_{2}{N}$ graph pooling layers regarding $\rm N$ EEG channels, the
pooling operation reduced the dimensionality, and the FC softmax layer derived
the final prediction. With regard to 64-channel EEG systems, the maximum
number of pooling layers was six. The implementation details were listed in
Table I, where $\rm N$ denotes the input size of the EEG signals, ${\rm
F}_{i}\in\left[{\rm F}_{1},{\rm F}_{2},{\rm F}_{3},{\rm F}_{4},{\rm
F}_{5},{\rm F}_{6}\right]$ donates the number of filters at each graph
convolutional layer, K denotes the polynomial order for filters, and O denotes
the number of MI tasks.

The hyperparameters in our work (e.g., learning rate, Dropout rate, and L2
regularization coefficient) were mainly empirically chosen during all the
experiments, which were not task-oriented tuned. The network parameters,
weights and biases, were updated via the Stochastic Gradient Descent (SGD)
algorithm implemented by the Adam optimizer (0.01 learning rate). Biases were
applied for every node of the graph. The batch size was 1,024 to maximize the
usage of GPU resources. For FC layers in Section III-A, a 50% dropout rate was
applied. Batch normalization (BN) has been employed for graph convolutions.
The non-linear Softplus (Smooth Rectified Linear Unit) activation function was
applied to the graph convolutional layers and FC layers to prevent gradient
vanishing.

$f(x)=\log\left(1+e^{x}\right)$ (5)

Softmax function was utilized to derive the final prediction, which denoted as
$\hat{y}\in[y^{1},y^{2},y^{3},y^{4}]$.

$\hat{y}=\operatorname{argmax}\left(\frac{e^{y^{i}}}{\sum_{i=1}^{4}e^{y^{i}}}\right)$
(6)

Where $\hat{y}$ was the output class. The cross-entropy with the L2
regularization was used as the loss function. Meanwhile, $1\times 10^{-6}$ was
applied as the coefficient $\lambda$ of the L2 regularization.

$\text{loss}=-\sum_{i=1}^{4}y_{i}\log\left(\hat{y}_{i}\right)+\lambda\|\sum_{j=1}^{n}{w_{j}^{2}+b_{j}^{2}}\|$
(7)

where $y_{i}$ and $n$ were the label and the number of input samples.

TABLE II: Performance comparison of the GCNs-Net

Model |  | Num of
---
Conv
Layers
| Num of
---
Pooling
Layers
Filter’s Account | Model Framework |  | GAA w.r.t.
---
$1^{st}order$
| GAA w.r.t.
---
$2^{nd}order$
| GAA w.r.t.
---
$3^{rd}order$
| GAA w.r.t.
---
$4^{th}order$
| GAA w.r.t.
---
$5^{th}order$
1 | 1 | 1 | 16 | C-P-S | 55.63% | 55.30% | 56.70% | 56.60% | 56.04%
2 | 2 | 1 | 16, 32 | C-C-P-S | 57.94% | 61.90% | 63.17% | 63.16% | 63.37%
3 | 2 | 2 | 16, 32 | (C-P)$\times$2-S | 60.04% | 62.32% | 62.55% | 62.96% | 62.08%
4 | 3 | 1 | 16, 32, 64 | C-C-C-P-S | 58.07% | 69.18% | 69.86% | 71.17% | 70.98%
5 | 3 | 2 | 16, 32, 64 | C-(C-P)$\times$2-S | 61.65% | 69.73% | 70.19% | 71.06% | 71.45%
6 | 3 | 3 | 16, 32, 64 | (C-P)$\times$3-S | 65.03% | 70.50% | 69.12% | 70.36% | 71.30%
7 | 4 | 2 | 16, 32, 64, 128 | C-C-(C-P)$\times$2-S | 63.27% | 77.09% | 77.04% | 78.42% | 78.15%
8 | 4 | 2 | 16, 32, 64, 128 | (C-C-P)$\times$2-S | 63.69% | 77.59% | 77.32% | 79.28% | 77.43%
9 | 4 | 3 | 16, 32, 64, 128 | C-(C-P)$\times$3-S | 67.50% | 77.63% | 77.67% | 79.60% | 78.36%
10 | 4 | 4 | 16, 32, 64, 128 | (C-P)$\times$4-S | 71.22% | 77.61% | 78.50% | 78.22% | 78.26%
11 | 5 | 3 | 16, 32, 64, 128, 256 | C-C-(C-P)$\times$3-S | 70.11% | 83.06% | 82.78% | 84.29% | 84.19%
12 | 5 | 3 | 16, 32, 64, 128, 256 | (C-C-P)$\times$2-C-P-S | 70.12% | 83.05% | 82.49% | 84.26% | 83.45%
13 | 5 | 4 | 16, 32, 64, 128, 256 | C-(C-P)$\times$4-S | 75.93% | 84.17% | 83.90% | 84.74% | 84.57%
14 | 5 | 5 | 16, 32, 64, 128, 256 | (C-P)$\times$5-S | 77.79% | 84.39% | 84.30% | 83.90% | 85.08%
15 | 6 | 3 | 16, 32, 64, 128, 256, 512 | C-C-C-(C-P)$\times$3-S | 70.73% | 86.77% | 86.52% | 87.62% | 87.30%
16 | 6 | 3 | 16, 32, 64, 128, 256, 512 | (C-C-P)$\times$3-S | 73.59% | 87.09% | 86.71% | 87.83% | 87.21%
17 | 6 | 4 | 16, 32, 64, 128, 256, 512 | C-C-(C-P)$\times$4-S | 77.69% | 87.63% | 87.18% | 88.18% | 87.94%
18 | 6 | 4 | 16, 32, 64, 128, 256, 512 | C-P-(C-C-P)$\times$2-C-P-S | 78.38% | 87.65% | 87.70% | 87.80% | 87.36%
19 | 6 | 5 | 16, 32, 64, 128, 256, 512 | C-(C-P)$\times$5-S | 81.89% | 88.60% | 88.08% | 88.45% | 88.97%
20 | 6 | 6 | 16, 32, 64, 128, 256, 512 | (C-P)$\times$6-S | 84.88% | 88.85% | 88.25% | 88.37% | 87.90%

(a) GAA of C1 P1 w.r.t. Polynomial Order

(b) GAA of C2 P2 w.r.t. Polynomial Order

(c) GAA of C3 P3 w.r.t. Polynomial Order

(d) GAA of C4 P4 w.r.t. Polynomial Order

(e) GAA of C5 P5 w.r.t. Polynomial Order

(f) GAA of C6 P6 w.r.t. Polynomial Order

(g) GAA of C6 P3 w.r.t. Polynomial Order

(h) GAA of C6 P5 w.r.t. Polynomial Order

Figure 3: GAA of the same amount of graph convolutional and pooling layers
w.r.t. different polynomial approximation order.

###  II-E Evaluation Metrics

To evaluate performances, multiple metrics were adopted including the global
average accuracy (GAA), Cohen’s Kappa coefficient (Kappa), single class
accuracy on each task, Macro-averaged precision, recall, F1 score, Receiver
Operating Characteristic Curve (ROC curve), and the Area Under ROC Curve
(AUC). Plus, to validate whether the performance differences were significant,
the pair-wise t-test was applied. In this work, the significance level of the
t-test was set as 0.05[52].

##  III Results and Discussion

###  III-A Novel Deep Learning Framework of the GCNs

To explore an optimal model for the GCNs-Net, as detailed in Table II,
multiple structures were investigated to grasp the impacts of the decoding
performance through changing the network hyperparameters, such as the amount
of graph convolutional (Conv) and pooling layers, the polynomial order of
Chebyshev polynomials for filters, and the number of convolutional filters. C
denotes a graph convolutional layer, P denotes a graph pooling layer, F
denotes an FC layer, and S denotes a Softmax layer. The PhysioNet Dataset was
chosen to explore different architectures’ performance, as it contained the
largest number of participants in the field of EEG MI. 20 subjects’ dataset
(S1$\sim$S20) with 1,075,200 samples was used to train and evaluate different
models.

First of all, while holding the number of graph Conv and pooling layers,
experiments were carried out by changing the Chebyshev polynomial order from
$1^{st}$ to $5^{th}$ described in Table II. Figure 3 displays that when there
was only one graph Conv layer followed by a graph pooling layer, the order of
Chebyshev did not make a difference. The GAAs regarding each order were less
than 58%. They overlapped with each other and fluctuated during training.
Additionally, when the number of graph Conv layers was greater than one, the
GAA of the model with the $1^{st}$ order approximation ascended. In the later
epochs, it entered a period of dormancy. In particular, the GAA regarding
models with the $1^{st}$ order witnessed a rugged and abrupt ascent when there
were more graph Conv layers. Apart from the number of graph Conv layers, the
GAAs of models moved upward smoothly with the increasing number of pooling
layers.

As illustrated in Figure 3, the GAA regarding models with the $1^{st}$
polynomial order was unsatisfactory. By contrast, the GAAs of models regarding
the $2^{nd}$ to $5^{th}$ polynomial order were doing a different climb. But
they almost overlapped and paralleled with each other during training. It
indicated that when the order of polynomial approximation was greater than
one, there was a minor impact on the EEG MI decoding. As a result, $2^{nd}$
order Chebyshev approximation for filters was employed to not only capture a
superior result but also reduce the model complexity.

Besides, the impacts of performance by changing the number of graph Conv and
pooling layers were reviewed at a specific polynomial approximation order.

(a) GAA of Different Models
w.r.t. $1^{st}$ Order Chebyshev Polynomial

(b) GAA of Different Models
w.r.t. $2^{nd}$ Order Chebyshev Polynomial

(c) GAA of Different Models
w.r.t. $5^{th}$ Order Chebyshev Polynomial

(d) GAA of Top Ten Models

Figure 4: GAA of different models while applying the same polynomial order.

Figure 4 demonstrates that when the number of graph Conv layers increased, the
GAAs took a steep climb. Notably, as for the $2^{nd}$ polynomial order in
Figure 4, it illustrated that the number of Conv layers did affect the
performance. While applying a deeper model, including extra Conv layers,
features can be better extracted from EEG signals. Meanwhile, the effects of
the number of graph pooling layers were also investigated. The number of
pooling layers promoted and enhanced the decoding performance, but with a
modest increment. As detailed in Figure 4, when the polynomial order was 2nd,
the GAAs were 88.6040% (Model C6-P5-K2), and 88.8523% (Model C6-P6-K2).

TABLE III: Model performance comparison by changing the Num. of FC layers and
filter’s amount

Model | Filter’s Amount |  | Num of
---
FC Layers
| Num of
---
Neurons at
FC Layer
Model Framework | GAA
Base | 16, 32, 64, 128, 256, 512 | 1 | 4 | (C-P)$\times$6-S | 88.85%
1 | 32, 64, 128, 256, 512, 1024 | 1 | 4 | (C-P)$\times$6-S | 90.60%
2 | 64, 128, 256, 512, 1024, 1536 | 1 | 4 | (C-P)$\times$6-S | 90.89%
3 | 16, 32, 64, 128, 256, 512 | 2 | 64, 4 | (C-P)$\times$6-F-S | 88.08%
4 | 16, 32, 64, 128, 256, 512 | 2 | 512, 4 | (C-P)$\times$6-F-S | 88.64%
5 | 16, 32, 64, 128, 256, 512 | 3 | 512, 64, 4 | (C-P)$\times$6-F$\times$2-S | 87.36%
6 | 16, 32, 64, 128, 256, 512 | 3 | 512, 256, 4 | (C-P)$\times$6-F$\times$2-S | 88.35%
7 | 32, 64, 128, 256, 512, 1024 | 3 | 512, 64, 4 | (C-P)$\times$6-F$\times$3-S | 90.45%

Furthermore, based on the optimal C6-P6-K2 model which contained six graph
Conv layers with the $2^{nd}$ polynomial order for filters, and six graph
pooling layers, the influence on performance by changing the number of Conv
filters at every Conv layer and the number of FC layers was also explored in
Table III.

(a) GAA of models with different Num. of
FC layers and convolutional filter’s amount.

(b) Loss of models with different Num. of
FC layers and convolutional filter’s amount.

Figure 5: GAA and loss of models with different Num. of FC layers and
convolutional filter’s amount.

Figure 5 shows that when there were more filters at the graph Conv layers, the
GAAs ascended marginally. However, as indicated in Figure 5, the loss value
rose slightly after a fall. It meant that the model with more filters had
caused overfitting. The reason was that the model structure became much more
complicated while applying more filters. Consequently, 16, 32, 64, 128, 256,
and 512 filters were used for the six-layer GCNs-Net to prevent overfitting.
Meanwhile, there was a gentle descendent when added more FC layers. As a
result, a softmax layer was directly implemented without applying extra FC
layers.

###  III-B Subject-level Validation

The GCNs-Net was validated on 10 subjects from the PhysioNet Dataset, each
with 53,760 samples. The decoding performance (accuracy) was as follows: S1
(97.08%), S2 (90.70%), S3 (97.92%), S4 (96.86%), S5 (80.49%), S6 (89.55%), S7
(84.82%), S8 (97.40%), S9 (97.02%), S10 (98.72%).

(a) PCC Matrix for Subject 10

(b) Absolute PCC Matrix for Subject 10

(c) Adjacency Matrix for Subject 10

(d) Graph Laplacian for Subject 10

(e) PCC Matrix for Subject 5

(f) Absolute PCC Matrix for Subject 5

(g) Adjacency Matrix for Subject 5

(h) Graph Laplacian for Subject 5

Figure 6: The PCC, Absolute PCC, adjacency matrices, and graph Laplacian for
Subject 10 and 5 from the PhysioNet Dataset.

According to Figure 6, the PCC matrix, absolute PCC matrix, adjacency matrix,
and graph Laplacian of Subject 10 and 5 were shown, which achieved the highest
and the lowest GAA. There were quite a lot of variations underlying inter-
subject EEG signals. For each one of the 10 subjects, 98.72% maximum accuracy
was achieved. As for the model of Subject 10, the AUC was 0.991. The single
class accuracy on L, R, B, and F for subject 10 was 99.92%, 97.96%, 98.08%,
and 98.93%, respectively. For 10 subjects, the highest F1 Score was 98.71%,
and the lowest was 80.19%.

Meanwhile, the presented GCNs-Net was validated on the High Gamma Dataset. The
model containing 14 subjects was separately trained and evaluated. Since there
were 44 electrodes[22], the maximum number of pooling layers was 2. We used
(C-C-C-P)$\times$2 architecture of the GCNs-Net to decode EEG tasks. From
Subject 1 to Subject 14, 96.43%, 95.63%, 93.04%, 99.18%, 98.65%, 94.77%,
93.49%, 97.91%, 95.48%, 96.77%, 98.55%, 98.69%, 98.34%, and 90.43% accuracies
were achieved, separately. The mean accuracy was 96.24%. The results indicated
that the GCNs-Net manages to handle individual variability due to its
robustness and effectiveness.

###  III-C Classification at the Group Level

Next, the GCNs-Net was evaluated at a group of 20 subjects (S1$\sim$S20) from
the PhysioNet Dataset. The GAA, Kappa, precision, recall, and F1 Score were
88.35%, 84.47%, 88.39%, 88.35%, and 88.34%, respectively. Further, the single
class accuracies on L, R, B, and F were 83.45%, 86.72%, 83.96%, and 99.42%.
The method performed well on the classes F and B. The AUC was 0.922.

Besides, it was evaluated on the High Gamma Dataset. 14 subjects’ data was
used in the experiment, with 12,504,375 samples in total. 80.89% GAA and
80.78% F1 Score were achieved, separately. The reason for the outcome was that
the GCNs-Net could converge for group-wise predictions, and succeed in
extracting relevant features from EEG signals.

###  III-D 10-fold Cross-validation for Reliability

The GCNs-Net was trained at the group level (20 subjects (S1$\sim$S20) from
the PhysioNet Dataset) based on the 10-fold cross-validation to validate the
stability and reliability. In each experiment, the dataset has been randomly
shuffled, following which it was divided into ten pieces, where the first nine
parts were used as the training set, and the last one served as the test set.
The procedure has been carried out ten times, thus provided us ten results.

(a) GAA of repetitive experiments
for 20 subjects comparison.

(b) F1 Score of repetitive experiments
for 20 subjects comparison.

Figure 7: GAA and F1 Score of repetitive experiments for 20 subjects from the
PhysioNet Dataset.

With the result from Section III-A (88.852% GAA, Model C6-P6-K2), 11 results
were listed in Figure 7. 89.387% maximum GAA was achieved, and the lowest was
87.896%. The averaged GAAs and F1 Score were both 88.57%. At the group level,
the performance was stably reproducible through the repetitive experiments for
cross-validation, showing reliability and stability of the GCNs-Net.

###  III-E Robustness to Data Size

It has also been trained and evaluated on different amounts of participants.
50 subjects’ dataset (S1$\sim$S50) with 2,688,000 samples and 100 subjects’
dataset (S1$\sim$S100) with 5,376,000 samples from the PhysioNet Dataset were
used. GAAs and loss were illustrated in Figure 8.

(a) GAA w.r.t. various sizes of dataset.

(b) Loss w.r.t. various sizes of dataset.

Figure 8: GAA and loss w.r.t. various sizes of dataset (20, 50, 100 subjects)
from the PhysioNet Dataset.

The GCNs-Net has accomplished 89.75% (test set) and 94.99% (training set)
accuracies for a group of 50 subjects. Further, it can also achieve 88.14%
(test set) and 93.24% (training set) accuracies for 100 subjects. The results
showed that it could learn generalized features from subjects at a larger
scale, which had a lower signal-to-noise ratio. The reason was that the GCNs-
Net could handle a larger amount of subjects, suggested adaptability and
robustness to individual variability. The presented method was much more
effective and robust to copy with graph-structured EEG-based MI signals since
it has considered the functional topological relationship of EEG electrodes.

###  III-F Comparison with State-of-the-art

Two levels were applied to compare the performance of the competitive models,
i.e., either from the subject or the group level. The decoding performance
mainly measured by the maximum GAA[18], averaged (mean) GAA[34] on two
datasets. First of all, multiple approaches were performed on the PhysioNet
Dataset[20, 53, 18]. The performance was compared in Table IV.

TABLE IV: Performance comparison on the PhysioNet Dataset

Related Work | Max. GAA | Avg. GAA |  _p_ -value | Level | Approach | Num of Subjects
---|---|---|---|---|---|---
Dose _et al._ (2018)[20] | - | 58.58% | $-$ | Group | CNNs | 105
80.38% | 68.51% | $<0.05$ | Subject | 1
Ma _et al._ (2018)[53] | 82.65% | 68.20% | $-$ | Group | RNNs | 12
Hou _et al._ (2020)[18] | 94.50% | $-$ | $-$ | Group | ESI-CNNs | 10
96.00% | $-$ | $>0.05$ | Subject | 1
Author | 89.387% | 88.57% | $-$ | Group | GCNs-Net |  | 20
---
88.14% | - | 100
98.72% | 93.056% | Subject | 1

The ESI$-$CNNs approach has achieved 94.50% maximum accuracy at the group
level (10 subjects)[18]. By contrast, the GCNs-Net has obtained competing
performance, 89.387% maximum accuracy for a group of 20 participants, and
88.14% for 100 participants. Meanwhile, at the subject level, the $p$-value
between the GCNs-Net and the CNNs model (Dose _et al._ (2018)[20]) was
significantly less than 0.05 ($p<0.05$). It indicated a significant difference
in predictive performance between two models, and the GCNs-Net was superior to
predict EEG tasks, with a 30.21% maximum accuracy increment. However, compared
with the ESI$-$CNNs approach[18], there was no significant difference in
classification performance at a 95% confidence interval as the $p$-value was
greater than 0.05 ($p>0.05$). Furthermore, the GCNs-Net attained the best
state-of-the-art performance at the hundred-subject level on the PhysioNet
Dataset, which far exceeded current studies. The reason for the outcome was
that the presented approach maintained robust and effective on the dataset
with a larger amount of participants, regardless of the inter-trial and inter-
subject variations.

We compared the classification performances of the most recently published
literature on the High Gamma Dataset[22, 54, 27].

TABLE V: Performance comparison on the High Gamma Dataset

Related Work | Avg. GAA |  _p_ -value | Level | Approach | Dataset
---|---|---|---|---|---
Schirrmeister _et al._ (2017)[22] | 92.50% | $<0.05$ | Subject | CNNs |  1 subjects
Li _et al._ (2019)[54] | 93.70% | $<0.05$ | CP-MixedNet
Tang _et al._ (2019)[27] | 95.30% | $>0.05$ | DAN
Author | 80.89% | $-$ | Group | GCNs-Net |  | 14 subject
---
96.24% | Subject |  | 1 subject
---

Evaluated on the High Gamma Dataset, the $p$-values were both less than 0.05
($p<0.05$) when comparing the introduced method with the CNN-based methods
(Schirrmeister _et al._ (2017)[22]) and the CP-MixedNet (Li _et al._
(2019)[54]). The performance was statistically significantly different among
the models. The GCNs-Net successfully predicts MI tasks with dominant
performances (99.18% maximum accuracy, and 96.24% averaged accuracy). The
$p$-value compared with the DAN approach (Tang _et al._ (2019)[27]) was
greater than 0.05 ($p>0.05$). The performance of the two models was less
different statistically, and both models already achieved very high
performance. Last but not least, dominant classification accuracy has verified
the robustness and effectiveness of our presented GCNs-Net.

##  IV Conclusion

In general, the GCNs-Net, a novel deep learning framework based on the GCNs,
was presented to distinguish four-class MI intentions by cooperating with the
functional topological relationship of EEG electrodes. The introduced method
has been proven to converge for both personalized and group-wise predictions.
Trained with individual data, the approach has achieved an averaged accuracy
of 93.056% (PhysioNet Dataset) and 96.24% (High Gamma Dataset) on predicting
the independent trials of the same participant, which is dominant in existing
studies, indicating that the GCNs-Net converges well for individuals.
Moreover, it has reached the uppermost accuracy on numerous sizes of group-
level prediction on the PhysioNet Dataset, i.e., with 89.387% accuracy for 20
subjects, 89.75% for 50 subjects, and 88.14% for 100 subjects, which implies
that it is considerably robust to individual variability. Further, it held an
averaged accuracy of 88.57% after 10-fold cross-validation showing reliability
and stability. On the other hand, it predicted all of the four MI tasks with
superior accuracy, the best among which was the two feet prediction with an
accuracy of 99.42%. It indicated that the introduced method was able to build
a generalized representation against both personalized and group-wise
variations. It can be performed to decode any EEG MI signals, and other EEG-
based graph-structured data towards developing effective and efficient BCI
systems.

## References

  * [1] G. Santhanam, S. I. Ryu, M. Y. Byron, A. Afshar, and K. V. Shenoy, “A high-performance brain–computer interface,” Nature, vol. 442, no. 7099, p. 195, 2006.
  * [2] M. A. Silver, D. Ress, and D. J. Heeger, “Topographic maps of visual spatial attention in human parietal cortex,” Journal of neurophysiology, vol. 94, no. 2, pp. 1358–1371, 2005.
  * [3] D. L. Schomer and F. L. Da Silva, Niedermeyer’s electroencephalography: basic principles, clinical applications, and related fields.  Lippincott Williams & Wilkins, 2012.
  * [4] J. Hubbard, A. Kikumoto, and U. Mayr, “Eeg decoding reveals the strength and temporal dynamics of goal-relevant representations,” Scientific reports, vol. 9, no. 1, p. 9051, 2019.
  * [5] D. J. McFarland, L. A. Miner, T. M. Vaughan, and J. R. Wolpaw, “Mu and beta rhythm topographies during motor imagery and actual movements,” Brain topography, vol. 12, no. 3, pp. 177–186, 2000.
  * [6] H. K. Lee and Y.-S. Choi, “A convolution neural networks scheme for classification of motor imagery eeg based on wavelet time-frequecy image,” in 2018 International Conference on Information Networking (ICOIN), pp. 906–909, IEEE, 2018.
  * [7] D. S. Bassett and O. Sporns, “Network neuroscience,” Nature neuroscience, vol. 20, no. 3, p. 353, 2017.
  * [8] J. Lv, X. Jiang, X. Li, D. Zhu, S. Zhang, S. Zhao, H. Chen, T. Zhang, X. Hu, J. Han, et al., “Holistic atlases of functional networks and interactions reveal reciprocal organizational architecture of cortical function,” IEEE Transactions on Biomedical Engineering, vol. 62, no. 4, pp. 1120–1131, 2014.
  * [9] J. Lv, V. T. Nguyen, J. van der Meer, M. Breakspear, and C. C. Guo, “N-way decomposition: Towards linking concurrent eeg and fmri analysis during natural stimulus,” in International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 382–389, Springer, 2017.
  * [10] S. Taran, V. Bajaj, D. Sharma, S. Siuly, and A. Sengur, “Features based on analytic imf for classifying motor imagery eeg signals in bci applications,” Measurement, vol. 116, pp. 68–76, 2018.
  * [11] S. Chaudhary, S. Taran, V. Bajaj, and S. Siuly, “A flexible analytic wavelet transform based approach for motor-imagery tasks classification in bci applications,” Computer Methods and Programs in Biomedicine, p. 105325, 2020.
  * [12] J. Kevric and A. Subasi, “Comparison of signal decomposition methods in classification of eeg signals for motor-imagery bci system,” Biomedical Signal Processing and Control, vol. 31, pp. 398–406, 2017.
  * [13] M. T. Sadiq, X. Yu, Z. Yuan, Z. Fan, A. U. Rehman, G. Li, and et al., “Motor imagery eeg signals classification based on mode amplitude and frequency components using empirical wavelet transform,” IEEE Access, vol. 7, pp. 127678–127692, 2019.
  * [14] S. Taran and V. Bajaj, “Motor imagery tasks-based eeg signals classification using tunable-q wavelet transform,” Neural Computing and Applications, pp. 1–8, 2018.
  * [15] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” Nature, vol. 521, no. 7553, p. 436, 2015.
  * [16] V. J. Lawhern, A. J. Solon, N. R. Waytowich, S. M. Gordon, C. P. Hung, and B. J. Lance, “Eegnet: a compact convolutional neural network for eeg-based brain–computer interfaces,” Journal of neural engineering, vol. 15, no. 5, p. 056013, 2018.
  * [17] O. Faust, Y. Hagiwara, T. J. Hong, O. S. Lih, and U. R. Acharya, “Deep learning for healthcare applications based on physiological signals: A review,” Computer methods and programs in biomedicine, vol. 161, pp. 1–13, 2018.
  * [18] Y. Hou, L. Zhou, S. Jia, and X. Lun, “A novel approach of decoding eeg four-class motor imagery tasks via scout esi and cnn,” Journal of neural engineering, vol. 17, no. 1, p. 016048, 2020.
  * [19] S. U. Amin, M. Alsulaiman, G. Muhammad, M. A. Mekhtiche, and M. S. Hossain, “Deep learning for eeg motor imagery classification based on multi-layer cnns feature fusion,” Future Generation Computer Systems, vol. 101, pp. 542–554, 2019.
  * [20] H. Dose, J. S. Møller, H. K. Iversen, and S. Puthusserypady, “An end-to-end deep learning approach to mi-eeg signal classification for bcis,” Expert Systems with Applications, vol. 114, pp. 532–542, 2018.
  * [21] S. Chaudhary, S. Taran, V. Bajaj, and A. Sengur, “Convolutional neural network based approach towards motor imagery tasks eeg signals classification,” IEEE Sensors Journal, vol. 19, no. 12, pp. 4494–4500, 2019.
  * [22] R. T. Schirrmeister, J. T. Springenberg, L. D. J. Fiederer, M. Glasstetter, K. Eggensperger, M. Tangermann, and et al., “Deep learning with convolutional neural networks for eeg decoding and visualization,” Human brain mapping, vol. 38, no. 11, pp. 5391–5420, 2017.
  * [23] C. J. Ortiz-Echeverri, S. Salazar-Colores, J. Rodríguez-Reséndiz, and R. A. Gómez-Loenzo, “A new approach for motor imagery classification based on sorted blind source separation, continuous wavelet transform, and convolutional neural network,” Sensors, vol. 19, no. 20, p. 4541, 2019\.
  * [24] D. Li, J. Wang, J. Xu, and X. Fang, “Densely feature fusion based on convolutional neural networks for motor imagery eeg classification,” IEEE Access, vol. 7, pp. 132720–132730, 2019.
  * [25] R. Zhang, Q. Zong, L. Dou, and X. Zhao, “A novel hybrid deep learning scheme for four-class motor imagery classification,” Journal of neural engineering, vol. 16, no. 6, p. 066004, 2019.
  * [26] R. Alazrai, M. Abuhijleh, H. Alwanni, and M. I. Daoud, “A deep learning framework for decoding motor imagery tasks of the same hand using eeg signals,” IEEE Access, vol. 7, pp. 109612–109627, 2019.
  * [27] X. Tang and X. Zhang, “Conditional adversarial domain adaptation neural network for motor imagery eeg decoding,” Entropy, vol. 22, no. 1, p. 96, 2020.
  * [28] A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, and et al., “Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals,” Circulation, vol. 101, no. 23, pp. e215–e220, 2000.
  * [29] D. Zhang, L. Yao, K. Chen, S. Wang, X. Chang, and Y. Liu, “Making sense of spatio-temporal preserving representations for eeg-based human intention recognition,” IEEE transactions on cybernetics, 2019.
  * [30] C. Brunner, R. Leeb, G. Müller-Putz, A. Schlögl, and G. Pfurtscheller, “Bci competition 2008–graz data set a,” Institute for Knowledge Discovery (Laboratory of Brain-Computer Interfaces), Graz University of Technology, vol. 16, 2008.
  * [31] M. Tangermann, K.-R. Müller, A. Aertsen, N. Birbaumer, C. Braun, C. Brunner, and et al., “Review of the bci competition iv,” Frontiers in neuroscience, vol. 6, p. 55, 2012.
  * [32] D. Zhang, L. Yao, K. Chen, and J. Monaghan, “A convolutional recurrent attention model for subject-independent eeg signal analysis,” IEEE Signal Processing Letters, vol. 26, no. 5, pp. 715–719, 2019.
  * [33] D. Zhang, K. Chen, D. Jian, and L. Yao, “Motor imagery classification via temporal attention cues of graph embedded eeg signals,” IEEE Journal of Biomedical and Health Informatics, 2020.
  * [34] D. Zhang, L. Yao, K. Chen, S. Wang, P. D. Haghighi, and C. Sullivan, “A graph-based hierarchical attention model for movement intention detection from eeg signals,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 27, no. 11, pp. 2247–2253, 2019.
  * [35] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in Advances in neural information processing systems, pp. 3844–3852, 2016.
  * [36] T. Zhang, X. Wang, X. Xu, and C. P. Chen, “Gcb-net: Graph convolutional broad network and its application in emotion recognition,” IEEE Transactions on Affective Computing, 2019.
  * [37] Z. Wang, Y. Tong, and X. Heng, “Phase-locking value based graph convolutional neural networks for emotion recognition,” IEEE Access, vol. 7, pp. 93711–93722, 2019.
  * [38] T. Song, W. Zheng, S. Peng, and C. Zhen, “Eeg emotion recognition using dynamical graph convolutional neural networks,” IEEE Transactions on Affective Computing, vol. PP, no. 99, pp. 1–1, 2018.
  * [39] X.-h. Wang, T. Zhang, X.-m. Xu, L. Chen, X.-f. Xing, and C. P. Chen, “Eeg emotion recognition using dynamical graph convolutional neural networks and broad learning system,” in 2018 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), pp. 1240–1244, IEEE, 2018.
  * [40] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” in Advances in Neural Information Processing Systems, pp. 1024–1034, 2017.
  * [41] F. Monti, D. Boscaini, J. Masci, E. Rodola, J. Svoboda, and M. M. Bronstein, “Geometric deep learning on graphs and manifolds using mixture model cnns,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 5115–5124, 2017.
  * [42] M. Niepert, M. Ahmed, and K. Kutzkov, “Learning convolutional neural networks for graphs,” in International conference on machine learning, pp. 2014–2023, 2016.
  * [43] H. Gao, Z. Wang, and S. Ji, “Large-scale learnable graph convolutional networks,” in Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 1416–1424, ACM, 2018.
  * [44] J. Bruna, W. Zaremba, A. Szlam, and Y. LeCun, “Spectral networks and locally connected networks on graphs,” arXiv preprint arXiv:1312.6203, 2013.
  * [45] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” IEEE signal processing magazine, vol. 30, no. 3, pp. 83–98, 2013.
  * [46] Z. Jiao, X. Gao, Y. Wang, J. Li, and H. Xu, “Deep convolutional neural networks for mental load classification based on eeg data,” Pattern Recognition, vol. 76, pp. 582–595, 2018.
  * [47] R. Chai, S. H. Ling, P. P. San, G. R. Naik, T. N. Nguyen, Y. Tran, A. Craig, and H. T. Nguyen, “Improving eeg-based driver fatigue classification using sparse-deep belief networks,” Frontiers in neuroscience, vol. 11, p. 103, 2017.
  * [48] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and P. S. Yu, “A comprehensive survey on graph neural networks,” arXiv preprint arXiv:1901.00596, 2019\.
  * [49] D. K. Hammond, P. Vandergheynst, and R. Gribonval, “Wavelets on graphs via spectral graph theory,” Applied and Computational Harmonic Analysis, vol. 30, no. 2, pp. 129–150, 2011.
  * [50] I. S. Dhillon, Y. Guan, and B. Kulis, “Weighted graph cuts without eigenvectors a multilevel approach,” IEEE transactions on pattern analysis and machine intelligence, vol. 29, no. 11, pp. 1944–1957, 2007.
  * [51] J. Shi and J. Malik, “Normalized cuts and image segmentation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 22, no. 8, pp. 888–905, 2000.
  * [52] P. Zhang, X. Wang, W. Zhang, and J. Chen, “Learning spatial–spectral–temporal eeg features with recurrent 3d convolutional neural networks for cross-task mental workload assessment,” IEEE Transactions on neural systems and rehabilitation engineering, vol. 27, no. 1, pp. 31–42, 2018.
  * [53] X. Ma, S. Qiu, C. Du, J. Xing, and H. He, “Improving eeg-based motor imagery classification via spatial and temporal recurrent neural networks,” in 2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 1903–1906, IEEE, 2018.
  * [54] Y. Li, X.-R. Zhang, B. Zhang, M.-Y. Lei, W.-G. Cui, and Y.-Z. Guo, “A channel-projection mixed-scale convolutional neural network for motor imagery eeg decoding,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 27, no. 6, pp. 1170–1180, 2019.


# Boosted dark matter from diffuse supernova neutrinos

Anirban Das<EMAIL_ADDRESS>SLAC National Accelerator Laboratory,
Stanford University, 2575 Sand Hill Road, Menlo Park, CA 94025, USA  Manibrata
Sen<EMAIL_ADDRESS>Department of Physics, University of California
Berkeley, Berkeley, California 94720, USA  Northwestern University, Department
of Physics & Astronomy, 2145 Sheridan Road, Evanston, IL 60208, USA

(June 11, 2021)

###### Abstract

The XENON collaboration recently reported an excess of electron recoil events
in the low energy region with a significance of around $3.3\sigma$. An
explanation of this excess in terms of thermal dark matter seems challenging.
We propose a scenario where dark matter in the Milky Way halo gets boosted as
a result of scattering with the diffuse supernova neutrino background. This
interaction can accelerate the dark-matter to semi-relativistic velocities,
and this flux, in turn, can scatter with the electrons in the detector,
thereby providing a much better fit to the data. We identify regions in the
parameter space of dark-matter mass and interaction cross-section which
satisfy the excess. Furthermore, considering the data only hypothesis, we also
impose bounds on the dark-matter scattering cross-section, which are
competitive with bounds from other experiments.

††preprint: SLAC-PUB-17591, N3AS-21-005

##  I Introduction

Dark Matter (DM) constitutes a whopping one-fourth of the energy budget of the
Universe, yet we remain in the dark about its origin and nature. The vanilla
DM candidate is expected to have little or no interactions with ordinary
matter beyond gravity, and hence poses a challenge to detect. Battling the
odds, the scientific community has made vast progress in the past few decades,
both in terms of theory, as well as detection techniques [1]. While no
conclusive evidence has emerged, in their recent Science Run 1 (SR1) data, the
XENON1T experiment has seen an excess of electron-recoil events in the low-
energy data $(2-3\,{\rm keV})$ [2]. More specifically, in the energy range
$1-7{\rm\ keV}$, 285 events have been observed compared to an $232\pm 15$
events expected from the background-only model. This $3.3\sigma$ excess in the
electron-recoil events has stirred up quite an excitement.

The XENON collaboration has offered three possible beyond-the-Standard-Model
(BSM) explanations of this excess– (i) solar axions, (ii) anomalous neutrino
magnetic moments, and (iii) bosonic dark matter [2]. However, the first two
solutions run into trouble with astrophysical cooling bounds [3], while the
third solution produces a mono-energetic peak instead of an excess in a few
bins. Other possible origin of the excess is attributed to unresolved
$\beta$-decays of tritium [2]. A host of other papers followed suit with
various intriguing BSM explanations [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 6, 18, 19, 20, 21, 22, 23, 24]. All these works emphasized that an
explanation in terms of the usual thermal DM seems difficult, owing to its
non-relativistic nature.

So, is a DM interpretation of the signal completely ruled out? Recent works
have demonstrated that the XENON1T excess can be vindicated by a component of
dark matter that gets boosted towards the earth [25, 26, 27]. The idea of
boosted dark matter (BDM) is not new; a number of mechanisms have been put
forward to render a non-thermal component to the DM velocity by scattering
with cosmic rays, or annihilation of heavier dark candidates [28, 29, 30, 31,
32, 33, 34]. Such a BDM can lead to more energetic signatures in current as
well as upcoming direct detection and neutrino experiments [32, 35, 36, 37].

DM boosted by neutrinos is, on the other hand, a relatively new concept.
Neutrinos interact very weakly with other SM particles, and are most elusive
among all particles. This naturally makes one wonder if neutrinos have any
secret interactions with DM. If true, then such interaction could alter the
energy spectrum of the DM particles. Implications of such processes have been
studied in a few cosmological and astrophysical contexts [38, 39, 40, 24]. It
is also interesting to study scenarios where DM interacts only with the
leptons in SM [41, 42, 43].

Our universe is bathed in a sea of MeV-neutrinos, emerging from massive stars
going supernova, right from the epoch of first star formation. The diffuse
supernova neutrino background (DSNB) is an isotropic flux of neutrinos and
antineutrinos of all flavors, emitted from all core-collapse supernovae
(CCSNe) in our Universe [44, 45]. The DSNB provides a perfect astrophysical
laboratory to test fundamental neutrino physics [46].

In this paper, we consider the novel idea that the DM in the Milky Way (MW)
halo experiences scattering with the DSNB, and gets boosted to semi-
relativistic velocities. Without alluding to any specific models, we envision
a scenario where DM couples to neutrinos and electrons with equal strength.
Such supernova-boosted-dark-matter (SnBDM ) can leave interesting signatures
in low-energy recoil experiments, and even account for the excess observed in
XENON1T. Furthermore, considering the XENON1T low energy data, we impose
constraints on the DM-electron cross-section. Implications of such
interactions of lepotophilic DM for underground detectors are also discussed.

##  II DSNB spectra

Estimating the DSNB flux at Earth requires a precise knowledge of the rate of
CCSN, as well as the flavor-dependent neutrino spectra emitted from a SN. The
CCSN rate (denoted by $R_{\rm CCSN}(z)$) is proportional to the rate of star
formation (SFR), and has been well estimated by astronomical surveys [47]. We
follow [46] (and references therein) to estimate the CCSN rate. The neutrino
emission from a SN follows an approximate thermal distribution, and can be
parameterized by non-degenerate Fermi-Dirac distribution [44],

$F_{\nu}(E)=\frac{E_{\nu}^{\rm
tot}}{6}\,\frac{120}{7\pi^{4}}\,\frac{E_{\nu}^{2}}{T_{\nu}^{4}}\,\frac{1}{e^{E_{\nu}/T_{\nu}}+1}\,,$
(1)

where $T_{\nu}$ is the flavor dependent neutrino temperature, and
$E_{\nu}^{\rm tot}$ is the total energy emitted in the form of neutrinos. As
neutrino emission is dominated by the cooling phase, lasting for at least
$10\,$s post-bounce, it is reasonable to expect that the total emitted energy
is equipartitioned among all the flavors.

The DSNB flux can be calculated as [48, 44]

$\Phi_{\nu}(E)=\int_{0}^{z_{\rm max}}\frac{dz}{H(z)}R_{\rm
CCSN}(z)F_{\nu}(E)\,,$ (2)

where $H(z)=H_{0}\,\sqrt{\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}}$ is the Hubble
function with $H_{0}=67.36\,{\rm km\,s}^{-1}\,{\rm Mpc}^{-1}$, and
$\Omega_{m}$ and $\Omega_{\Lambda}$ give the matter and vacuum contribution to
the energy density [49]. The energy observed at the earth is redshifted such
that $E=E^{\prime}(1+z)$, while the maximum redshift of star-formation is
taken to be $z_{\rm max}\sim 5$.

Figure 1: DSNB flux for each neutrino species. The temperatures are taken to
be $T_{\nu_{e}}=6{\rm\ MeV}$,$T_{\overline{\nu}_{e}}=7{\rm\ MeV}$ and
$T_{\nu_{x}}=10{\rm\ MeV}$, consistent with the bounds imposed by Super-
Kamiokande [50]. The shaded regions represent the uncertainties in the DSNB
spectra. The range of values for $R_{\rm CCSN}(z)$ is taken from [48].

We show the DSNB flux for all three flavors in Fig. 1, with the uncertainties
stemming from the SFR. Due to charged current interactions,
$\nu_{e},\overline{\nu}_{e}$ have a lower temperature than $\nu_{x=\mu,\tau}$,
leading to a larger sensitivity for the electron flavor neutrinos. Note that
we only assume adiabatic Mikheyev-Smirnov-Wolfenstein resonant oscillations in
effect inside the SN [51, 52], and neglect the more uncertain collective
oscillations [53, 54] and shock-related effects [55]. In the energy range we
are interested in, the neutrino fluxes are fairly similar, and hence outcomes
due to collective oscillations would not change our results by more than
$\mathcal{O}(10\%)$ [45]. Given the uncertainty in the SFR is of
$\mathcal{O}(40\%)$ [48], neglecting collective oscillations is a fair
assumption. We operate in the normal mass ordering, where all the $\nu_{e}$
exit the CCSN as a $\nu_{3}$, while the other flavors exit as incoherent
superpositions of the remaining mass states.

##  III Boosted Dark Matter

Without referring to any particular model, we assume that the DM particle
$\chi$, of mass $m_{\chi}$, couples to electrons and neutrinos with a cross-
section $\sigma_{\chi e}$. This type of scenario can arise in several
leptophilic neutrino-portal models of particle DM [42, 56, 41, 43, 57, 58, 59,
60]. For concreteness, we assume that the mediator is heavier than typical
energy transfer of a few ${\rm\ keV}$ to the electrons, but lighter than
$m_{\chi}$. This allows for large values of the DM-electron cross-section for
reasonably small couplings.

The mechanism of boosting DM particles with neutrinos is similar to that of
cosmic-ray boosted DM, which has been studied extensively in the literature
[32, 33]. We follow their work to compute the SnBDM flux. The DSNB neutrinos
have energies of order $E_{\nu}\sim\mathcal{O}(10{\rm\ MeV})$. The DM
particles in the halo can be assumed to be at rest relative to the neutrinos.
This implies that the initial virial velocity of DM is not relevant since
collisions with neutrinos will impart a much higher velocity to it. We shall
also neglect the neutrino mass throughout this work. With these assumptions,
the energy $T_{\chi}$ transferred to $\chi$ by a neutrino in a single
scattering is given by

$\displaystyle
T_{\chi}=T_{\chi}^{\mathrm{max}}\left(\frac{1-\cos\theta_{\mathrm{CM}}}{2}\right)\,,$
(3) $\displaystyle
T_{\chi}^{\mathrm{max}}=\frac{E_{\nu}^{2}}{E_{\nu}+m_{\chi}/2}\,.$

Here $\theta_{\mathrm{CM}}$ is the scattering angle in the center-of-mass
frame, and $T_{\chi}^{\mathrm{max}}$ is the maximum energy that can be
transferred by a neutrino of energy $E_{\nu}$.

Figure 2: The SnBDM flux and velocity distribution for $m_{\chi}=300{\rm\
MeV}$ and $\sigma_{\chi e}=4\times 10^{-30}{\rm\,cm^{2}}$. The rescaled
Standard Halo Model (SHM) velocity distribution is shown in dashed gray line.

The DM density $\rho_{\chi}(r)$ in the MW halo is assumed to be given by a
Navarro-Frenk-White (NFW) profile where $r$ is the radial distance from the
Galactic center (GC) [61]. Looking from the earth, the DM density is a
function of both line-of-sight (los) distance $l$ and its direction $\Omega$.
Therefore, the total upscattered flux of the DM integrated over the whole MW
halo is given by,

$\frac{d\phi_{\chi}}{dE_{\nu}}=\int\frac{d\Omega}{4\pi}\int
dl\frac{\rho_{\chi}}{m_{\chi}}\Phi_{\nu}(E_{\nu})\sigma_{\chi
e}=D_{\mathrm{halo}}\Phi_{\nu}(E_{\nu})\frac{\sigma_{\chi e}}{m_{\chi}}\,.$
(4)

Here, $\Phi_{\nu}(E_{\nu})$ is the DSNB flux computed in the previous section.
In the last step, we have factored the whole quantity into three mutually
independent parts: $D_{\mathrm{halo}}$ encodes the los and the angular
integrals over the DM density $\rho_{\chi}(r)$ and depends only on the halo
profile, $\Phi_{\nu}(E_{\nu})$ is the DSNB flux, and $\sigma_{\chi
e},\,m_{\chi}$ represent the new BSM parameters of the theory. We find
$D_{\mathrm{halo}}=2.02\times 10^{25}{\rm\ MeV}\,\mathrm{cm^{-2}}$. Using
Eq.(4), the boosted DM flux can be written as

$\frac{d\phi_{\chi}}{dT_{\chi}}=\int
dE_{\nu}\frac{d\phi_{\chi}}{dE_{\nu}}\frac{1}{T_{\chi}^{\mathrm{max}}(E_{\nu})}\Theta\left[T_{\chi}^{\mathrm{max}}(E_{\nu})-T_{\chi}\right]\,.$
(5)

Here, the Heaviside step function
$\Theta\left[T_{\chi}^{\mathrm{max}}(E_{\nu})-T_{\chi}\right]$ ensures that
for a given neutrino energy, there is a maximum energy upto which the DM can
be boosted. For a heavy-mediator interaction, the differential scattering
cross-section yields a flat distribution. Fig. 2 shows an example flux of the
SnBDM for $m_{\chi}=300\,{\rm\ MeV}$ and $\sigma_{\chi e}=4\times
10^{-30}{\rm\,cm^{2}}$. It is quite evident that the DM flux gets boosted by
orders of magnitude, compared to that in the Standard Halo Model (SHM).

The boosted DM, after traveling through the MW, scatters with the electrons in
the detector at the earth. The formalism for DM-$e$ scattering is similar to
scattering with neutrinos with the following replacements:
$\nu\rightarrow\chi$, and $\chi\rightarrow e$. As mentioned earlier, we assume
the scattering cross-section to be the same as in case of DM-$\nu$ scattering.
Phenomenological implications of this will be discussed later.

Following Eq.(5), we can write the electron recoil spectrum

$\displaystyle\dfrac{d\Gamma}{dT_{e}}=\dfrac{Z_{\mathrm{Xe}}}{m_{\mathrm{Xe}}}\int
dT_{\chi}\dfrac{d\phi_{\chi}}{dT_{\chi}}\dfrac{1}{T^{\mathrm{max}}_{e}(T_{\chi})}\sigma_{\chi
e}\,.$ (6)

Here, $Z_{\mathrm{Xe}}=40$ is the _effective atomic number_ of xenon and takes
care of the fact that all electrons in a xenon atom are not available for
ionization, $m_{\mathrm{Xe}}$ is the mass of a xenon atom, and
$d\phi_{\chi}/dT_{\chi}$ is given by Eq.(5) [2]. We further convolute this
electron recoil spectrum with a Gaussian detector response function with a
width

$\sigma(E)=a\sqrt{E}+bE\,,$ (7)

with $a=0.31\sqrt{\text{keV}},b=0.0037$ [2]. The XENON1T SR1 data and the
background model spectrum $B_{0}$ are taken from Ref. [2].

Using these details, the corresponding recoil spectrum is shown in Fig. 3 for
a sample value of $m_{\chi}=300\,{\rm\ MeV}$ and $\sigma_{\chi e}=4\times
10^{-30}{\rm\,cm^{2}}$. Clearly, the recoil spectra predicted by the SnBDM
model+$B_{0}$ (red line) presents a much better fit to the XENON1T data than
the $B_{0}$-only hypothesis. In the following section, we perform a detailed
scan of the parameter space of SnBDM to identify regions which can explain the
excess.

Figure 3: The electron recoil spectrum (red, solid) in XENON1T for
$m_{\chi}=300{\rm\ MeV}$ and $\sigma_{\chi e}=4\times 10^{-30}{\rm\,cm^{2}}$.
The black dashed curve represents the contribution from SnBDM only. The black
dotted curve represents the background model $B_{0}$ [2].

##  IV Methodology & Results

We utilize the SR1 data in two ways. Given the $\sim 3\sigma$ excess in
electron recoil events in the $2-5{\rm\ keV}$ region, we try to explain the
excess in terms of the SnBDM, and identify regions in the
$m_{\chi}-\sigma_{\chi e}$ parameter space that yields a better fit to the
data than the $B_{0}$. Furthermore, using only the XENON data, we also do a
general exclusion analysis. These are explained in detail below***The codes
used for the calculations in this paper and the jupyter notebooks used to
generate the plots can be found at this https url..

We define the following test statistic

$\chi^{2}=\sum_{i}\dfrac{\left(D_{i}-M_{i}(m_{\chi},\sigma_{\chi
e})\right)^{2}}{\sigma_{i}^{2}}\,,$ (8)

where, $D_{i}$ is the XENON data, $M_{i}$ denotes the total number of events
predicted by the the SnBDM model and the background, and $\sigma_{i}^{2}\equiv
M_{i}+\sigma_{Di}^{2}$ is the combined uncertainty in model and data, with $i$
running over all energy bins in the electron recoil spectrum. For simplicity,
we neglect the uncertainty arising from SFR in this as well as subsequent
analyses. However, we have explicitly checked that the $40\%$ uncertainty in
DSNB changes the best-fit $\chi^{2}$ by about $0.4$. The resulting $\chi^{2}$
plots corresponding to $1\sigma$ and $2\sigma$ confidence limits are shown in
blue contours in Fig. 4. The best-fit point with $m_{\chi}=286{\rm\ MeV}$ and
$\sigma_{\chi e}=3\times 10^{-30}{\rm\,cm^{2}}$, marked with a white asterisk,
corresponds to a $\Delta\chi^{2}=-5.1$ relative to the background model
$B_{0}$. The contours show a positive correlation between $m_{\chi}$ and
$\sigma_{\chi e}$. This can be understood as follows: heavier DM mass
decreases the probability of electron recoil, so larger cross section is
needed to reproduce the observed spectrum.

We also compute the values of the cross section $\sigma_{\chi e}$ that can be
excluded given the observed electron recoil. To this end, we take the
following conservative approach to construct a test statistic $\chi_{>}^{2}$
[62],

$\chi_{>}^{2}=\sum_{i}\dfrac{\left(D_{i}-M_{i}(\sigma_{\chi
e})\right)^{2}}{\sigma_{i}^{2}}\,.$ (9)

Here, we sum over only those energy bins where the model $M_{i}$ alone
predicts more events than what is observed. This is a conservative method and
stays completely agnostic about any background model. For each value of
$m_{\chi}$, we scan over $\sigma_{\chi e}$ to compute $\chi_{>}^{2}$. This
allows us to compute an exclusion contour at $95\%$ confidence limit when
$\chi_{>}^{2}\geq 4$. This limit is shown by a red contour in Fig. 4. All
cross sections above this line are excluded at $2\sigma$ significance. By
construction, this does not have an overlap with the best-fit contours at
$2\sigma$.

The detector of the XENON experiment sits in an underground facility about 1.4
km below Earth’s surface. Therefore, if the DM-electron cross section is too
large, the DM particles could lose most of its energy and be stopped in the
earth before it reaches the detector. Following Ref. [63], we perform a simple
energy loss analysis to estimate this. If a DM particle has kinetic energy
$T_{\chi}$ and it scatters with the electrons, then its energy loss rate is
given by

$\dfrac{dT_{\chi}}{dx}=-n_{\mathrm{e}}\int_{0}^{E_{R}^{\mathrm{max}}}\text{d}E_{R}E_{R}\dfrac{d\sigma_{\chi
e}}{dE_{R}}\,.$ (10)

Here, $n_{\mathrm{e}}$ is the electron number density in Earth, $E_{R}$ is the
electron recoil energy in a single collision. This relation can be simplified
by performing the $E_{R}$ integration in the heavy mediator limit,
$d\sigma_{\chi e}/dE_{R}=\sigma_{\chi e}/E_{R}^{\mathrm{max}}$. This leads to
the following expression for the final energy of the DM particle after
traveling a distance $x$,

$T_{\chi}(x)=T_{\chi}^{\text{ini}}\exp\left[-\dfrac{2n_{\mathrm{e}}\sigma_{\chi
e}m_{\mathrm{e}}x}{m_{\chi}}\right]\,.$ (11)

The XENON experiment has a low energy threshold of $\simeq 2{\rm\ keV}$ for
electron recoil [64]. Now Eq.(11) can be inverted to find out the threshold
cross section $\sigma_{\chi e}^{\text{th}}$ for which a DM particle can travel
a distance $x$ before its energy falls below the XENON detector threshold.

Figure 4: Contours in the $m_{\chi}-\sigma_{\chi e}$ plane, depicting the
regions which satisfy the XENON1T excess, as well as exclusion contours. The 1
and $2\sigma$ confidence intervals satisfying the excess are shown in blue
contours, with the best-fit point marked with a white asterisk. The region
above the red line is excluded at $95\%$ confidence level. The dashed gray
lines denote the cross sections for which a DM particle travels a distance of
$100{\rm\ km}$ and $1000{\rm\ km}$, respectively, before its energy falls
below the threshold of XENON1T experiment. We also show the exclusion limits
from SENSEI (dashed blue), EDELWEISS (dashed orange), and PandaX (dashed
green) experiments.

Fig. 4 shows two gray dashed lines for $\sigma_{\chi e}^{\text{th}}$
corresponding to distances $x=100$ km, $1000$ km. A more detailed analysis
would involve integrating $\sigma_{\chi e}^{\text{th}}$ over all energies
weighted by the flux spectrum. However, we assumed $T_{\chi}^{\text{ini}}\sim
1{\rm\ MeV}$ for the sake of simplicity as $\sigma_{\chi e}^{\text{th}}$ has
only logarithmic dependence on $T_{\chi}^{\text{ini}}$. We see that the
traveling distance of the DM particles inside the earth is far greater than
the depth of the XENON experiment for $\sigma_{\chi e}\sim
10^{-30}-10^{-29}{\rm\,cm^{2}}$. Nevertheless, this will imprint a
directionality in the DM flux. The flux coming through the bulk of the earth
will be more attenuated compared to the flux coming from above. In Fig. 4, we
also show the exclusion limits from SENSEI (dashed blue), EDELWEISS (dashed
orange), and PandaX (dashed green) experiments [65, 66, 24].

##  V Daily Modulation of Dark Matter Signal

As the Earth rotates about its axis, the boosted DM flux coming from a
specific direction in the sky has to travel a varying distance through the
Earth. This will cause a daily modulation in the DM scattering event rate in a
detector. Such daily modulation in DM signal has been calculated before, for
example, in Ref.[67]. In this section, we shall provide a rough estimate of
such a daily modulation in the SnBDM signal.

Let us consider a reference frame where the $z$-axis is along the north-south
pole of the Earth. The DM wind is assumed to be coming from a direction with a
latitude $\approx 42^{\circ}$. Therefore, $\theta\approx 48^{\circ}$ is the
corresponding polar angle. The velocity can be written as

$\vec{v}_{\chi}=v_{\chi}(\hat{x}\sin\theta\cos\phi+\hat{y}\sin\theta\sin\phi+\hat{z}\cos\theta)\,.$
(12)

We assume that the DM particles retain this directionality even after
scattering with the DSNB because the interaction mediator is heavier than
typical momentum exchange, resulting in isotropic scattering. We shall ignore
the velocity of the Earth compared to that of the DM particles, which are
boosted to $\sim\mathcal{O}(1{\rm\ MeV})$ energy after scattering. If
$\hat{n}$ denotes the direction of the XENON detector, then

$\hat{n}=\hat{x}\cos{\theta_{d}}\cos{\omega
t}+\hat{y}\cos{\theta_{d}}\sin{\omega t}+\hat{z}\sin{\theta_{d}}\,,$ (13)

where $\theta_{d}$ is the latitude of the detector and $\omega$ is the angular
velocity of the Earth. The angle $\psi$ between the detector and the direction
of DM wind is given as

$\cos\psi=\hat{v}_{\chi}\cdot\hat{n}$ (14)

Using the angle $\psi$, we can find the distance travelled by the DM particle
through the Earth,

$L_{\chi}=(R_{E}-L_{d})\cos\psi+\sqrt{(R_{E}-L_{d})^{2}\cos^{2}\psi-
L_{d}^{2}+2R_{E}L_{d}}\,,$ (15)

where $R_{E}$ is the radius of the Earth and $L_{d}=1.4{\rm\ km}$ is the depth
of the location of the detector.

DM particles lose energy due to scattering with the electrons in the medium
while travelling inside the Earth. When the interaction mediator is heavier
than the typical momentum exchanged (heavy mediator limit), the final energy
of the DM particle after travelling a distance $L_{\chi}$ is given by

$T_{\chi}=T_{\chi}^{\text{ini}}\exp\left[-\frac{2n_{e}\sigma_{\chi
e}m_{\mathrm{e}}L_{\chi}}{m_{\chi}}\right]\,.$ (16)

A daily modulation in $L_{\chi}$ will impart a modulation in DM kinetic energy
$T_{\chi}$, which will in turn imply a 24-hour period variation in the event
rate. In Fig. 5, we show this time variation in the total expected event rate
and DM energy after passing through the Earth. Clearly, one can notice a
significant daily modulation in the event rate, as well as the final DM
energy, and this can be used to put constraints on the DM-electron coupling.

Figure 5: (Left) The total electron-recoil rate variation with time. (Right)
The time variation of the final DM energy at the XENON detector assuming an
initial energy at the surface of the Earth $E_{\mathrm{ini}}=1{\rm\ MeV}$.

##  VI Conclusion & Outlook

The observation of the electronic recoil excess in the XENON1T experiment has
renewed the interest in boosted DM. In this work, we present a scenario where
DM couples to neutrinos and electrons. The Universe is permeated by a constant
source of MeV-energy neutrinos originated from all past supernovae, which may
boost the DM through scattering. This can be a lucrative explanation for the
excess of events observed by the XENON1T experiment. We also perform a
conservative analysis to constrain large DM-electron cross sections without
referring to the excess as such.

DM-electron interactions also lead to a shielding effect as the DM travels
through the earth. The flux is attenuated more when it travels a longer
distance through the earth to reach the detector. This creates a diurnal
modulation in the event rate in the detector [26]. The resulting modulation in
the total event rate is computed. If such a signal is observed, the amplitude
of the modulation could be used to determine the DM interaction strength with
ordinary matter irrespective of the total number of recoil events. This will
be an interesting direction for future investigation.

Furthermore, such BSM interaction with leptons could have other observable
effects. For example, the scattering with DM will change the DSNB energy
spectrum. The neutrinos will lose energy by upscattering the DM particles, and
the DSNB spectrum as shown in Fig. 1 will shift towards lower energy. We found
that change in DSNB would be pronounced if the DM and the mediator masses
$\lesssim\mathcal{O}(\text{few}{\rm\ MeV})$. However, models with such small
DM masses and large interaction cross section will inevitably run into
conflict with BBN. We leave this exercise for a future analysis.

Neutrino telescopes are another class of experiments that are promising to
detect SnBDM . The typical detection threshold for these detectors are above
$\sim 10{\rm\ MeV}$ which is not suitable for detecting the vanilla galactic
halo DM. However, if the DM is boosted to semi-relativistic speeds, then they
can scatter with the electrons in the neutrino detectors and deposit an energy
above threshold. The expected event rate has been computed in a few previous
works and competitive constraints on DM-electron cross section have been
obtained [32]. One way to test the solution presented here would be to compute
the expected number of events in the MiniBooNE or BOREXINO experiment and
compare with the data [68, 69, 70]. The DM particles will also undergo
collisions with the electrons in the sun, and get scattered isotropically in
all directions. As such, the sun will act as a reflector of the DM particles.
This could add extra directionality in the observed DM flux at Earth.

DM-electron scattering before and during recombination would affect the cosmic
microwave background (CMB) anisotropy spectrum. Previous works have
constrained DM-baryon cross section using the CMB data [71, 72]. The bounds
can be translated to the present case by re-scaling the relevant interaction
term with appropriate modifications. We have checked that this limit does not
constrain our best-fit DM mass and cross-section. Finally, we note that such
large DM-electron cross-section as presented in this work, and all other works
on boosted DM, has been excluded by SENSEI, EDELWEISS, and CRESST-II
experiments as shown in Fig. 4 [65, 66, 73]. The requirement of a mediator
lighter than the electron brings up the possibility of its thermalization
before BBN. However, we have explicitly checked that with mediator mass
$m_{\rm med}\lesssim 0.01{\rm\ MeV}$ and electron coupling $g_{e}\lesssim
10^{-7}$, it is possible to reproduce the required DM-electron scattering
cross section, and not violate the BBN constraints at the same time.

Finally, note that neutrino-boosted DM has been studied before in Refs. [40,
24] in the context of solar neutrinos, neutrinos from Sun-like stars, and
cosmic-ray electrons. The constraints arising on DM-electron cross-section
from these studies are typically stronger, due to the larger flux of neutrinos
as compared to the DSNB. On the other hand, the DSNB represents an isotropic
source of neutrinos present from a much older epoch $z\sim 5$. As a result,
the constraints from this study are independent of the solar constraints.

## Acknowledgments

We thank Basudeb Dasgupta and Andre de Gouvea for useful suggestions and
comments on the manuscript. AD was supported by the U.S. Department of Energy
under contract number DE-AC02-76SF00515. MS acknowledges support from the
National Science Foundation, Grant PHY-1630782, and to the Heising-Simons
Foundation, Grant 2017-228.

## References

  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005), arXiv:hep-ph/0404175 .
  * Aprile  _et al._ [2020] E. Aprile, J. Aalbers, F. Agostini, M. Alfonsi, L. Althueser, F. Amaro, V. Antochi, E. Angelino, J. Angevaare, F. Arneodo, and et al., Physical Review D 102, 10.1103/physrevd.102.072004 (2020).
  * Arceo Díaz  _et al._ [2019] S. Arceo Díaz, K.-P. Schröder, K. Zuber, D. Jack, and E. E. Bricio Barrios, arXiv e-prints , arXiv:1910.10568 (2019), arXiv:1910.10568 [astro-ph.SR] .
  * Smirnov and Beacom [2020] J. Smirnov and J. F. Beacom, Phys. Rev. Lett. 125, 131301 (2020), arXiv:2002.04038 [hep-ph] .
  * Bell  _et al._ [2020] N. F. Bell, J. B. Dent, B. Dutta, S. Ghosh, J. Kumar, and J. L. Newstead, Phys. Rev. Lett. 125, 161803 (2020), arXiv:2006.12461 [hep-ph] .
  * Davoudiasl  _et al._ [2020] H. Davoudiasl, P. B. Denton, and J. Gehrlein, Phys. Rev. D 102, 091701 (2020), arXiv:2007.04989 [hep-ph] .
  * Boehm  _et al._ [2020] C. Boehm, D. G. Cerdeno, M. Fairbairn, P. A. N. Machado, and A. C. Vincent, Phys. Rev. D 102, 115013 (2020), arXiv:2006.11250 [hep-ph] .
  * Arcadi  _et al._ [2021] G. Arcadi, A. Bally, F. Goertz, K. Tame-Narvaez, V. Tenorth, and S. Vogl, Phys. Rev. D 103, 023024 (2021), arXiv:2007.08500 [hep-ph] .
  * Borah  _et al._ [2020] D. Borah, S. Mahapatra, D. Nanda, and N. Sahu, Phys. Lett. B 811, 135933 (2020), arXiv:2007.10754 [hep-ph] .
  * Chala and Titov [2020] M. Chala and A. Titov, JHEP 09, 188, arXiv:2006.14596 [hep-ph] .
  * Lindner  _et al._ [2020] M. Lindner, Y. Mambrini, T. B. de Melo, and F. S. Queiroz, Phys. Lett. B 811, 135972 (2020), arXiv:2006.14590 [hep-ph] .
  * Di Luzio  _et al._ [2020] L. Di Luzio, M. Fedele, M. Giannotti, F. Mescia, and E. Nardi, Phys. Rev. Lett. 125, 131804 (2020), arXiv:2006.12487 [hep-ph] .
  * Bally  _et al._ [2020] A. Bally, S. Jana, and A. Trautner, Phys. Rev. Lett. 125, 161802 (2020), arXiv:2006.11919 [hep-ph] .
  * Gao  _et al._ [2020] C. Gao, J. Liu, L.-T. Wang, X.-P. Wang, W. Xue, and Y.-M. Zhong, Phys. Rev. Lett. 125, 131806 (2020), arXiv:2006.14598 [hep-ph] .
  * Bramante and Song [2020] J. Bramante and N. Song, Phys. Rev. Lett. 125, 161805 (2020), arXiv:2006.14089 [hep-ph] .
  * Cai  _et al._ [2020] C. Cai, H. H. Zhang, M. T. Frandsen, M. Rosenlyst, and G. Cacciapaglia, Phys. Rev. D 102, 075018 (2020), arXiv:2006.16267 [hep-ph] .
  * Miranda  _et al._ [2020] O. G. Miranda, D. K. Papoulias, M. Tórtola, and J. W. F. Valle, Phys. Lett. B 808, 135685 (2020), arXiv:2007.01765 [hep-ph] .
  * Karozas  _et al._ [2021] A. Karozas, S. F. King, G. K. Leontaris, and D. K. Papoulias, Phys. Rev. D 103, 035019 (2021), arXiv:2008.03295 [hep-ph] .
  * Zu  _et al._ [2021] L. Zu, R. Foot, Y.-Z. Fan, and L. Feng, JCAP 01, 070, arXiv:2007.15191 [hep-ph] .
  * Farzan and Rajaee [2020] Y. Farzan and M. Rajaee, Phys. Rev. D 102, 103532 (2020), arXiv:2007.14421 [hep-ph] .
  * Khan [2020] A. N. Khan, Phys. Lett. B 809, 135782 (2020), arXiv:2006.12887 [hep-ph] .
  * Dutta  _et al._ [2021] K. Dutta, A. Ghosh, A. Kar, and B. Mukhopadhyaya, (2021), arXiv:2103.14664 [hep-ph] .
  * McKeen  _et al._ [2020] D. McKeen, M. Pospelov, and N. Raj, Phys. Rev. Lett. 125, 231803 (2020), arXiv:2006.15140 [hep-ph] .
  * Jho  _et al._ [2021] Y. Jho, J.-C. Park, S. C. Park, and P.-Y. Tseng, (2021), arXiv:2101.11262 [hep-ph] .
  * Kannike  _et al._ [2020] K. Kannike, M. Raidal, H. Veermäe, A. Strumia, and D. Teresi, Phys. Rev. D 102, 095002 (2020), arXiv:2006.10735 [hep-ph] .
  * Fornal  _et al._ [2020] B. Fornal, P. Sandick, J. Shu, M. Su, and Y. Zhao, Physical Review Letters 125, 10.1103/physrevlett.125.161804 (2020).
  * Alhazmi  _et al._ [2020] H. Alhazmi, D. Kim, K. Kong, G. Mohlabeng, J.-C. Park, and S. Shin, (2020), arXiv:2006.16252 [hep-ph] .
  * D’Eramo and Thaler [2010] F. D’Eramo and J. Thaler, JHEP 06, 109, arXiv:1003.5912 [hep-ph] .
  * Agashe  _et al._ [2014] K. Agashe, Y. Cui, L. Necib, and J. Thaler, JCAP 10, 062, arXiv:1405.7370 [hep-ph] .
  * Berger  _et al._ [2015] J. Berger, Y. Cui, and Y. Zhao, JCAP 02, 005, arXiv:1410.2246 [hep-ph] .
  * Ema  _et al._ [2019] Y. Ema, F. Sala, and R. Sato, Phys. Rev. Lett. 122, 181802 (2019), arXiv:1811.00520 [hep-ph] .
  * Bringmann and Pospelov [2019] T. Bringmann and M. Pospelov, Phys. Rev. Lett. 122, 171801 (2019), arXiv:1810.10543 [hep-ph] .
  * Cappiello  _et al._ [2019] C. V. Cappiello, K. C. Y. Ng, and J. F. Beacom, Phys. Rev. D 99, 063004 (2019), arXiv:1810.07705 [hep-ph] .
  * Yin [2019] W. Yin, EPJ Web Conf. 208, 04003 (2019), arXiv:1809.08610 [hep-ph] .
  * McKeen and Raj [2019] D. McKeen and N. Raj, Phys. Rev. D 99, 103003 (2019), arXiv:1812.05102 [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] .
  * Ema  _et al._ [2021] Y. Ema, F. Sala, and R. Sato, SciPost Phys. 10, 072 (2021), arXiv:2011.01939 [hep-ph] .
  * Olivares-Del Campo  _et al._ [2018] A. Olivares-Del Campo, C. Bœhm, S. Palomares-Ruiz, and S. Pascoli, Phys. Rev. D 97, 075039 (2018), arXiv:1711.05283 [hep-ph] .
  * Ghosh  _et al._ [2020] S. Ghosh, R. Khatri, and T. S. Roy, Phys. Rev. D 102, 123544 (2020), arXiv:1908.09843 [hep-ph] .
  * Zhang [2020] Y. Zhang, (2020), arXiv:2001.00948 [hep-ph] .
  * Fox and Poppitz [2009] P. J. Fox and E. Poppitz, Physical Review D 79, 10.1103/physrevd.79.083528 (2009).
  * Bernabei  _et al._ [2008] R. Bernabei, P. Belli, F. Montecchia, F. Nozzoli, F. Cappella, A. Incicchitti, D. Prosperi, R. Cerulli, C. J. Dai, H. L. He, and et al., Physical Review D 77, 10.1103/physrevd.77.023506 (2008).
  * Chang  _et al._ [2014] S. Chang, R. Edezhath, J. Hutchinson, and M. Luty, Physical Review D 90, 10.1103/physrevd.90.015011 (2014).
  * Beacom [2010] J. F. Beacom, Ann. Rev. Nucl. Part. Sci. 60, 439 (2010), arXiv:1004.3311 [astro-ph.HE] .
  * Lunardini [2016] C. Lunardini, Astropart. Phys. 79, 49 (2016), arXiv:1007.3252 [astro-ph.CO] .
  * De Gouvêa  _et al._ [2020] A. De Gouvêa, I. Martinez-Soler, Y. F. Perez-Gonzalez, and M. Sen, Phys. Rev. D 102, 123012 (2020), arXiv:2007.13748 [hep-ph] .
  * Hopkins and Beacom [2006] A. M. Hopkins and J. F. Beacom, Astrophys. J. 651, 142 (2006), arXiv:astro-ph/0601463 [astro-ph] .
  * Horiuchi  _et al._ [2009] S. Horiuchi, J. F. Beacom, and E. Dwek, Phys. Rev. D 79, 083013 (2009), arXiv:0812.3157 [astro-ph] .
  * Aghanim  _et al._ [2018] N. Aghanim _et al._ (Planck), (2018), arXiv:1807.06209 [astro-ph.CO] .
  * Zhang  _et al._ [2015] H. Zhang _et al._ (Super-Kamiokande), Astropart. Phys. 60, 41 (2015), arXiv:1311.3738 [hep-ex] .
  * Wolfenstein [1978] L. Wolfenstein, Phys. Rev. D 17, 2369 (1978).
  * Mikheev and Smirnov [1985] S. P. Mikheev and A. Yu. Smirnov, Sov. J. Nucl. Phys. 42, 913 (1985), [Yad. Fiz.42,1441(1985)].
  * Duan  _et al._ [2006] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys. Rev. D74, 105014 (2006), arXiv:astro-ph/0606616 [astro-ph] .
  * Hannestad  _et al._ [2006] S. Hannestad, G. G. Raffelt, G. Sigl, and Y. Y. Y. Wong, Phys. Rev. D74, 105010 (2006), [Erratum: Phys. Rev.D76,029901(2007)], arXiv:astro-ph/0608695 [astro-ph] .
  * Dasgupta and Dighe [2007] B. Dasgupta and A. Dighe, Phys. Rev. D 75, 093002 (2007), arXiv:hep-ph/0510219 .
  * Falkowski  _et al._ [2009] A. Falkowski, J. Juknevich, and J. Shelton, Dark matter through the neutrino portal (2009), arXiv:0908.1790 [hep-ph] .
  * Lindner  _et al._ [2010] M. Lindner, A. Merle, and V. Niro, Physical Review D 82, 10.1103/physrevd.82.123529 (2010).
  * Macias and Wudka [2015] V. G. Macias and J. Wudka, Effective theories for dark matter interactions and the neutrino portal paradigm (2015), arXiv:1506.03825 [hep-ph] .
  * Blennow  _et al._ [2019] M. Blennow, E. Fernandez-Martinez, A. Olivares-Del Campo, S. Pascoli, S. Rosauro-Alcaraz, and A. V. Titov, The European Physical Journal C 79, 10.1140/epjc/s10052-019-7060-5 (2019).
  * Fileviez Pérez  _et al._ [2019] P. Fileviez Pérez, C. Murgui, and A. D. Plascencia, Phys. Rev. D 100, 035041 (2019), arXiv:1905.06344 [hep-ph] .
  * Navarro  _et al._ [1996] J. F. Navarro, C. S. Frenk, and S. D. M. White, The Astrophysical Journal 462, 563 (1996).
  * Ackermann  _et al._ [2015] M. Ackermann _et al._ (Fermi-LAT), JCAP 09, 008, arXiv:1501.05464 [astro-ph.CO] .
  * Starkman  _et al._ [1990] G. D. Starkman, A. Gould, R. Esmailzadeh, and S. Dimopoulos, Phys. Rev. D 41, 3594 (1990).
  * Aprile  _et al._ [2019] E. Aprile _et al._ (XENON), Phys. Rev. D 99, 112009 (2019), arXiv:1902.11297 [physics.ins-det] .
  * Barak  _et al._ [2020] L. Barak _et al._ (SENSEI), Phys. Rev. Lett. 125, 171802 (2020), arXiv:2004.11378 [astro-ph.CO] .
  * Arnaud  _et al._ [2020] Q. Arnaud _et al._ (EDELWEISS), Phys. Rev. Lett. 125, 141301 (2020), arXiv:2003.01046 [astro-ph.GA] .
  * Kouvaris and Shoemaker [2014] C. Kouvaris and I. M. Shoemaker, Phys. Rev. D 90, 095011 (2014), arXiv:1405.1729 [hep-ph] .
  * Aguilar-Arevalo  _et al._ [2018] A. Aguilar-Arevalo, M. Backfish, A. Bashyal, B. Batell, B. Brown, R. Carr, A. Chatterjee, R. Cooper, P. deNiverville, R. Dharmapalan, and et al., Physical Review D 98, 10.1103/physrevd.98.112004 (2018).
  * Aguilar-Arevalo  _et al._ [2017] A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. Lett. 118, 221803 (2017), arXiv:1702.02688 [hep-ex] .
  * Bellini  _et al._ [2014] G. Bellini, J. Benziger, D. Bick, G. Bonfini, D. Bravo, M. Buizza Avanzini, B. Caccianiga, L. Cadonati, F. Calaprice, P. Cavalcante, and et al., Physical Review D 89, 10.1103/physrevd.89.112007 (2014).
  * Gluscevic and Boddy [2018] V. Gluscevic and K. K. Boddy, Physical Review Letters 121, 10.1103/physrevlett.121.081301 (2018).
  * Xu  _et al._ [2018] W. L. Xu, C. Dvorkin, and A. Chael, Physical Review D 97, 10.1103/physrevd.97.103530 (2018).
  * Abdelhameed  _et al._ [2019] A. H. Abdelhameed _et al._ (CRESST), Phys. Rev. D 100, 102002 (2019), arXiv:1904.00498 [astro-ph.CO] .