# Nonparametric Feature Selection by Random Forests and Deep Neural Networks

Xiaojun Mao
School of Mathematical Sciences,
Shanghai Jiao Tong University, Shanghai, China

Liuhua Peng
School of Mathematics and Statistics,
The University of Melbourne, Victoria, Australia
and
Zhonglei Wang
Wang Yanan Institute for Studies in Economics,
Xiamen University, Xiamen, Fujian, China Correspondence to<EMAIL_ADDRESS>
###### Abstract

Random forests are a widely used machine learning algorithm, but their
computational efficiency is undermined when applied to large-scale datasets
with numerous instances and useless features. Herein, we propose a
nonparametric feature selection algorithm that incorporates random forests and
deep neural networks, and its theoretical properties are also investigated
under regularity conditions. Using different synthetic models and a real-world
example, we demonstrate the advantage of the proposed algorithm over other
alternatives in terms of identifying useful features, avoiding useless ones,
and the computation efficiency. Although the algorithm is proposed using
standard random forests, it can be widely adapted to other machine learning
algorithms, as long as features can be sorted accordingly.

Key words: Feature importance, Maximum mean discrepancy, Reproducing kernel
Hilbert space.

##  1 Introduction

Random forests (Breiman, 2001, RFs) are a widely used machine learning
algorithm (Caruana and Niculescu-Mizil, 2006; Criminisi et al., 2012;
Fernández-Delgado et al., 2014). However, their computational efficiency is
compromised when they are applied to large-scale datasets with numerous
useless features.

Since the landmark of Breiman (2001), research on RFs has been active in
different scientific fields. Xu and Jelinek (2005) used RFs for structured
language learning, and they showed that their method outperformed its
competitors in terms of perplexity and error rates. Based on a synthetic
dataset generated from a certain reference distribution, Shi and Horvath
(2006) applied RFs to obtain dissimilarities among the original unlabeled
data. Payet and Todorovic (2010) proposed an RFs algorithm to directly
estimate the ratio of the proposal and posterior distributions
nonparametrically for the Metropolis-Hastings algorithm, and the corresponding
theoretical properties were also investigated. Using a Mondrian process,
Lakshminarayanan et al. (2014) proposed a computationally efficient RFs
algorithm for online learning, but splits were made independent of the
response of interest. Inspired by the local polynomial regression, Li and
Martin (2017) proposed a robust RFs algorithm incorporating different loss
functions and showed that their method generalized the standard and quantile
RFs (Meinshausen, 2006). Haghiri et al. (2018) proposed a comparison-based RFs
algorithm for the case when the sample was not representative and when it was
difficult to measure the distance between instances. Siblini et al. (2018)
proposed an algorithm for extreme multi-label learning using a tree-based
method and demonstrated that the computation was more efficient than other
competitors in parallel. Scornet et al. (2015), Mentch and Hooker (2016), and
Wager and Athey (2018) investigated theoretical properties of RFs algorithms.
Refer to Criminisi et al. (2012) and Goel et al. (2017) for a more
comprehensive review of RFs.

Albeit it is common to have numerous features in practice, only a limited
portion contributes to the response of interest (Fan and Lv, 2008). Because
existing RFs algorithms do not identify useful features before growing trees,
the corresponding computation efficiency is undermined, especially when most
features are useless. For example, it is well known that local polynomial
regression suffers from the curse of dimensionality; thus, the estimation
efficiency of Li and Martin (2017) is questionable. In addition, the online
RFs (Lakshminarayanan et al., 2014) may also lead to inefficient estimation
when the number of useless features is large, since splits are made based on
features regardless of the response of interest. Thus, feature selection is
essential for a high dimensional dataset.

Even though RFs have been actively investigated under different scientific
fields, feature selection by RFs does not receive so much attention. There are
two main approaches for feature selection by RFs. One approach is permutation-
based, using feature importance (Breiman, 2001). Strobl et al. (2008) proposed
a permutation algorithm to compute conditional feature importance, but they
did not provide a general guidance for feature selection. Kursa and Rudnicki
(2010) used a new set of shadow features to debias the feature importance, and
a feature selection procedure was also proposed based on a “Z score”; also see
Sandri and Zuccolotto (2008) for a similar approach. Altmann et al. (2010)
proposed to permute the response vector to get the “null importance”, and it
was used for feature selection heuristically. Genuer et al. (2010) proposed a
two-step procedure for feature section based on feature importance, and their
method worked well for highly correlated features. However, the theoretical
properties of existing permutation-based methods have not been rigorously
investigated. The other approach is based on the minimum depth proposed by
Ishwaran et al. (2010), and it has been widely applied in survival analysis
(Ishwaran et al., 2008; Twyman-Saint Victor et al., 2015; Benci et al., 2016).
The corresponding theoretical properties were investigated by Ishwaran et al.
(2010) under strong conditions. For example, the split for a node should
always be the median of the values, but such an assumption may undermine the
estimation efficiency.

Herein, we propose a nonparametric feature selection algorithm that
incorporates RFs and deep neural networks (NFSRD). Specifically, we adopt
nonparametric two-sample tests using deep neural networks (Liu et al., 2020)
to select _useful_ features, and the corresponding theoretical properties are
investigated under regularity conditions. To improve the computational
efficiency of the NFSRD, subsampling is adopted. Experiments reveal that the
NFSRD outperforms its alternatives in terms of detecting useful features,
avoiding useless ones, and the computation efficiency. Another advantage of
subsampling is that it saves the computer memory while retaining the desired
accuracy for feature selection; refer to Section 4 for details.

The NFSRD differs from existing works in the following aspects. First, we
propose the use of nonparametric two-sample tests to select useful features
incorporating RFs and deep neural networks, and the corresponding theoretical
properties can be rigorously established under regularity conditions. In
addition, subsampling is adopted for feature selection; thus, the NFSRD is
more computationally efficient than existing RFs-based feature selection
algorithms. Moreover, we do not make any strong assumption for splits, as in
Ishwaran et al. (2010). Rather, the only crucial assumption we make is that a
limited portion of the features contributes to the response of interest, and
such an assumption is also widely adopted for feature selection (Fan and Lv,
2008). Besides, the NFSRD adopts a forward feature selection procedure based
on the features sorted by their importance, so it can be widely adapted to
existing algorithms mentioned in the preceding paragraph and other machine
learning algorithms (Chen and Guestrin, 2016; Wager and Athey, 2018), as long
as the features can be sorted accordingly.

The remainder of this paper is organized as follows. The model setup is
introduced in Section 2. The detailed algorithm for the NFSRD is presented in
Section 3, and the corresponding theoretical properties are also investigated.
Simulation studies are presented in Section 4. Section 5 describes an
application of the NFSRD to identify useful features based on a
superconductivity dataset. The conclusions are provided in Section 6.

##  2 Model Setup

Let $\bm{x}_{i}=(x_{i1},\ldots,x_{ip})\in\mathcal{X}$ be a $p$-dimensional
feature vector and $y_{i}\in\mathbb{R}$ be the corresponding response of
interest, where $i\in\mathcal{I}_{n}$, $\mathcal{I}_{n}=\\{1,\ldots,n\\}$ is
the sample index set of size $n$, and $\mathcal{X}\subset\mathbb{R}^{p}$.
Consider the following regression model:

$y_{i}=f(\bm{x}_{i})+\epsilon_{i}\quad(i\in\mathcal{I}_{n}),$ (1)

where $f(\bm{x}_{i})$ is a smooth function involving $K_{0}\ll p$ useful
features (Fan and Lv, 2008), and $\epsilon_{i}$ is white noise. We are
interested in identifying the $K_{0}$ useful features based on the sample
$\mathcal{D}_{n}=\\{(\bm{x}_{i},y_{i}):i\in\mathcal{I}_{n}\\}$. Without loss
of generality, lower cases denote observed data, and upper cases denote the
associated random variables. The vectors are column-wise, unless explicitly
explained otherwise.

Before presenting the NFSRD algorithm, we briefly introduce a bagged-tree
learner and RFs. A bagged-tree learner is

$\widehat{\theta}_{TR}(\bm{x})=\mathbb{E}_{*}\\{t(\bm{x};\bm{Z}^{*}_{1},\ldots,\bm{Z}^{*}_{n})\\},$

where $\bm{z}_{i}=(\bm{x}_{i},y_{i})$, and
$\\{\bm{Z}^{*}_{i}:i\in\mathcal{I}_{n}\\}$ is a random sample of size $n$
generated by an empirical distribution
$\mathbb{P}^{*}(\bm{z})=n^{-1}\sum_{i\in\mathcal{I}_{n}}\mathds{I}{(\bm{z}=\bm{z}_{i})}$,
$\mathds{I}{(\bm{z}=\bm{z}_{i})}=1$ if $\bm{z}=\bm{z}_{i}$ and 0 otherwise,
and $\mathbb{E}_{*}(\cdot)$ is the expectation with respect to
$\mathbb{P}^{*}(\bm{z})$. For ease of notation, the sample size $n$ is omitted
for $\widehat{\theta}_{TR}(\bm{x})$ and other statistics.

Notably, RFs extend bagged-tree learners by allowing additional randomness
within trees to reduce the correlation among them. Specifically, an RFs
learner is

$\widehat{\theta}_{RF}(\bm{x})=\mathbb{E}_{*}\\{t(\bm{x};\xi,\bm{Z}^{*}_{1},\ldots,\bm{Z}^{*}_{n})\\},$
(2)

where $\xi\sim\Xi$, and $\Xi$ is a pre-specified distribution; see Breiman
(2001) for details. A popular choice for $\Xi$ is the random selection of
candidate features for the split of each node (Hastie et al., 2009, §15.2).
The RFs learner $\widehat{\theta}_{RF}(\bm{x})$ is a bagged-tree learner
$\widehat{\theta}_{TR}(\bm{x})$ if $\Xi$ is omitted. Bootstrapping is used to
approximate $\widehat{\theta}_{RF}(\bm{x})$ in (2) by

$\widehat{\theta}^{(B)}_{RF}(\bm{x})=\frac{1}{B}\sum_{b=1}^{B}t(\bm{x};\xi_{b},\bm{z}^{*}_{b1},\ldots,\bm{z}^{*}_{bn}),$
(3)

where $B$ is the number of bootstraps, $\xi_{b}{\sim}\Xi$, and
$\bm{z}^{*}_{bi}$ is a “realization” of $\bm{Z}^{*}_{i}$ for $b=1,\ldots,B$
and $i\in\mathcal{I}_{n}$. In practice, a large $B$ is suggested to make the
approximation error between (2) and (3) negligible .

##  3 NFSRD

The assumption $K_{0}\ll p$ validates the feature selection (Guyon and
Elisseeff, 2003). The NFSRD consists of two steps. The first step corresponds
to obtaining bias-corrected feature importance (BCFI) by shadow features, and
the second step is a forward feature selection based on the _ordered_ features
by BCFI. Two-sample tests using maximum mean discrepancy and deep neural
networks (MMD-D) are conducted sequentially for feature selection in the
second step. To improve the computational efficiency, subsampling is applied
in both steps.

Before diving into details, we briefly discuss the intuition for the NFSRD.
The proposed method is a forward-stepwise selection algorithm. To avoid
including useless features as much as possible, we first order the features
based on their “importance” for estimating the response of interest. Due to
their flexibility for nonparametric modeling, RFs are implemented to order the
features by a subset of the instances. Based on the ordered features, forward-
stepwise selection is conducted by sequential hypothesis tests. After fitting
a full model and a reduced model only retaining the first several important
features nonparametrically, the corresponding null hypothesis is that the
distributions of the residuals from a full model and a reduced model are
identical. If the null hypothesis holds for a certain reduced model, then we
treat the involved features as useful and the remaining ones as useless. We
still use RFs to train the full and reduced models, and other nonparametric
algorithms can be implemented as long as certain consistency results hold; see
Supplementary Material and Theorem 1 of Mentch and Hooker (2016) for details.

###  3.1 Bias-Corrected Feature Importance

Feature importance is a rudimentary indicator of the usefulness of features,
and it serves as a building block of the NFSRD.

Existing algorithms calculate the feature importance using the entire sample
$\mathcal{D}_{n}$. To achieve better computational efficiency, we propose
using subsampling. That is, the feature importance is obtained based on a
subsample $\\{\bm{z}_{i}:i\in\mathcal{A}_{0}\\}$, where $\mathcal{A}_{0}$ is a
subset of $\mathcal{I}_{n}$, and its size is $m_{0}<n$. The numerical results
reveal that $m_{0}$ should be large to guarantee good performance, and we
suggest $m_{0}\geq 400$ for practical guidance; see Section 4 for details. For
example, we can evaluate the feature importance by

$F_{k}=\frac{1}{B}\sum_{t_{b}}\sum_{j:s_{bj}=k}\Delta(s_{bj},t_{b})\quad(k=1,\ldots,p),$
(4)

where the first summation is with respect to the $B$ trees, the second
summation is with respect the splits made for the $b$th tree, $s_{bj}$ is the
feature used to split the $j$th node of $t_{b}$,
$\Delta(s_{bj},t_{b})=w_{bj}V_{bj}-w^{(l)}_{bj}V^{(l)}_{bj}-w^{(r)}_{bj}V^{(r)}_{bj}$
is the weighted decrease in variance, $w_{bj}=n_{bj}/m_{0}$, $n_{bj}$ is the
number of instances in the $j$th node of $t_{b}$, $w^{(l)}_{bj}$ and
$w^{(r)}_{bj}$ are the corresponding proportions on the left and right
subnodes after splitting the $j$th node by $s_{bj}$, $V_{bj}$ is the sample
variance of the response of interest in the $j$th node, and $V^{(l)}_{bj}$ and
$V^{(r)}_{bj}$ are variances of the two subnodes. See Breiman (2001) and
Sandri and Zuccolotto (2008) for details.

However, the feature importance in (4) is unfairly biased toward those with
numerous distinct values (White and Liu, 1994; Louppe et al., 2013). To
debias, we incorporate shadow features (Sandri and Zuccolotto, 2008; Kursa and
Rudnicki, 2010). Specifically, for $i\in\mathcal{A}_{0}$, a shadow feature
$\bm{x}^{\dagger}_{i}$ is randomly selected from
$\\{\bm{x}_{i}:i\in\mathcal{A}_{0}\\}$ without replacement. Thereafter, an RF
is trained based on the extended data
$\\{\bm{z}_{i}^{\dagger}:i\in\mathcal{A}_{0}\\}$ with
$\bm{z}_{i}^{\dagger}=(\bm{x}_{i},\bm{x}^{\dagger}_{i},y_{i})$; see Sandri and
Zuccolotto (2008) for details. The entire procedure is repeated $R$ times, and
the BCFI of $X_{k}$ is calculated as

$I_{k}=\frac{1}{R}\sum_{r=1}^{R}(F_{rk}-F_{rk^{\dagger}})\quad(k=1,\ldots,p),$
(5)

where $F_{rk}$ and $F_{rk^{\dagger}}$ correspond to $X_{k}$ and
$X_{k}^{\dagger}$, respectively, and $X_{k}^{\dagger}$ is the corresponding
shadow feature for $X_{k}$ for the $r$th repetition; Sandri and Zuccolotto
(2008) suggest that $R=100$ in practice. Algorithm 1 shows the algorithm for
BCFI.

Algorithm 1 Bias-Corrected Feature Importance (BCFI)

Input: $\mathcal{D}_{n}$, $m_{0}$, and $R$.

Select $\mathcal{A}_{0}$ randomly from $\mathcal{I}_{n}$.

for $r=1,\ldots,R$ do

for $i\in\mathcal{A}_{0}$ do

Generate $\bm{x}_{i}^{\dagger}$ from $\\{\bm{x}_{i}:i\in\mathcal{A}_{0}\\}$
randomly without replacement.

end for

Train an RF using
$\\{(\bm{x}_{i},\bm{x}_{i}^{\dagger},y_{i}):i\in\mathcal{A}_{0}\\}$.

Solicit $\\{(F_{rk},F_{rk^{\dagger}}):k=1,\ldots,p\\}$.

for $k=1,\ldots,p$ do

Calculate $I_{rk}=F_{rk}-F_{rk^{\dagger}}$.

end for

end for

for $k=1,\ldots,p$ do

Calculate $I_{k}=R^{-1}\sum_{r=1}^{R}I_{rk}$

end for

Output: $\\{I_{k}:k=1,\ldots,p\\}$.

######  Remark 1.

We should be cautious about algorithms using feature importance as the _only_
criterion for feature selection, especially as no rigorous theoretical
properties have been investigated. Thus, in the next subsection, we propose
the use of nonparametric two-sample tests for feature selection by ordered
features using BCFI.

######  Remark 2.

Other than the shadow features, another way to debias the feature importance
is to use the out of bag sample (OOB). Li et al. (2019) argued that the OOB-
based feature importance outperforms the traditional one in terms of AUC
scores. In addition, we can also consider the minimum depth (Ishwaran et al.,
2010) as an importance metric; see Section 4 for details. Thus, those two
metrics can be applied to obtain the feature importance instead.

###  3.2 Feature Selection by Deep Neural Network

Ideally, we want to identify a set of features, say $\bm{X}^{(K_{0})}$, such
that

$Y\perp\bm{X}^{(-K_{0})}\mid\bm{X}^{(K_{0})},$

where $X\perp Y\mid Z$ denotes that $X$ and $Y$ are conditionally independent
given $Z$, and $\bm{X}^{(-K_{0})}$ contains features other than
$\bm{X}^{(K_{0})}$. Equivalently, we want to identify $\bm{x}^{(K_{0})}$ such
that

$f(\bm{x})=f(\bm{x}^{(K_{0})})\quad(\bm{x}\in\mathcal{X}).$ (6)

Based on the ordered features by BCFI, a forward feature selection algorithm
is proposed by sequentially conducting nonparametric two-sample tests using
deep neural networks (FS-D). Let the ordered BCFI be $I^{(1)}\geq
I^{(2)}\geq\cdots\geq I^{(p)}$ and $x_{i}^{(k)}$ be the feature corresponding
to $I^{(k)}$ for $k=1,\ldots,p$. Consider the following hypothesis testing
problem:

$H_{0K}:f(\bm{x})=f(\bm{x}^{(K)})\quad(\bm{x}\in\mathcal{X}),$ (7)

where $\bm{x}^{(K)}=(x^{(1)},\ldots,x^{(K)})$ contains the first $K$ features
with respect to the ordered BCFI, and $K=1,\ldots,p$.

Although Scornet et al. (2015) investigated theoretical properties of an RFs
estimator, the convergence rate of $g(\bm{x})-f(\bm{x})$ is still an open
question, where $g(\bm{x})=\mathbb{E}\\{\widehat{\theta}_{RF}(\bm{x})\\}$.
Thus, it is difficult to work with (7) directly. Instead, we consider
$\eta=y-g(\bm{x})$ and $\eta^{(K)}=y-g(\bm{x}^{(K)})$, where
$g(\bm{x}^{(K)})=\mathbb{E}\\{\widehat{\theta}_{RF}(\bm{x}^{(K)})\\}$. Thus,
we have

$\displaystyle\eta^{(K)}$ $\displaystyle=$ $\displaystyle
y-g(\bm{x})+[g(\bm{x})-g(\bm{x}^{(K)})]$ $\displaystyle=$
$\displaystyle\eta+[g(\bm{x})-g(\bm{x}^{(K)})].$

Instead of $H_{0K}$ in (7), we consider

$H_{0K}^{\prime}:\mathbb{P}=\mathbb{P}^{(K)},$ (8)

where $\mathbb{P}$ and $\mathbb{P}^{(K)}$ are the distributions of $\eta$ and
$\eta^{(K)}$, respectively. If a feature is useful, it also contributes to
$g(\bm{x})$, and vice versa. Thus, it is valid to use (8) for feature
selection.

######  Remark 3.

Because the standard RFs (Breiman, 2001) is widely implemented in practice, we
apply it to obtain BCFI and approximate $g(\bm{x})$ and $g(\bm{x}^{(K)})$
using (3). If we consider an honest tree (Wager and Athey, 2018), we can
directly use $g(\bm{x})$ to approximate $f(\bm{x})$ according to the
consistency result in Theorem 3 of Wager and Athey (2018).

To test $H_{0K}^{\prime}$ in (8), we adopt a nonparametric two-sample test
using the MMD-D (Liu et al., 2020). The maximum mean discrepancy between
$\mathbb{P}$ and $\mathbb{P}^{(K)}$ is

$\displaystyle\mbox{MMD}(\mathbb{P},\mathbb{P}^{(K)};\mathcal{H}_{\kappa})$
$\displaystyle=$ $\displaystyle\sup_{f\in\mathcal{H}_{\kappa};\lVert
f\rVert_{\mathcal{H}_{\kappa}}\leq 1}\lvert
E\\{f(\bm{X})\\}-E\\{f(\bm{X}^{(K)})\\}\rvert$ $\displaystyle=$
$\displaystyle\\{E[\kappa(\eta,\eta)]+E[\kappa(\eta^{(K)},\eta^{(K)})]-2E[\kappa(\eta,\eta^{(K)})]\\}^{1/2},$

where $\kappa:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ is the kernel for a
reproducing kernel Hilbert space (RKHS) $\mathcal{H}_{\kappa}$, $\lVert
f\rVert_{\mathcal{H}_{\kappa}}$ is the corresponding norm; see Supplementary
Material for a brief introduction on RKHS. For a characteristic kernel
$\kappa$ (Fukumizu et al., 2007; Gretton et al., 2012),
$\mbox{MMD}(\mathbb{P},\mathbb{P}^{(K)};\mathcal{H}_{\kappa})=0$ is equivalent
to $\mathbb{P}=\mathbb{P}^{(K)}$. Recall that a kernel $\kappa$ is
characteristic if the map $\mathbb{Q}\to m_{\mathbb{Q}}$ is one-to-one, where
$m_{\mathbb{Q}}=\mathbb{E}_{\zeta\sim\mathbb{Q}}\kappa(\cdot,\zeta)\in\mathcal{H}_{\kappa}$
for $\mathbb{Q}\in\mathcal{P}_{\mathbb{R}}$, $\mathcal{P}_{\mathbb{R}}$ is the
set of probability measures on the measurable space
$(\mathbb{R},\mathcal{B})$, $\mathcal{B}$ is the Borel $\sigma$-algebra on
$\mathbb{R}$, and $\mathbb{E}_{\zeta\sim\mathbb{Q}}\kappa(\cdot,\zeta)$ is the
expectation of $\kappa(\cdot,\zeta)$ with respect to $\zeta\sim\mathbb{Q}$.

######  Remark 4.

The test statistic is based on (LABEL:eq:_test_stat), and its intuition is
briefly discussed. Evidently, the value
$\mbox{MMD}(\mathbb{P},\mathbb{P}^{(K)};\mathcal{H})$ is determined by the
functional space $\mathcal{H}$. On the one hand, the functional space
$\mathcal{H}$ should be sufficiently large to distinguish two different
distributions $\mathbb{P}$ and $\mathbb{P}^{(K)}$ by the supremum of $\lvert
E\\{f(\bm{X})\\}-E\\{f(\bm{X}^{(K)})\\}\rvert$ for $f\in\mathcal{H}$. On the
other hand, the functional space $\mathcal{H}$ should also be restricted, such
that the estimator of $\mbox{MMD}(\mathbb{P},\mathbb{P}^{(K)};\mathcal{H})$
converges duly to guarantee good statistical properties. Thus, a unit ball of
an RKHS, $\\{f\in\mathcal{H}_{\kappa};\lVert
f\rVert_{\mathcal{H}_{\kappa}}\leq 1\\}$, is a good choice for $\mathcal{H}$.
See Gretton et al. (2012) for details on (LABEL:eq:_test_stat).

If $\\{\eta_{i}:i\in\mathcal{A}_{1}\\}$ and
$\\{\eta_{i}^{(K)}:i\in\mathcal{A}_{2}\\}$ were observed, an estimator of
$\mbox{MMD}(\mathbb{P},\mathbb{P}^{(K)};\mathcal{H}_{\kappa})$ could be
obtained by a U-statistic based on two subsamples $\mathcal{A}_{1}$ and
$\mathcal{A}_{2}$:

$\widehat{\mbox{MMD}}_{u}^{2}(\widehat{\mathbb{P}},\widehat{\mathbb{P}}^{(K)};\kappa)=\frac{1}{m_{1}(m_{1}-1)}\sum_{i\neq
j}H_{ij}^{(K)},$ (10)

where
$\widehat{\mathbb{P}}(\eta)=m_{1}^{-1}\sum_{i\in\mathcal{A}_{1}}\mathds{I}{(\eta=\eta_{j})}$
is the empirical distribution of $\\{\eta_{i}:i\in\mathcal{A}_{1}\\}$,
$\widehat{\mathbb{P}}^{(K)}$ is the one of
$\\{\eta_{i}^{(K)}:i\in\mathcal{A}_{2}\\}$, $m_{1}$ is the size of both
$\mathcal{A}_{1}$ and $\mathcal{A}_{2}$,
$\mathcal{A}_{1}\cap\mathcal{A}_{2}=\emptyset$, and
$H_{ij}^{(K)}=\kappa(\eta_{i},\eta_{j})+\kappa(\eta_{i}^{(K)},\eta_{j}^{(K)})-\kappa(\eta_{i},\eta_{j}^{(K)})-\kappa(\eta_{i}^{(K)},\eta_{j})$.
The disjoint condition between $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$
guarantees independence between the two error sets. The numerical results
reveal that the sample sizes of $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ should
be large, and we suggest that these sizes be larger than 400 for practical
guidance ; see Section 4 for details.

However, $\eta_{i}$ and $\eta_{i}^{(K)}$ are unavailable; thus, we use
$\widehat{\eta}_{ni}=y_{i}-\widehat{\theta}_{RF}^{(B)}(\bm{x}_{i})$ for
$i\in\mathcal{A}_{1}$ and
$\widehat{\eta}_{ni}^{(K)}=y_{i}-\widehat{\theta}_{RF}^{(B)}(\bm{x}_{i}^{(K)})$
$i\in\mathcal{A}_{2}$ instead. The following theorem validates this choice:

######  Theorem 1.

Under mild conditions, in Supplementary Material,

$\widetilde{\mbox{MMD}}_{u}^{2}(\widetilde{\mathbb{P}},\widetilde{\mathbb{P}}^{(K)};\kappa)=\frac{1}{m_{1}(m_{1}-1)}\sum_{i\neq
j}\widehat{H}_{ij}^{(K)}$

has the same limiting distribution as
$\widehat{\mbox{MMD}}_{u}^{2}(\mathbb{P},\mathbb{P}^{(K)};k)$ in (10), where
$\widetilde{\mathbb{P}}$ and $\widetilde{\mathbb{P}}^{(K)}$ are the empirical
distributions of $\\{\widehat{\eta}_{ni}:i\in\mathcal{A}_{1}\\}$ and
$\\{\widehat{\eta}_{ni}^{(K)}:i\in\mathcal{A}_{2}\\}$,
$\\{\widehat{\eta}_{ni}:i\in\mathcal{A}_{1}\\}$ is independent of
$\\{\widehat{\eta}_{ni}^{(K)}:i\in\mathcal{A}_{2}\\}$, and
$\widehat{H}_{ij}^{(K)}=\kappa(\widehat{\eta}_{i},\widehat{\eta}_{j})+\kappa(\widehat{\eta}_{i}^{(K)},\widehat{\eta}_{j}^{(K)})-\kappa(\widehat{\eta}_{i},\widehat{\eta}_{j}^{(K)})-\kappa(\widehat{\eta}_{i}^{(K)},\widehat{\eta}_{j})$.

The proof of Theorem 1 is relegated to Supplementary Material. The difference
between (10) and the one in Theorem 1 is that the estimated residuals are
considered instead. Theorem 1 validates that instead of the true residuals, we
can use the estimated values, $\\{\widehat{\eta}_{ni}:i\in\mathcal{A}_{1}\\}$
and $\\{\widehat{\eta}_{ni}^{(K)}:i\in\mathcal{A}_{2}\\}$, for analysis. Thus,
$\widetilde{\mbox{MMD}}_{u}^{2}(\widetilde{\mathbb{P}},\widetilde{\mathbb{P}}^{(K)};\kappa)$
serves as the test statistic for the hypothesis testing problem (8).

To guarantee independence between
$\\{\widehat{\eta}_{ni}:i\in\mathcal{A}_{1}\\}$ and
$\\{\widehat{\eta}_{ni}^{(K)}:i\in\mathcal{A}_{2}\\}$, we train
$\widehat{\theta}_{RF}^{(B)}(\bm{x})$ and
$\widehat{\theta}_{RF}^{(B)}(\bm{x}^{(K)})$ by
$\\{(\bm{x}_{i},y_{i}):i\in\mathcal{A}_{3}\\}$ and
$\\{(\bm{x}_{i}^{(K)},y_{i}):i\in\mathcal{A}_{4}\\}$ separately, where
$\mathcal{A}_{3}$ and $\mathcal{A}_{4}$ are two subsamples with the same
sample size $m_{2}$, such that $\mathcal{A}_{1},\ldots,\mathcal{A}_{4}$ are
mutually disjoint. Moreover, to avoid correlation, we suggest that
$\mathcal{A}_{1},\ldots,\mathcal{A}_{4}$ are generated from
$mathcal{I}_{n}\setminus\mathcal{A}_{0}$. For practical guidance, we suggest
that $m_{2}\geq 400$; see Section 4 for details.

It is hard for simple kernels to distinguish two distributions with complex
structures. For example, a translation-invariant Gaussian kernel requires a
large sample to distinguish two distributions, since it cannot identify
“direction” information around each mode in a multivariate setup; see Figure 1
and its discussion of Liu et al. (2020) for details. To avoid the drawbacks of
traditional parametric kernels, Liu et al. (2020) proposed obtaining one using
deep neural networks:

$\widehat{\kappa}_{\omega}=\operatorname*{arg\,max}_{\kappa_{\omega}}\frac{\widetilde{\mbox{MMD}}_{u}^{2}(\widehat{\mathbb{P}},\widehat{\mathbb{P}}^{(K)};\kappa_{\omega})}{\widetilde{\sigma}_{1,\lambda}(\widehat{\mathbb{P}},\widehat{\mathbb{P}}^{(K)};\kappa_{\omega})},$

where
$\kappa_{\omega}=\\{(1-\eta)\kappa[\phi_{\omega}(x),\phi_{\omega}(y)]+\eta\\}q(x,y)$,
$\eta>0$ is predefined, $\phi_{\omega}(x)$ is a deep neural network with
parameter $\omega$ that extracts features, $\kappa(x,y)$ and $q(x,y)$ are
Gaussian kernels with lengthscales $\sigma_{\phi}$ and $\sigma_{q}$,
respectively, and
$\widetilde{\sigma}_{1,\lambda}^{2}(\widehat{\mathbb{P}},\widehat{\mathbb{P}}^{(K)};\kappa_{\omega})=4m_{1}^{-3}\sum_{i=1}^{m_{1}}(\sum_{j=1}^{m_{1}}H_{ij}^{(K)})^{2}-4m_{1}^{-4}(\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{1}}H_{ij}^{(K)})^{2}+\lambda$.

The limiting distributions of the U-statistic
$\widetilde{\mbox{MMD}}_{u}^{2}(\widehat{\mathbb{P}},\widehat{\mathbb{P}}^{(K)};\mathcal{H}_{\kappa})$
are established below for the null and alternative hypotheses.

######  Lemma 1.

Under the null hypothesis $H_{0K}^{\prime}:\mathbb{P}=\mathbb{P}^{(K)}$ and
regularity conditions in Supplementary Material, we have

$m_{1}\widetilde{\mbox{MMD}}_{u}^{2}(\widehat{\mathbb{P}},\widehat{\mathbb{P}}^{(K)};\kappa_{\omega})\to\sum_{i=1}^{\infty}\zeta_{i}(Z_{i}^{2}-2)$
(11)

in distribution, where $\zeta_{i}$ are the eigenvalues satisfying

$\int\kappa_{\omega}(\bm{x},\bm{z})\Psi_{i}(\bm{x})\mathbb{P}(\mathrm{d}\bm{x})=\zeta_{i}\Psi_{i}(\bm{z})$

for $i=1,2,\ldots,$, $\\{\Psi_{i}(\bm{x}):i=1,2,\ldots\\}$ are the
eigenfunctions, $Z_{i}\sim N(0,2)$, and $\mathcal{N}(\mu,\sigma^{2})$ is a
normal distribution with mean $\mu$ and variance $\sigma^{2}$.

Under the alternative hypothesis
$H_{aK}^{\prime}:\mathbb{P}\neq\mathbb{P}^{(K)}$, we have

$\sqrt{m_{1}}[\widehat{\mbox{MMD}}_{u}^{2}-MMD^{2}]\to\mathcal{N}(0,\sigma^{2}_{1})$

in the distribution, where
$\sigma^{2}_{1}=4[\mathbb{E}(H_{12}H_{13})-\mathbb{E}(H_{12})^{2}]$.

By Theorem 1, we can prove Lemma 1 in a manner similar to Theorem 12 of
Gretton et al. (2012), so its proof is omitted. Instead of deriving the
asymptotic distributions in Lemma 1, Liu et al. (2020) suggested permutation
for hypothesis testing. Algorithm 2 shows the routine of FS-D, and the
detailed algorithm for the MMD-D step is relegated to Supplementary Material.

Algorithm 2 Feature Selection by Deep Neural Network (FS-D)

Input:$\mathcal{D}_{n}$, $m_{2}$, $\\{I_{k}:k=1,\ldots,p\\}$, and $\alpha$.

Generates mutually disjoint sets $\mathcal{A}_{1},\ldots,\mathcal{A}_{4}$ from
$\mathcal{I}_{n}\setminus\mathcal{A}_{0}$.

Train $\widehat{\theta}_{RF}^{(B)}(\bm{x})$ by
$\\{(\bm{x}_{i},y_{i}):i\in\mathcal{A}_{3}\\}$.

Obtain $\\{\widehat{\eta}_{ni}:i\in\mathcal{A}_{1}\\}$.

for $k=1,\ldots,p$ do

if $K=1$ then

Train KRR by $\\{(\bm{x}_{i}^{(1)},y_{i}):i\in\mathcal{A}_{4}\\}$.

Obtain $\\{\widehat{\eta}_{ni}^{(K)}:i\in\mathcal{A}_{2}\\}$.

else

Train $\widehat{\theta}_{RF}^{(B)}(\bm{x}^{(K)})$ by
$\\{(\bm{x}_{i}^{(K)},y_{i}):i\in\mathcal{A}_{3}\\}$.

Obtain $\\{\widehat{\eta}_{ni}^{(K)}:i\in\mathcal{A}_{2}\\}$.

end if

Conduct MMD-D with significant level $\alpha$.

if $H_{0K}^{\prime}$ is rejected then

Continue.

else

Denote $\widehat{K}_{0}=K$.

Break the for loop.

end if

end for

Output: Selected feature set $\mathcal{P}$.

######  Remark 5.

The basic idea of the NFSRD is to apply nonparametric two-sample tests
sequentially based on ordered features. Thus, the NFSRD is also applicable to
other machine learning algorithms, such as XGBoost (Chen and Guestrin, 2016)
and the causal tree (Wager and Athey, 2018), as long as the features can be
sorted accordingly.

##  4 Simulation

In this section, we conduct Monte Carlo simulations to compare the performance
of the NFSRD with its alternatives under different model setups. As mentioned
in Section 1, there are two main approaches for feature selection by random
forests. Thus, we compare the NFSRD with a feature-importance-based algorithm
and a minimum-depth-based one.

We use 19 Xeon Cascade Lake (2.5 GHz) CPUs to train RFs in parallel, and an
NVIDIA Tesla T4 GPU is used for FS-D. Five-fold cross validation is conducted
to tune the model parameters for RFs; see Supplementary Material for more
computational details.

###  4.1 Independent Features

Table 1 shows the setups for six synthetic models. Model 1 is linear and is
widely used for feature selection (Tibshirani, 1996). Model 2 represents a
nonlinear model and only involves one useful feature. Both Model 1 and Model 2
are additive, and Model 3 is more complex and is non-additive. Besides,
interaction is involved in Model 3. The difference between Models 1–3 and
Models 4–6 is the distribution to generate features. Specifically, features
are generated independently from a uniform distribution over $[1,10]$ for
Models 1–3, and a _skewed_ beta distribution with shape parameter (2,4) is
used for Models 4–6. The distribution parameters are chosen such that the
features have approximately the same variance in different synthetic models,
and the regression parameters are selected such that the signal-to-noise ratio
(SNR) ranges from 1 to 3 approximately, where SNR is obtained by
$\\{\mathrm{var}{[f(\bm{X})]}/\mathrm{var}{(\epsilon)}\\}^{1/2}$, and
$\mathrm{var}{(X)}$ is the variance of a random variable $X$. For each model,
$\epsilon\sim\mathcal{N}(0,1)$. Furthermore, we consider two feature sizes,
including $p=200$ and $p=400$. Although the number of features is large, the
useful ones are limited. The sample size is $n=10\,000$. For each synthetic
model, we conduct $200$ independent Monte Carlo simulations.

Table 1: Setups for synthetic models. “MI” stands for the model index, and “Dist” represents the distribution to generate features. “Uniform” corresponds to $X_{k}\sim\mbox{U}(1,10)$, and “Beta” relates to $X_{k}\sim 14.5\mbox{Beta}(2,4)+1$, where $k=1,\ldots,p$. The error term is $\epsilon\sim\mathcal{N}(0,1)$.  MI | Dist | Model Setup | SNR
---|---|---|---
1 | Uniform | $Y=0.3X_{1}+0.3X_{2}+\epsilon$ | 1.1
2 | Uniform | $Y=3\sin(X_{1})+\epsilon$ | 2.1
3 | Uniform | $Y=5\sin(X_{1}/10)\sqrt{X_{2}}+\epsilon$ | 3.0
4 | Beta | $Y=0.3X_{1}+0.3X_{2}+\epsilon$ | 1.1
5 | Beta | $Y=3\sin(X_{1})+\epsilon$ | 2.1
6 | Beta | $Y=5\sin(X_{1}/10)\sqrt{X_{2}}+\epsilon$ | 2.8

We compare the following methods in terms of feature selection under the
significance level $\alpha=0.05$:

  1. 1.

Boruta (Kursa and Rudnicki, 2010, BRT). BRT calculates the feature importance
by a set of shadow features, and a “Z score” is used for feature selection;
see Supplementary Material for a brief description of BRT. BRT is implemented
using the R package `Boruta`.

  2. 2.

BRT-N. NFSRD is conducted based on the features sorted by BCFI, as obtained
from BRT.

  3. 3.

Minimal depth variable selection (Ishwaran et al., 2010, MVS). MVS selects
features using the minimum depth (Ishwaran et al., 2010); see Supplementary
Material for a brief description of MVS. This method is implemented using the
R package `randomForestSRC`.

  4. 4.

MVS-N. NFSRD is conducted based on the features sorted by the minimum depth
obtained from MVS.

Although `VSURF` (Genuer et al., 2015) also conducts feature selection by RFs,
it is not considered because of its heavy computational burden. For both BRT-N
and MVS-N, we consider two scenarios: $m_{0}=m_{1}=m_{2}=200$ and
$m_{0}=m_{1}=m_{2}=400$. For a fair comparison, $5m_{0}$ instances are used
for feature selection by BRT and MVS.

The four methods are compared in terms of the accuracy of identifying useful
features, $\mu_{c}$, and the average number of useless features included,
$n_{w}$, where
$\mu_{c}=(200K_{0})^{-1}\sum_{l=1}^{200}\sum_{k=1}^{K_{0}}\mathds{I}_{lk}$,
$n_{w}=200^{-1}\sum_{l=1}^{200}n_{lw}$, $\mathds{I}_{lk}=1$ is the $k$th
useful feature that is correctly identified in the $l$th Monte Carlo
simulation, and 0 otherwise, $n_{lw}$ is the number of useless features
identified in the $l$th Monte Carlo simulation. The simulation results are
summarized in Table 2. When the subsample sizes are small, $\mu_{c}$ is less
than 1 for Model 1, indicating that it is difficult for the two NFSRD
algorithms to identify some useful features even for a linear model when the
SNR is small. This is because $\widehat{\theta}_{RF}^{(B)}(\bm{x})-g(\bm{x})$
or $\widehat{\theta}_{RF}^{(B)}(\bm{x}^{(K)})-g(\bm{x}^{(K)})$ is not
negligible compared with the residual; thus, we cannot test (8) correctly.
However, as the SNR increases, both BRT-N and MVS-N identify more useful
features on average, even for small subsamples and more complex models. On the
one hand, even when $m_{0}$ is small, BRT and MVS, on the other hand, can
correctly select useful features with $5m_{0}$ instances for different models
because $5m_{0}$ instances are used for feature selection. As $m_{1}$ and
$m_{2}$ increase, useful features can be correctly identified by BRT-N and
MVS-N as well as BRT and MVS for different models. However, BRT-N and MVS-N
identify far less useless features compared with BRT and MVS, especially when
$m_{1}$ is large. For example, when $m_{1}=200$ and $p=400$, BRT and MVS
identify 2.17 and 34.88 useless features on average for a linear model,
respectively, but the average numbers are only 0.4 and 0.3 for BRT-N and
MVS-N. The same conclusion applies to the other setups. As $m_{1}$ and $m_{2}$
increase to 400, BRT-N and MVS-N still outperform their alternatives in terms
of avoiding useless features.

Table 2: Summary of feature selection the four methods, including BRT, BRT-N, MVS, and MVS-N, based on 200 Monte Carlo simulations, where the features are generated independently. “MI” stands for the model index. $m_{0}$, $m_{1}$ and $m_{2}$ are the sizes of subsamples for NFSRD algorithms. “BRT” represent Boruta, “MVS” stands for the minimum-depth variable selection algorithm, and “BRT-N” and “MVS-N” are NFSRD algorithms with features sorted by BRT and MVS, respectively. $\mu_{c}$ is proportion that useful features are correctly identified, and $n_{w}$ is the average length of useless features are included. The best result is highlighted by underline.  $p$ | MI |  | $m_{0}=m_{1}=m_{2}=200$ |  | $m_{0}=m_{1}=m_{2}=400$
---|---|---|---|---|---
BRT | BRT-N | MVS | MVS-N |  | BRT | BRT-N | MVS | MVS-N
200 | 1 | $\mu_{c}$ | 1.00 | 0.90 | 1.00 | 0.90 |  | 1.00 | 0.98 | 1.00 | 0.98
$n_{w}$ | 1.75 | 0.38 | 4.91 | 0.34 |  | 1.35 | 0.40 | 0.02 | 0.46
2 | $\mu_{c}$ | 1.00 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
$n_{w}$ | 2.16 | 0.60 | 5.17 | 0.59 |  | 1.42 | 0.54 | 0.12 | 0.49
3 | $\mu_{c}$ | 1.00 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
$n_{w}$ | 1.52 | 0.83 | 1.21 | 0.95 |  | 0.98 | 0.82 | 0.00 | 0.88
4 | $\mu_{c}$ | 1.00 | 0.78 | 1.00 | 0.81 |  | 1.00 | 0.90 | 1.00 | 0.88
$n_{w}$ | 2.21 | 0.26 | 6.83 | 0.32 |  | 1.43 | 0.32 | 0.04 | 0.43
| 5 | $\mu_{c}$ | 1.00 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
| $n_{w}$ | 1.98 | 0.66 | 4.42 | 0.71 |  | 1.30 | 0.54 | 0.08 | 0.40
| 6 | $\mu_{c}$ | 1.00 | 0.94 | 1.00 | 0.93 |  | 1.00 | 1.00 | 1.00 | 1.00
| $n_{w}$ | 1.57 | 0.78 | 1.61 | 0.75 |  | 1.13 | 0.66 | 0.00 | 0.74
400 | 1 | $\mu_{c}$ | 1.00 | 0.86 | 1.00 | 0.85 |  | 1.00 | 0.98 | 1.00 | 0.98
$n_{w}$ | 2.19 | 0.37 | 23.72 | 0.24 |  | 1.44 | 0.42 | 5.50 | 0.46
2 | $\mu_{c}$ | 1.00 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
$n_{w}$ | 2.35 | 0.64 | 19.60 | 0.56 |  | 1.34 | 0.58 | 5.43 | 0.44
3 | $\mu_{c}$ | 1.00 | 0.99 | 1.00 | 0.99 |  | 1.00 | 1.00 | 1.00 | 1.00
$n_{w}$ | 1.79 | 1.00 | 3.98 | 1.16 |  | 1.12 | 0.66 | 1.21 | 0.68
4 | $\mu_{c}$ | 1.00 | 0.90 | 1.00 | 0.90 |  | 1.00 | 0.99 | 1.00 | 0.99
$n_{w}$ | 2.17 | 0.40 | 34.88 | 0.30 |  | 1.48 | 0.52 | 7.03 | 0.61
| 5 | $\mu_{c}$ | 1.00 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
| $n_{w}$ | 2.30 | 0.73 | 28.79 | 0.47 |  | 1.38 | 0.47 | 4.75 | 0.58
| 6 | $\mu_{c}$ | 1.00 | 0.87 | 1.00 | 0.88 |  | 1.00 | 1.00 | 1.00 | 1.00
| $n_{w}$ | 1.58 | 0.61 | 17.16 | 0.66 |  | 0.99 | 0.81 | 1.28 | 0.76

We also compare the four methods in terms of computational efficiency, and
Figure 1 shows boxplots of the computation time based on the 200 Monte Carlo
simulations under different setups. When both $m_{0}$ and $p$ are small, the
computational efficiencies of the four methods are comparable, but, as shown
in Table 2, BRT and MVS select more useless features than BRT-N and MVS-N.
Thus, BRT and MVS are still less favored. As $m_{0}$ or $p$ increases, the
computational efficiency of BRT-N and MVS-N is generally much better,
especially when both $m_{0}$ and $p$ are large. For example, when $m_{0}=400$
and $p=400$, the computation time required by BRT-N and MVS-N is less than 50
seconds in general, but it is more than 125 seconds and approximately 75
seconds for BRT and MVS, respectively. In addition, it is far less efficient
to use BRT and MVS to select features for large datasets with numerous useless
features. However, BRT-N and MVS-N can be used to solve this problem. Compared
with MVS-N, BRT-N is slightly more computationally efficient in general, but
the difference between the two methods is minor.

Figure 1: Boxplots for the computation time of different methods based on 200
Monte Carlo simulations, where the features are independently generated. The
horizontal segments within the boxplots represent the median computation time.
“BRT” represent Boruta, “MVS” stands for the minimum-depth variable selection
algorithm, and “BRT-N” and “MVS-N” are NFSRD algorithms with features sorted
by BRT and MVS, respectively.

######  Remark 6.

It is common to have correlated features in practice, and we conduct a
simulation study for this case as well. We still consider the same setups in
Table 1, but the features are generated differently. First, generate
$X_{1}^{\prime},\ldots,X_{p}^{\prime}$ independently by a uniform or skewed
beta distribution, as shown in Table 1. Then, let $X_{1}=X_{1}^{\prime}$ and
$X_{k}=0.7X_{k}^{\prime}+0.3X_{k-1}^{\prime}$ for $k=2,\ldots,p$. The features
consist of $X_{1},\ldots,X_{p}$, and we still use the six models in Table 1 to
generate the responses of interest. The simulation results are similar to the
aforementioned outcomes, and we relegate them to Supplementary Material.

Ishwaran et al. (2010) compared the MVS with some commonly used feature
selection algorithms, including the adaptive lasso (Zou, 2006) and the
$l_{1}$-regularized regression model (Park and Hastie, 2007). Their simulation
results demonstrated that the MVS outperforms those two in terms of the false
discovery rate and the false nondiscovery rate, regardless of the correlation
among features. Since the MVS-N is generally more preferable than MVS, we do
not consider those two feature selection algorithms in the simulation study.

##  5 Application

The superconductivity dataset Hamidieh (2018) is used to test the performance
of the NFSRD. Superconducting materials have wide applications in practice,
such as magnetic resonance imaging systems in hospitals and superconducting
coils in the Large Hadron Collider at CERN. The accurate prediction of the
superconducting critical temperature is important because the corresponding
superconductor can only conduct current without resistance at or below this
temperature. There are 21 263 instances in the superconductivity dataset with
81 features extracted for each instance, and the goal is select useful
features for the critical temperature; see Hamidieh (2018) for details about
the technical information of the features.

For BRT-N, we consider two significance levels: $\alpha=0.05$ and
$\alpha=0.01$, and the size of $\mathcal{A}_{0},\ldots,\mathcal{A}_{4}$ is
$1\,000$. Table 3 lists the selected useful features. As the significance
level decreases from 0.05 to 0.01, three more features are selected. In
addition, all the selected features are among the top 20 most important
features by XGBoost (Chen and Guestrin, 2016); see Table 5 of Hamidieh (2018)
for details. In contrast, BRT blindly identifies all the features to be
useful, even for $\alpha=0.05$. Although BRT works reasonably for synthetic
models in Section 4, it fails to identify useful features for
superconductivity data. A similar conclusion holds for the two minimum-depth-
based methods, and, thus, we omit them.

Table 3: The selected features by BRT-N under different significance levels for the superconductivity data. Selected features are indicated by “$\checkmark$”, and “Feature” shows the corresponding names. “$\alpha$” is the significant level. Feature | $\alpha$
---|---
$0.01$ | $0.05$
`range_ThermalConductivity` | $\checkmark$ | $\checkmark$
`wtd_std_ThermalConductivity` | $\checkmark$ | $\checkmark$
`range_atomic_radius` | $\checkmark$ | $\checkmark$
`wtd_entropy_atomic_mass` | $\checkmark$ | $\checkmark$
`wtd_entropy_Valence` | $\checkmark$ | $\checkmark$
`wtd_mean_Valence` | $\checkmark$ | $\checkmark$
`wtd_mean_ThermalConductivity` | $\checkmark$ | $\checkmark$
`std_ThermalConductivity` | $\checkmark$ | $\checkmark$
`wtd_gmean_Valence` | $\checkmark$ |
`wtd_gmean_ThermalConductivity` | $\checkmark$ |
`wtd_std_ElectronAffinity` | $\checkmark$ |

##  6 Concluding Remarks

In this study, we propose the NFSRD to identify useful features. Feature
selection is conducted by nonparametric two-sample tests using deep neural
networks, and the theoretical properties are also investigated. Experiments
show that the NFSRD outperforms its alternatives in terms of identifying
useful features, avoiding useless ones and the computation efficiency.

Although the NFSRD is proposed using the standard RFs (Breiman, 2001), the
same idea can be easily adapted to other machine learning methods. In
addition, other nonparametric tests can also be used; these include other
nonparametric kernel-based methods (Scholkopf and Smola, 2018) and traditional
statistical tests. However, we should pay attention to some widely used tests,
such as the $t$-test and the Kolmogorov–Smirnov test, because they may suffer
from model misspecification.

##  7 Acknowledgements

The authors thank the editor, associate editor, and two anonymous reviewers
for their detailed and constructive comments.

This work was supported by the National Natural Science Foundation of China
(NSFC) (No. 12001109, No. 92046021, No. 11901487, No. 72033002), Science and
Technology Commission of Shanghai Municipality grant (No. 20dz1200600) and
Fundamental Research Funds for the Central Universities(No. 20720191056).

## References

  * (1)
  * Altmann et al. (2010) Altmann, A., Toloşi, L., Sander, O. and Lengauer, T. (2010).  Permutation importance: a corrected feature importance measure, Bioinformatics 26(10): 1340–1347.
  * Benci et al. (2016) Benci, J. L., Xu, B., Qiu, Y., Wu, T. J., Dada, H., Twyman-Saint Victor, C., Cucolo, L., Lee, D. S., Pauken, K. E., Huang, A. C. et al. (2016).  Tumor interferon signaling regulates a multigenic resistance program to immune checkpoint blockade, Cell 167(6): 1540–1554.
  * Breiman (2001) Breiman, L. (2001).  Random forests, Machine Learning 45(1): 5–32.
  * Caruana and Niculescu-Mizil (2006) Caruana, R. and Niculescu-Mizil, A. (2006).  An empirical comparison of supervised learning algorithms, Proceedings of the 23rd International Conference on Machine Learning, pp. 161–168.
  * Chen and Guestrin (2016) Chen, T. and Guestrin, C. (2016).  Xgboost: A scalable tree boosting system, Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 785–794.
  * Criminisi et al. (2012) Criminisi, A., Shotton, J. and Konukoglu, E. (2012).  Decision forests: A unified framework for classification, regression, density estimation, manifold learning and semi-supervised learning, Foundations and Trends® in Computer Graphics and Vision 7(2–3): 81–227.
  * Fan and Lv (2008) Fan, J. and Lv, J. (2008).  Sure independence screening for ultrahigh dimensional feature space, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70(5): 849–911.
  * Fernández-Delgado et al. (2014) Fernández-Delgado, M., Cernadas, E., Barro, S. and Amorim, D. (2014).  Do we need hundreds of classifiers to solve real world classification problems?, Journal of Machine Learning Research 15(1): 3133–3181.
  * Fukumizu et al. (2007) Fukumizu, K., Gretton, A., Sun, X. and Schölkopf, B. (2007).  Kernel measures of conditional dependence., Advances in Neural Information Processing Systems, Vol. 20, pp. 489–496.
  * Genuer et al. (2010) Genuer, R., Poggi, J.-M. and Tuleau-Malot, C. (2010).  Variable selection using random forests, Pattern Recognition Letters 31(14): 2225–2236.
  * Genuer et al. (2015) Genuer, R., Poggi, J.-M. and Tuleau-Malot, C. (2015).  VSURF: an R package for variable selection using random forests, The R journal 7(2): 19–33.
  * Goel et al. (2017) Goel, E., Abhilasha, E., Goel, E. and Abhilasha, E. (2017).  Random forest: A review, International Journal of Advanced Research in Computer Science and Software Engineering 7(1): 251–257.
  * Gretton et al. (2012) Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B. and Smola, A. (2012).  A kernel two-sample test, The Journal of Machine Learning Research 13(1): 723–773.
  * Guyon and Elisseeff (2003) Guyon, I. and Elisseeff, A. (2003).  An introduction to variable and feature selection, Journal of Machine Learning Research 3(Mar): 1157–1182.
  * Haghiri et al. (2018) Haghiri, S., Garreau, D. and von Luxburg, U. (2018).  Comparison-based random forests, Proceedings of the 35th International Conference on Machine Learning, pp. 1871–1880.
  * Hamidieh (2018) Hamidieh, K. (2018).  A data-driven statistical model for predicting the critical temperature of a superconductor, Computational Materials Science 154: 346–354.
  * Hastie et al. (2009) Hastie, T., Tibshirani, R. and Friedman, J. (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction, 2nd edn, Springer, New York.
  * Ishwaran et al. (2008) Ishwaran, H., Kogalur, U. B., Blackstone, E. H. and Lauer, M. S. (2008).  Random survival forests, Annals of Applied Statistics 2(3): 841 – 860.
  * Ishwaran et al. (2010) Ishwaran, H., Kogalur, U. B., Gorodeski, E. Z., Minn, A. J. and Lauer, M. S. (2010).  High-dimensional variable selection for survival data, Journal of the American Statistical Association 105(489): 205–217.
  * Kursa and Rudnicki (2010) Kursa, M. B. and Rudnicki, W. R. (2010).  Feature selection with the Boruta package, Journal of Statistical Software 36(1): 1–13.
  * Lakshminarayanan et al. (2014) Lakshminarayanan, B., Roy, D. M. and Teh, Y. W. (2014).  Mondrian forests: Efficient online random forests, Advances in Neural Information Processing Systems, pp. 3140–3148.
  * Li and Martin (2017) Li, A. H. and Martin, A. (2017).  Forest-type regression with general losses and robust forest, Proceedings of the 34th International Conference on Machine Learning, pp. 2091–2100.
  * Li et al. (2019) Li, X., Wang, R., Basu, S., Kumbier, K. and Yu, B. (2019).  A debiased MDI feature importance measure for random forests, Advances in Neural Information Processing Systems, Vol. 32, pp. 1–22.
  * Liu et al. (2020) Liu, F., Xu, W., Lu, J., Zhang, G., Gretton, A. and Sutherland, D. J. (2020).  Learning deep kernels for non-parametric two-sample tests, Proceeding of the 37th International Conference on Machine Learning, pp. 6316–6326.
  * Louppe et al. (2013) Louppe, G., Wehenkel, L., Sutera, A. and Geurts, P. (2013).  Understanding variable importances in forests of randomized trees, Advances in Neural Information Processing Systems, pp. 431–439.
  * Meinshausen (2006) Meinshausen, N. (2006).  Quantile regression forests, Journal of Machine Learning Research 7(Jun): 983–999.
  * Mentch and Hooker (2016) Mentch, L. and Hooker, G. (2016).  Quantifying uncertainty in random forests via confidence intervals and hypothesis tests, Journal of Machine Learning Research 17(1): 841–881.
  * Park and Hastie (2007) Park, M. Y. and Hastie, T. (2007).  L1-regularization path algorithm for generalized linear models, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 69(4): 659–677.
  * Payet and Todorovic (2010) Payet, N. and Todorovic, S. (2010).  (RF)^2 – random forest random field, Advances in Neural Information Processing Systems, pp. 1885–1893.
  * Sandri and Zuccolotto (2008) Sandri, M. and Zuccolotto, P. (2008).  A bias correction algorithm for the Gini variable importance measure in classification trees, Journal of Computational and Graphical Statistics 17(3): 611–628.
  * Scholkopf and Smola (2018) Scholkopf, B. and Smola, A. J. (2018).  Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond, MIT Press, London.
  * Scornet et al. (2015) Scornet, E., Biau, G. and Vert, J.-P. (2015).  Consistency of random forests, Annals of Statistics 43(4): 1716–1741.
  * Shi and Horvath (2006) Shi, T. and Horvath, S. (2006).  Unsupervised learning with random forest predictors, Journal of Computational and Graphical Statistics 15(1): 118–138.
  * Siblini et al. (2018) Siblini, W., Kuntz, P. and Meyer, F. (2018).  CRAFTML, an efficient clustering-based random forest for extreme multi-label learning, Proceedings of the 35th International Conference on Machine Learning, pp. 4671–4680.
  * Strobl et al. (2008) Strobl, C., Boulesteix, A.-L., Kneib, T., Augustin, T. and Zeileis, A. (2008).  Conditional variable importance for random forests, BMC Bioinformatics 9(1): 307.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society: Series B (Methodological) 58(1): 267–288.
  * Twyman-Saint Victor et al. (2015) Twyman-Saint Victor, C., Rech, A. J., Maity, A., Rengan, R., Pauken, K. E., Stelekati, E., Benci, J. L., Xu, B., Dada, H., Odorizzi, P. M. et al. (2015).  Radiation and dual checkpoint blockade activate non-redundant immune mechanisms in cancer, Nature 520(7547): 373–377.
  * Wager and Athey (2018) Wager, S. and Athey, S. (2018).  Estimation and inference of heterogeneous treatment effects using random forests, Journal of the American Statistical Association 113(523): 1228–1242.
  * White and Liu (1994) White, A. P. and Liu, W. Z. (1994).  Bias in information-based measures in decision tree induction, Machine Learning 15(3): 321–329.
  * Xu and Jelinek (2005) Xu, P. and Jelinek, F. (2005).  Using random forests in the structured language model, Advances in Neural Information Processing Systems, pp. 1545–1552.
  * Zou (2006) Zou, H. (2006).  The adaptive lasso and its oracle properties, Journal of the American Statistical Association 101(476): 1418–1429.


provide a gravitational-analogue for particle creation by considering an
initial half-filling condition that first populates the lower band of the
lattice models. Then, a period of de Sitter expansion between two asymptotic
flat vacua is simulated by the real-time evolution of the Fermi gas loaded in
the Raman lattice with a time-dependent detuning (83). During the expansion,
not only the state of the Fermi gas will change, but also the Hamiltonian
itself, and consequently the band structure of the system. Hence, if a
fermionic atom is initially in a specific Bloch state, after the expansion, it
will generally be left in a superposition of Bloch states corresponding to the
two different bands, which is the analogue of Eq. (27) provided a particle-
hole transformation is applied. The number of particles produced for a certain
quasi-momentum, given by the squared modulus of the Bogoliubov coefficient
$|\beta_{\rm k}(t+\eta_{0})|^{2}$ is then understood, according to Eq. (42),
as the probability of measuring the evolved state in the upper band of the
final Hamiltonian for a certain value of the quasi-momentum. Since we can
measure the population of each band for each quasi-momentum with the band-
mapping technique Tarruell  _et al._ (2012), this Bogoliubov coefficient can
thus be accessed experimentally. Going further, combining the band-mapping
technique with spin-resolved measurements Alba  _et al._ (2011), one could
even characterize the band topology of the expanding lattice field theory
Cooper  _et al._ (2019).

##  V Conclusions and outlook

In this work, we have developed the theory of particle production for a
fermionic Dirac field in a (1+1)-dimensional FRW spacetime, both from the
perspective of the usual continuous theory and two types of discretizations on
the lattice. We have shown that the phenomenon of fermion production for a de
Sitter expansion admits an exact solution in terms of a pair of decoupled
Bessel differential equations. To avoid problems with the interpretation of
the instantaneous vacua, we added an adiabatic switching that connects the de
Sitter expansion to a pair of asymptotic Minkowski spacetimes, both of which
have a well-defined notion of vacuum. We have numerically shown that the extra
switching periods do not modify the particle production, which can still be
described in terms of the analytically-solvable Bessel equation.

To pave the way for a quantum simulation of this phenomenon, we have
considered two possible lattice discretizations of the Hamiltonian field
theory associated to Dirac fermions in a curved spacetime: a naïve- and a
Wilson-fermion discretization, which allow us to discuss universes with
spatial boundaries. Focusing on the bulk of the lattice, we have shown that
the naïve-fermion discretization reproduces the spectrum of particle
production predicted by the continuum QFT for small momenta. However, as one
approaches the edge of the Brillouin zone, the phenomenon of fermion doubling
leads to a mirror-image of the spectrum, resulting in a doubled number of the
produced fermions with respect to the continuum prediction. Turning into the
Wilson-fermion discretization, which sends the fermion doubler to the ultra-
violet cutoff of the QFT, we have shown that the spectrum and total number of
particles produced at the bulk matches nicely the continuum expressions.

On the other hand, the asymptotic vacua of the Wilson-fermion lattice field
theories can actually correspond to a couple of SPT groundstates characterized
by a non-zero topological invariant, which has a boundary correspondence in
terms of the groundstate degeneracy by the appearance of topological zero-
energy modes localized to the spatial boundaries of the FRW universe. This has
allowed us to study the role of topology in the phenomenon of particle
production and, in particular, to show that there can be gravitational
particle production of zero-energy fermions localized to the edges of the
spacetime. A similar phenomenon would also occur in situations where the mass
of the Dirac fermion has some solitonic profile that changes its sign, leading
to particle production in the form of domain-wall fermions. By numerically
solving a set of coupled differential equations for the case of open boundary
conditions, we have calculated the spectrum of fermion production, which is
now a function of the energy of the produced particles as momentum is no
longer a good quantum number in the absence of translational symmetry. We have
shown that this spectrum allows one to identify clearly the production of
fermions at the boundary of the expanding FRW universe, as one finds
production of particles for energies below the mass gap.

To conclude, we have shown that current experiments of ultra-cold alkaline-
earth Fermi gases in Raman optical lattices would be an ideal platform for the
quantum simulation of this curved quantum field theories, provided that one
exploits the simplifications that arise when working with conformal time. In
connection to those experiments, our model requires to simulate the expanding
spacetime by encoding its effect in a time-dependent mass that depends on the
scale factor of the expansion, which corresponds in the experiment to changing
the Raman-beam detuning as a function of time. This would allow for an
experimental simulation of fermion production in a FRW spacetime, including
its interplay with topology and the boundaries of such an effective universe.

###### Acknowledgements.

We thank S. Hands and C. Sabín for useful discussions. We acknowledge support
from PGC2018-099169- BI00 (MCIU/AEI/FEDER, UE), and PID2021- 127726NB-I00
(MCIU/AEI/FEDER, UE), from the Grant IFT Centro de Excelencia Severo Ochoa
CEX2020-001007-S, funded by MCIN/AEI/10.13039/501100011033, from the grant
QUITEMAD+ S2013/ICE-2801, and from the CSIC Research Platform on Quantum
Technologies PTI-001.

##  Appendix A Analytical solution for a de Sitter expansion

In this Appendix, we present the details of the derivation of Eq. (43). From
the single-particle Dirac Hamiltonian (38), we have

$\displaystyle\mathrm{i}\dot{u}_{1}=\frac{m}{H\eta}u_{1}-\mathrm{i}\mathrm{k}u_{2},$
(84a)
$\displaystyle\mathrm{i}\dot{u}_{2}=\mathrm{i}\mathrm{k}u_{1}-\frac{m}{H\eta}u_{2},$
(84b)

where we omit the subscript $k$ in the mode functions to ease the notation,
and use dots for the derivatives with respect to conformal time. Taking
derivatives on the first equation again

$\mathrm{i}\ddot{u}_{1}=-\frac{m}{H\eta^{2}}u_{1}+\frac{m}{H\eta}\dot{u}_{1}-\mathrm{i}\mathrm{k}\dot{u}_{2},$
(85)

such that one can now insert Eq.(84b) in (85) to get

$\mathrm{i}\ddot{u}_{1}=-\frac{m}{H\eta^{2}}u_{1}+\frac{m}{H\eta}\dot{u}_{1}-\mathrm{i}\mathrm{k}^{2}u_{1}+\frac{mk}{H\eta}u_{2}.$
(86)

We can find $u_{2}$ as a function of $u_{1}$ from Eq. (84a)

$u_{2}=\frac{1}{\mathrm{k}}\left(-\dot{u}_{1}-\mathrm{i}\frac{m}{H\eta}u_{1}\right),$
(87)

which, after substitution in Eq. (86), leads to

$\ddot{u}_{1}+\left(\mathrm{k}^{2}+\frac{m^{2}}{H^{2}\eta^{2}}-\mathrm{i}\frac{m}{H\eta^{2}}\right)u_{1}=0.$
(88)

Repeating this with $u_{2}$ we arrive at an analogous equation

$\ddot{u}_{2}+\left(\mathrm{k}^{2}+\frac{m^{2}}{H^{2}\eta^{2}}+\mathrm{i}\frac{m}{H\eta^{2}}\right)u_{2}=0.$
(89)

We can solve these equations by making the replacement
$u_{1,2}(\eta)=\sqrt{\eta}f_{1,2}(\eta)$. The functions $f_{(1,2)}$ satisfy
the ODEs

$\displaystyle\eta^{2}\ddot{f}_{1}+\eta\dot{f}_{1}+\left(\mathrm{k}^{2}\eta^{2}-\left(\textstyle\frac{1}{2}+\mathrm{i}\frac{m}{H}\right)^{2}\right)f_{1}=0,$
(90)
$\displaystyle\eta^{2}\ddot{f}_{2}+\eta\dot{f}_{1}+\left(\mathrm{k}^{2}\eta^{2}-\left(\textstyle\frac{1}{2}-\mathrm{i}\frac{m}{H}\right)^{2}\right)f_{2}=0.$
(91)

Finally, after performing the change of variable $\eta\rightarrow
z=\mathrm{k}\eta$, we arrive at

$\displaystyle
z^{2}f^{\prime\prime}_{1}+zf^{\prime}_{1}+\left(z^{2}-\left(\textstyle\frac{1}{2}+\mathrm{i}\frac{m}{H}\right)^{2}\right)f_{1}=0,$
(92) $\displaystyle
z^{2}f^{\prime\prime}_{2}+zf^{\prime}_{1}+\left(z^{2}-\left(\textstyle\frac{1}{2}-\mathrm{i}\frac{m}{H}\right)^{2}\right)f_{2}=0,$
(93)

where the prime denotes differentiation with respect to $z$. We have finally
arrived at a pair of decoupled Bessel equations for each of the spinor
components. Defining $\nu_{+}=\textstyle\frac{1}{2}+\mathrm{i}\frac{m}{H}$ for
(92) and $\nu_{-}=\textstyle\frac{1}{2}-\mathrm{i}\frac{m}{H}$ for (93), a
well-known basis of solutions is that of the so-called Hankel functions
Gradshteyn and Ryzhik (2007)

$\displaystyle
f_{1}=C_{(1,1)}H^{(1)}_{\nu_{+}}(\mathrm{k}\eta)+C_{(1,2)}H^{(2)}_{\nu_{+}}(\mathrm{k}\eta),$
(94) $\displaystyle
f_{2}=C_{(2,1)}H^{(1)}_{\nu_{-}}(\mathrm{k}\eta)+C_{(2,2)}H^{(2)}_{\nu_{-}}(\mathrm{k}\eta).$
(95)

We choose this basis instead of, for example, the more commonly used Bessel
functions of the first- $J_{\nu}(z)$ and second-kind $Y_{\nu}(z)$, because of
their asymptotic behaviour. In particular, the Hankel functions behave
asymptotically as plane waves for $z\rightarrow\infty$, which will be useful
for us when connecting to the vacuum of the remote past. As a consistency
check, let us note that our solution can be shown to be equivalent to the one
found by V. M. Villalba Villalba (1995), which was expressed in terms the so-
called cylinder functions. Although not apparent at first sight, this
equivalence can be proven by using the properties of the Hankel functions
$H^{(1)}_{-\nu}(z)={\rm e}^{{\rm i}\nu\pi}H^{(1)}_{\nu}(z)$,
$H^{(2)}_{-\nu}(z)={\rm e}^{-{\rm i}\nu\pi}H^{(2)}_{\nu}(z)$, and their
relation to the aforementioned cylinder functions
$Z_{\nu}(z)=aH^{(1)}_{\nu}(z)+bH^{(2)}_{\nu}(z)$,
$Z_{\nu-1}(z)=cH^{(1)}_{\nu-1}(z)+dH^{(2)}_{\nu-1}(z)$, where $a,b,c,d$ are
related to the constants of integration $C_{(i,j)}$ in (94)-(95).

##  Appendix B Particle production for infinitely-long expansions

In this Appendix, we present the details of the derivation of the Bogoliubov
parameter in Eq. (45). Let us begin from the limiting form of the
instantaneous eigenstates $\mathsf{v}_{k}^{-}(\eta_{\rm out}\rightarrow
0^{-})=(0,\,1)^{\rm t}$, and use the asymptotic solution of the Dirac spinor
in Eq.(44). We have that $|\beta_{\mathrm{k}}(\eta)|^{2}=|u_{k}^{\rm
t}(\eta)\cdot\mathsf{v}_{k}^{-}(\eta)|^{2}=|C_{\mathrm{k}}\sqrt{\eta}{\rm
e}^{\mathrm{i}\pi\nu_{-}}H_{\nu_{-}}^{(1)}(k\eta)|^{2}$, which is only valid
when $m\neq 0$ and for $\eta\approx 0^{-}$, since we used the limit
$\eta\rightarrow 0^{-}$. Let us recall that
$C_{\mathrm{k}}=\frac{\sqrt{\pi\mathrm{k}}}{2}{\rm
e}^{-\mathrm{i}(\mathrm{k}\eta_{\rm in}+\frac{\pi\nu_{-}}{2}-\frac{\pi}{4})}$,
and that the expansion of the Hankel function is function is
$H_{\nu}^{(1)}(z)\to-\frac{1}{\pi}\Gamma\left(\nu\right)\left(\textstyle\frac{1}{2}z\right)^{-\nu}$
as $z\to 0$. Using this, and substituting $\nu_{-}=\textstyle\frac{1}{2}-{\rm
i}\frac{m}{H}$, we find

$\displaystyle|\beta_{k}(\eta)|^{2}=\frac{1}{2\pi}e^{-\frac{\pi
m}{H}}\left|\Gamma\left(\textstyle\frac{1}{2}-\mathrm{i}\frac{m}{H}\right)\right|^{2}.$
(96)

Using the property of the gamma function
$\left|\Gamma\left(\textstyle\frac{1}{2}+\mathrm{i}b\right)\right|^{2}=\frac{\pi}{\text{cosh}(\pi
b)}$, we finally arrive at

$\displaystyle|\beta_{k}(\eta)|^{2}$ $\displaystyle=\frac{e^{\frac{-\pi
m}{H}}}{2\pi}\frac{\pi}{\text{cosh}\left(\frac{\pi
m}{H}\right)}=\frac{1}{1+{\rm e}^{\frac{2\pi m}{H}}}.$ (97)

which is the result that we stated in Eq. (45). Our derivation differs from
that found in Villalba (1995). There, no restriction on the mass was assumed,
which would imply particle production in the massless conformal limit $m=0$.
Additionally, our derivation also differs from that of Haro and Elizalde
(2008), which was only valid for $m/H\gg 0$. We note that, for the massless
case $m=0$, the previous asymptotic eigenvector is no longer valid, and must
be substituted by $\mathsf{v}_{\mathrm{k}}^{-}(\eta\rightarrow
0^{-})=\left(\mathrm{i},\,1\right)^{\rm t}/\sqrt{2}$. The corresponding
Bogoliubov parameter vanishes, as it should be, since in the massless case the
fermionic field is conformally coupled, and hence there is no particle
production.

## References

  * Peskin and Schroeder (1995) M. E. Peskin and D. V. Schroeder,  _An Introduction to quantum field theory_ (Addison-Wesley, Reading, USA, 1995).
  * Fradkin (2013) E. Fradkin,  _Field Theories of Condensed Matter Physics_, 2nd ed. (Cambridge University Press, 2013).
  * Sachdev (2011) S. Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, 2011).
  * Castro Neto  _et al._ (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Armitage  _et al._ (2018) N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Anderson (1972) P. W. Anderson, Science 177, 393 (1972).
  * Unruh (1981) W. G. Unruh, Phys. Rev. Lett. 46, 1351 (1981).
  * Carroll (2019) S. M. Carroll,  _Spacetime and Geometry: An Introduction to General Relativity_ (Cambridge University Press, 2019).
  * Barceló  _et al._ (2005) C. Barceló, S. Liberati, and M. Visser, Living Reviews in Relativity 8, 12 (2005).
  * Jacquet  _et al._ (2020) M. J. Jacquet, S. Weinfurtner, and F. König, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 378, 20190239 (2020).
  * Parker and Toms (2009) L. E. Parker and D. Toms,  _Quantum Field Theory in Curved Spacetime: Quantized Field and Gravity_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2009).
  * Hawking (1974) S. W. Hawking, Nature 248, 30 (1974).
  * Guth (1981) A. H. Guth, Phys. Rev. D 23, 347 (1981).
  * Belgiorno  _et al._ (2011) F. Belgiorno, S. L. Cacciatori, G. Ortenzi, L. Rizzi, V. Gorini, and D. Faccio, Phys. Rev. D 83, 024015 (2011).
  * Philbin  _et al._ (2008) T. G. Philbin, C. Kuklewicz, S. Robertson, S. Hill, F. König, and U. Leonhardt, Science 319, 1367 (2008).
  * Belgiorno  _et al._ (2010) F. Belgiorno, S. L. Cacciatori, M. Clerici, V. Gorini, G. Ortenzi, L. Rizzi, E. Rubino, V. G. Sala, and D. Faccio, Phys. Rev. Lett. 105, 203901 (2010).
  * Drori  _et al._ (2019) J. Drori, Y. Rosenberg, D. Bermudez, Y. Silberberg, and U. Leonhardt, Phys. Rev. Lett. 122, 010404 (2019).
  * Garay  _et al._ (2000) L. J. Garay, J. R. Anglin, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 85, 4643 (2000).
  * Lahav  _et al._ (2010) O. Lahav, A. Itah, A. Blumkin, C. Gordon, S. Rinott, A. Zayats, and J. Steinhauer, Phys. Rev. Lett. 105, 240401 (2010).
  * Steinhauer (2016) J. Steinhauer, Nature Physics 12, 959 (2016).
  * Kolobov  _et al._ (2021) V. I. Kolobov, K. Golubkov, J. R. Muñoz de Nova, and J. Steinhauer, Nature Physics 17, 362 (2021).
  * Unruh (1976) W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * Hu  _et al._ (2019) J. Hu, L. Feng, Z. Zhang, and C. Chin, Nature Physics 15, 785 (2019).
  * Parker (1968) L. Parker, Phys. Rev. Lett. 21, 562 (1968).
  * Parker (1969) L. Parker, Phys. Rev. 183, 1057 (1969).
  * Barceló  _et al._ (2003) C. Barceló, S. Liberati, and M. Visser, Phys. Rev. A 68, 053613 (2003).
  * Fedichev and Fischer (2004) P. O. Fedichev and U. R. Fischer, Phys. Rev. A 69, 033602 (2004).
  * Fey  _et al._ (2018) C. Fey, T. Schaetz, and R. Schützhold, Phys. Rev. A 98, 033407 (2018).
  * Eckel  _et al._ (2018) S. Eckel, A. Kumar, T. Jacobson, I. B. Spielman, and G. K. Campbell, Phys. Rev. X 8, 021021 (2018).
  * Wittemer  _et al._ (2019) M. Wittemer, F. Hakelberg, P. Kiefer, J.-P. Schröder, C. Fey, R. Schützhold, U. Warring, and T. Schaetz, Phys. Rev. Lett. 123, 180502 (2019).
  * Viermann  _et al._ (2022) C. Viermann, M. Sparn, N. Liebster, M. Hans, E. Kath, Á. Parra-López, M. Tolosa-Simeón, N. Sánchez-Kuntz, T. Haas, H. Strobel, S. Floerchinger, and M. K. Oberthaler, Nature 611, 260 (2022).
  * Tolosa-Simeón  _et al._ (2022) M. Tolosa-Simeón, A. Parra-López, N. Sánchez-Kuntz, T. Haas, C. Viermann, M. Sparn, N. Liebster, M. Hans, E. Kath, H. Strobel, M. K. Oberthaler, and S. Floerchinger, Phys. Rev. A 106, 033313 (2022).
  * Sanchez-Kuntz  _et al._ (2022) N. Sanchez-Kuntz, A. Parra-López, M. Tolosa-Simeón, T. Haas, and S. Floerchinger, Physical Review D 105 (2022), 10.1103/physrevd.105.105020.
  * Bloch  _et al._ (2008a) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008a).
  * Parker (1971) L. Parker, Phys. Rev. D 3, 346 (1971).
  * Kaluza (2018) T. Kaluza, International Journal of Modern Physics D 27, 1870001 (2018).
  * Klein (1926) O. Klein, Zeitschrift für Physik 37, 895 (1926).
  * Khanna  _et al._ (2014) F. Khanna, A. Malbouisson, J. Malbouisson, and A. Santana, Physics Reports 539, 135 (2014), quantum field theory on toroidal topology: Algebraic structure and applications.
  * Wald (2001) R. M. Wald, Living Reviews in Relativity 4, 6 (2001).
  * Casimir (1948) H. B. G. Casimir, Proc. Akad. Wet. Amsterdam 51, 793 (1948).
  * Moore (1970) G. T. Moore, Journal of Mathematical Physics 11, 2679 (1970).
  * Witten (1989) E. Witten, Communications in Mathematical Physics 121, 351 (1989).
  * Elitzur  _et al._ (1989) S. Elitzur, G. Moore, A. Schwimmer, and N. Seiberg, Nuclear Physics B 326, 108 (1989).
  * Balachandran  _et al._ (1995) A. Balachandran, L. Chándal, and E. Ercolessi, International Journal of Modern Physics A 10, 1969 (1995).
  * Golterman  _et al._ (1993) M. F. Golterman, K. Jansen, and D. B. Kaplan, Physics Letters B 301, 219 (1993).
  * Callan and Harvey (1985) C. Callan and J. Harvey, Nuclear Physics B 250, 427 (1985).
  * Kaplan (2009) D. B. Kaplan, in  _Les Houches Summer School: Session 93: Modern perspectives in lattice QCD: Quantum field theory and high performance computing_ (2009) pp. 223–272, arXiv:0912.2560 [hep-lat] .
  * Kaplan (1992) D. B. Kaplan, Physics Letters B 288, 342 (1992).
  * Feynman (1982) R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982).
  * Gross and Neveu (1974) D. J. Gross and A. Neveu, Phys. Rev. D 10, 3235 (1974).
  * Weinberg (1995) S. Weinberg,  _The Quantum Theory of Fields_, Vol. 1 (Cambridge University Press, 1995).
  * Freedman and Van Proeyen (2012) D. Z. Freedman and A. Van Proeyen,  _Supergravity_ (Cambridge University Press, 2012).
  * Kibble (1961) T. W. B. Kibble, Journal of Mathematical Physics 2, 212 (1961).
  * Shapiro (2016) I. L. Shapiro, (2016), arXiv:1611.02263 [gr-qc] .
  * Birrell and Davies (1984) N. D. Birrell and P. C. W. Davies,  _Quantum Fields in Curved Space_, Cambridge Monographs on Mathematical Physics (Cambridge Univ. Press, Cambridge, UK, 1984).
  * Mukhanov and Winitzki (2007) V. Mukhanov and S. Winitzki,  _Introduction to quantum effects in gravity_ (Cambridge University Press, 2007).
  * Brown (1988) J. D. Brown,  _Lower Dimensional Gravity_ (1988).
  * Teitelboim (1983) C. Teitelboim, Phys. Lett. B 126, 41 (1983).
  * Jackiw (1985) R. Jackiw, Nucl. Phys. B 252, 343 (1985).
  * Mann  _et al._ (1990) R. B. Mann, A. Shiekh, and L. Tarasov, Nucl. Phys. B 341, 134 (1990).
  * Sikkema and Mann (1991) A. E. Sikkema and R. B. Mann, Classical and Quantum Gravity 8, 219 (1991).
  * Sabín (2018) C. Sabín, New Journal of Physics 20, 053028 (2018).
  * García and Sabín (2019) J. F. García and C. Sabín, Phys. Rev. D 99, 025008 (2019).
  * So (1985) H. So, Progress of Theoretical Physics 73, 528 (1985).
  * Ford (2021) L. H. Ford, Reports on Progress in Physics 84, 116901 (2021).
  * Preskill (1990) J. Preskill, “Quantum field theory in curved spacetime,”  (1990).
  * Bunch and Davies (1978) T. S. Bunch and P. C. W. Davies, Proc. Roy. Soc. Lond. A 360, 117 (1978).
  * Barbero G.  _et al._ (2018) J. F. Barbero G., A. Ferreiro, J. Navarro-Salas, and E. J. S. Villaseñor, Phys. Rev. D 98, 025016 (2018).
  * Valatin (1958) J. G. Valatin, Il Nuovo Cimento 7, 843 (1958).
  * Bogolyubov (1958) N. N. Bogolyubov, Nuovo Cim. 7, 794 (1958).
  * Anderson  _et al._ (2018) P. R. Anderson, E. Mottola, and D. H. Sanders, Phys. Rev. D 97, 065016 (2018), arXiv:1712.04522 [gr-qc] .
  * Duncan (1978) A. Duncan, Phys. Rev. D 17, 964 (1978).
  * Fuentes  _et al._ (2010) I. Fuentes, R. B. Mann, E. Martín-Martínez, and S. Moradi, Phys. Rev. D 82, 045030 (2010).
  * Martín-Martínez and Menicucci (2012) E. Martín-Martínez and N. C. Menicucci, Classical and Quantum Gravity 29, 224003 (2012).
  * Grib  _et al._ (1994) A. A. Grib, V. M. Mostepanenko, and S. G. Mamayev (1994).
  * Haro and Elizalde (2008) J. Haro and E. Elizalde, Journal of Physics A: Mathematical and Theoretical 41, 372003 (2008).
  * Abramowitz  _et al._ (1988) M. Abramowitz, I. A. Stegun, and R. H. Romer, American Journal of Physics 56, 958 (1988).
  * Ema  _et al._ (2019) Y. Ema, K. Nakayama, and Y. Tang, Journal of High Energy Physics 2019, 60 (2019).
  * Borrajo Gutiérrez  _et al._ (2020) D. E. Borrajo Gutiérrez, J. A. R. Cembranos, L. J. Garay, and J. M. Sánchez Velázquez, Journal of High Energy Physics 2020, 69 (2020).
  * Villalba (1995) V. M. Villalba, Phys. Rev. D 52, 3742 (1995), arXiv:hep-th/9507021 .
  * Wilson and Kogut (1974) K. G. Wilson and J. Kogut, Physics Reports 12, 75 (1974).
  * Wilson (1974a) K. G. Wilson, Phys. Rev. D 10, 2445 (1974a).
  * Gattringer and Lang (2010) C. Gattringer and C. B. Lang,  _Quantum chromodynamics on the lattice_, Vol. 788 (Springer, Berlin, 2010).
  * Rothe (2012) H. J. Rothe,  _Lattice Gauge Theories : An Introduction (Fourth Edition)_, Vol. 43 (World Scientific Publishing Company, 2012).
  * Wilson (1977) K. G. Wilson, in  _New Phenomena in Subnuclear Physics_ (Springer US, 1977) pp. 69–142.
  * Bermudez  _et al._ (2010) A. Bermudez, L. Mazza, M. Rizzi, N. Goldman, M. Lewenstein, and M. A. Martin-Delgado, Phys. Rev. Lett. 105, 190404 (2010).
  * Kaplan and Sun (2012) D. B. Kaplan and S. Sun, Phys. Rev. Lett. 108, 181807 (2012).
  * Jünemann  _et al._ (2017) J. Jünemann, A. Piga, S.-J. Ran, M. Lewenstein, M. Rizzi, and A. Bermudez, Phys. Rev. X 7, 031057 (2017).
  * Bermudez  _et al._ (2018) A. Bermudez, E. Tirrito, M. Rizzi, M. Lewenstein, and S. Hands, Annals of Physics 399, 149 (2018).
  * Kuno (2019) Y. Kuno, Phys. Rev. B 99, 064105 (2019).
  * Tirrito  _et al._ (2019) E. Tirrito, M. Rizzi, G. Sierra, M. Lewenstein, and A. Bermudez, Phys. Rev. B 99, 125106 (2019).
  * Sen (2020) S. Sen, Phys. Rev. D 102, 094520 (2020).
  * Kaplan and Sen (2021) D. B. Kaplan and S. Sen, “Index theorems, generalized hall currents and topology for gapless defect fermions,”  (2021).
  * Liu  _et al._ (2016) C.-X. Liu, S.-C. Zhang, and X.-L. Qi, Annual Review of Condensed Matter Physics 7, 301 (2016).
  * Haldane (1988) F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * Haldane (2017) F. D. M. Haldane, Rev. Mod. Phys. 89, 040502 (2017).
  * Qi  _et al._ (2008) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B 78, 195424 (2008).
  * Landau (1937) L. Landau, Zh. Eksp. Teor. Fiz. 7, 19 [Phys. Z. Sowjetunion 11, 26 (1937)] (1937).
  * Wen (2007) X.-G. Wen,  _Quantum field theory of many-body systems: from the origin of sound to an origin of light and electrons_ (Oxford University Press, Oxford, 2007).
  * Thouless  _et al._ (1982) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * Nakahara (2017) M. Nakahara,  _Geometry, Topology and Physics_ (CRC Press, 2017).
  * Senthil (2015) T. Senthil, Annual Review of Condensed Matter Physics 6, 299 (2015).
  * Kitaev (2009) A. Kitaev, AIP Conference Proceedings 1134, 22 (2009).
  * Chiu  _et al._ (2016) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * Ryu  _et al._ (2010) S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, New Journal of Physics 12, 065010 (2010).
  * Su  _et al._ (1979) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698 (1979).
  * Asbóth  _et al._ (2016) J. K. Asbóth, L. Oroszlány, and A. Pályi,  _A Short Course on Topological Insulators_, Vol. 919 (2016).
  * Shankar (2018) R. Shankar, “Topological insulators – a review,”  (2018).
  * Zak (1989) J. Zak, Phys. Rev. Lett. 62, 2747 (1989).
  * Berry (1984) M. V. Berry, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 392, 45 (1984).
  * Xiao  _et al._ (2010) D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * Jackiw and Rebbi (1976) R. Jackiw and C. Rebbi, Phys. Rev. D 13, 3398 (1976).
  * Magnifico  _et al._ (2019a) G. Magnifico, D. Vodola, E. Ercolessi, S. P. Kumar, M. Müller, and A. Bermudez, Phys. Rev. D 99, 014503 (2019a).
  * Magnifico  _et al._ (2019b) G. Magnifico, D. Vodola, E. Ercolessi, S. P. Kumar, M. Müller, and A. Bermudez, Phys. Rev. B 100, 115152 (2019b).
  * Kogut and Susskind (1975) J. Kogut and L. Susskind, Phys. Rev. D 11, 395 (1975).
  * Susskind (1977) L. Susskind, Phys. Rev. D 16, 3031 (1977).
  * Creutz (1999a) M. Creutz, Phys. Rev. Lett. 83, 2636 (1999a).
  * Creutz (2001) M. Creutz, Rev. Mod. Phys. 73, 119 (2001).
  * Mazza  _et al._ (2012) L. Mazza, A. Bermudez, N. Goldman, M. Rizzi, M. A. Martin-Delgado, and M. Lewenstein, New Journal of Physics 14, 015007 (2012).
  * Velasco and Paredes (2019) C. G. Velasco and B. Paredes, “Classification of topological ladder models,” (2019), arXiv:1907.11460 .
  * Zurita  _et al._ (2021) J. Zurita, C. Creffield, and G. Platero, Quantum 5, 591 (2021).
  * Viyuela  _et al._ (2012) O. Viyuela, A. Rivas, and M. A. Martin-Delgado, Phys. Rev. B 86, 155140 (2012).
  * Viyuela  _et al._ (2014) O. Viyuela, A. Rivas, and M. A. Martin-Delgado, Phys. Rev. Lett. 112, 130401 (2014).
  * Bermudez  _et al._ (2009) A. Bermudez, D. Patanè, L. Amico, and M. A. Martin-Delgado, Phys. Rev. Lett. 102, 135702 (2009).
  * Jafari  _et al._ (2019) R. Jafari, H. Johannesson, A. Langari, and M. A. Martin-Delgado, Phys. Rev. B 99, 054302 (2019).
  * Sun and Lim (2017) N. Sun and L.-K. Lim, Phys. Rev. B 96, 035139 (2017).
  * Haller  _et al._ (2020) A. Haller, M. Rizzi, and M. Filippone, Phys. Rev. Research 2, 023058 (2020).
  * Chan  _et al._ (2022) S. M. Chan, B. Grémaud, and G. G. Batrouni, Phys. Rev. B 105, 024502 (2022).
  * Orito  _et al._ (2022) T. Orito, Y. Kuno, and I. Ichinose, Phys. Rev. B 105, 094201 (2022).
  * Orito  _et al._ (2021a) T. Orito, Y. Kuno, and I. Ichinose, Phys. Rev. B 103, L060301 (2021a).
  * Orito  _et al._ (2021b) T. Orito, Y. Kuno, and I. Ichinose, Phys. Rev. B 104, 094202 (2021b).
  * Zhang  _et al._ (2021) G.-Q. Zhang, L.-Z. Tang, L.-F. Zhang, D.-W. Zhang, and S.-L. Zhu, Phys. Rev. B 104, L161118 (2021).
  * González-Cuadra  _et al._ (2020) D. González-Cuadra, L. Tagliacozzo, M. Lewenstein, and A. Bermudez, Phys. Rev. X 10, 041007 (2020).
  * Tirrito  _et al._ (2022) E. Tirrito, M. Lewenstein, and A. Bermudez, Phys. Rev. B 106, 045147 (2022).
  * Wilson (1974b) K. G. Wilson, Phys. Rev. D 10, 2445 (1974b).
  * Kibble (1980) T. Kibble, Physics Reports 67, 183 (1980).
  * Zurek (1996) W. Zurek, Physics Reports 276, 177 (1996).
  * Zurek  _et al._ (2005) W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 95, 105701 (2005).
  * Creutz (1999b) M. Creutz, Phys. Rev. Lett. 83, 2636 (1999b).
  * Bloch  _et al._ (2012) I. Bloch, J. Dalibard, and S. Nascimbène, Nature Physics 8, 267 (2012).
  * Goldman  _et al._ (2016) N. Goldman, J. C. Budich, and P. Zoller, Nature Physics 12, 639 (2016).
  * Lewenstein  _et al._ (2007) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen(De), and U. Sen, Advances in Physics 56, 243 (2007).
  * Boada  _et al._ (2011) O. Boada, A. Celi, J. I. Latorre, and M. Lewenstein, New Journal of Physics 13, 035002 (2011).
  * Minář and Grémaud (2015) J. Minář and B. Grémaud, Journal of Physics A: Mathematical and Theoretical 48, 165001 (2015).
  * Rodríguez-Laguna  _et al._ (2017) J. Rodríguez-Laguna, L. Tarruell, M. Lewenstein, and A. Celi, Phys. Rev. A 95, 013627 (2017).
  * Celi (2017) A. Celi, The European Physical Journal Special Topics 226, 2729 (2017).
  * Kosior  _et al._ (2018) A. Kosior, M. Lewenstein, and A. Celi, SciPost Phys. 5, 061 (2018).
  * Wiese (2013) U.-J. Wiese, Annalen der Physik 525, 777 (2013).
  * Zohar  _et al._ (2015) E. Zohar, J. I. Cirac, and B. Reznik, Reports on Progress in Physics 79, 014401 (2015).
  * Bañuls  _et al._ (2020) M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik, E. Rico, L. Tagliacozzo, K. Van Acoleyen, F. Verstraete, U.-J. Wiese, M. Wingate, J. Zakrzewski, and P. Zoller, The European Physical Journal D 74, 165 (2020).
  * Bañuls and Cichy (2020) M. C. Bañuls and K. Cichy, Reports on Progress in Physics 83, 024401 (2020).
  * Aidelsburger  _et al._ (2022) M. Aidelsburger, L. Barbiero, A. Bermudez, T. Chanda, A. Dauphin, D. González-Cuadra, P. R. Grzybowski, S. Hands, F. Jendrzejewski, J. Jünemann, G. Juzeliūnas, V. Kasper, A. Piga, S.-J. Ran, M. Rizzi, G. Sierra, L. Tagliacozzo, E. Tirrito, T. V. Zache, J. Zakrzewski, E. Zohar, and M. Lewenstein, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 380, 20210064 (2022).
  * Klco  _et al._ (2022) N. Klco, A. Roggero, and M. J. Savage, Reports on Progress in Physics 85, 064301 (2022).
  * Ziegler  _et al._ (2022) L. Ziegler, E. Tirrito, M. Lewenstein, S. Hands, and A. Bermudez, Phys. Rev. Research 4, L042012 (2022).
  * Ziegler  _et al._ (2022) L. Ziegler, E. Tirrito, M. Lewenstein, S. Hands, and A. Bermudez, Annals of Physics 439, 168763 (2022), arXiv:2111.04485 .
  * Zhang and Liu (2018) L. Zhang and X.-J. Liu, in  _Synthetic Spin-Orbit Coupling in Cold Atoms_ (WORLD SCIENTIFIC, 2018) pp. 1–87.
  * Liu  _et al._ (2013) X.-J. Liu, Z.-X. Liu, and M. Cheng, Phys. Rev. Lett. 110, 076401 (2013).
  * Liu  _et al._ (2014a) X.-J. Liu, K. T. Law, and T. K. Ng, Phys. Rev. Lett. 112, 086401 (2014a).
  * Liu  _et al._ (2014b) X.-J. Liu, K. T. Law, and T. K. Ng, Phys. Rev. Lett. 113, 059901 (2014b).
  * Cazalilla and Rey (2014) M. A. Cazalilla and A. M. Rey, Reports on Progress in Physics 77, 124401 (2014).
  * Liang  _et al._ (2021) M.-C. Liang, Y.-D. Wei, L. Zhang, X.-J. Wang, H. Zhang, W.-W. Wang, W. Qi, X.-J. Liu, and X. Zhang, “Realization of qi-wu-zhang model in spin-orbit-coupled ultracold fermions,”  (2021).
  * Wu  _et al._ (2016) Z. Wu, L. Zhang, W. Sun, X.-T. Xu, B.-Z. Wang, S.-C. Ji, Y. Deng, S. Chen, X.-J. Liu, and J.-W. Pan, Science 354, 83 (2016).
  * Sun  _et al._ (2018) W. Sun, B.-Z. Wang, X.-T. Xu, C.-R. Yi, L. Zhang, Z. Wu, Y. Deng, X.-J. Liu, S. Chen, and J.-W. Pan, Phys. Rev. Lett. 121, 150401 (2018).
  * Song  _et al._ (2018) B. Song, L. Zhang, C. He, T. F. J. Poon, E. Hajiyev, S. Zhang, X.-J. Liu, and G.-B. Jo, Science Advances 4, eaao4748 (2018).
  * Jaksch  _et al._ (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
  * Hofstetter  _et al._ (2002) W. Hofstetter, J. I. Cirac, P. Zoller, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 89, 220407 (2002).
  * Jaksch and Zoller (2003) D. Jaksch and P. Zoller, New Journal of Physics 5, 56 (2003).
  * Ben Dahan  _et al._ (1996) M. Ben Dahan, E. Peik, J. Reichel, Y. Castin, and C. Salomon, Phys. Rev. Lett. 76, 4508 (1996).
  * Huang and Parker (2009) X. Huang and L. Parker, Phys. Rev. D 79, 024020 (2009).
  * Bloch  _et al._ (2008b) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008b).
  * Ket (2008) La Rivista del Nuovo Cimento 31, 247–422 (2008).
  * Greiner  _et al._ (2001) M. Greiner, I. Bloch, O. Mandel, T. W. Hänsch, and T. Esslinger, Phys. Rev. Lett. 87, 160405 (2001), arXiv:cond-mat/0105105 [cond-mat.soft] .
  * Köhl  _et al._ (2005) M. Köhl, H. Moritz, T. Stöferle, K. Günter, and T. Esslinger, Physical Review Letters 94 (2005), 10.1103/physrevlett.94.080403.
  * Kastberg  _et al._ (1995) A. Kastberg, W. D. Phillips, S. L. Rolston, R. J. C. Spreeuw, and P. S. Jessen, Phys. Rev. Lett. 74, 1542 (1995).
  * Tarruell  _et al._ (2012) L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, Nature (London) 483, 302 (2012), arXiv:1111.5020 [cond-mat.quant-gas] .
  * Alba  _et al._ (2011) E. Alba, X. Fernandez-Gonzalvo, J. Mur-Petit, J. K. Pachos, and J. J. Garcia-Ripoll, Phys. Rev. Lett. 107, 235301 (2011).
  * Cooper  _et al._ (2019) N. R. Cooper, J. Dalibard, and I. B. Spielman, Reviews of Modern Physics 91, 015005 (2019), arXiv:1803.00249 [cond-mat.quant-gas] .
  * Gradshteyn and Ryzhik (2007) I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ , seventh ed. (Elsevier/Academic Press, Amsterdam, 2007) pp. xlviii+1171, translated from the Russian, Translation edited and with a preface by Alan Jeffrey and Daniel Zwillinger, With one CD-ROM (Windows, Macintosh and UNIX).


# Upgrading gSeaGen: from MeV to PeV neutrinos

A. Garcia Soto,11footnotetext: Corresponding author. C. Distefano  P.
Kalaczyński

###### Abstract

The gSeaGen framework has been upgraded to simulate events detectable in
neutrino telescopes induced by neutrino or cosmic ray interactions. The new
version is well-suited to generate neutrino interactions at energies from a
few MeV to EeV, profiting from the latest GENIE extensions to lower and higher
energies. In addition, a brand new functionality to propagate leptons from
CORSIKA air showers has been developed. Novel features of gSeaGen will be
presented, related to simulation of PeV tau neutrinos and charged lepton
propagation, using the KM3NeT detectors as physics case.

##  1 Introduction

Under-water/ice Cherenkov neutrino detectors, also known as neutrino
telescopes, have a unique feature as they can detect neutrinos in a vast
energy range between few MeV and PeV. In particular, the KM3NeT collaboration
is building two different detectors (ORCA and ARCA) that will be sensitive to
neutrinos of MeV, GeV and PeV energies [1]. Hence, one of the most challenging
tasks for these experiments is to understand distinct aspects (e.g. neutrino
sources, transport through the planet and generation of events near the
sensitive volume) in very different energy regimes.

The gSeaGen framework [2] was designed to efficiently simulate neutrino
interactions in large volumes accounting for the aforementioned points, and
for this reason is a very suitable generator for neutrino telescopes. In its
early stages, the package was mainly developed for KM3NeT. Nevertheless, the
current version of the code is open source and can be easily customized to the
needs of other neutrino telescopes.

Fig. 1 shows a diagram illustrating the main steps used to generate events
with gSeaGen. One key feature of the code is that it allows the user to change
the properties of the media (e.g. ice, water, rock) as well as detector
locations (i.e. depth and geographic coordinates). In addition, the code
includes a driver to automatically calculate the volume in which an
interaction can produce detectable particles depending on the energy and
flavor of the incoming neutrino. It is also able to generate neutrinos coming
from diffuse, point-like and extended sources (both steady and transient).
Finally, the code relies on external packages like GENIE [3], MUSIC [4] or
PROPOSAL [5] to simulate the neutrino interactions and the propagation of
charged leptons.

Figure 1: Block diagram of gSeaGen. Blocks show the main steps that are
processed in the package and arrows illustrate the sequence to generate an
event in the detector.

In this work we describe new features that are available in the code,
particularly the extension to the high energy regime profiting from recent
developments in GENIE [6]. The generation of neutrino events above 1 TeV in a
KM3NeT ARCA-like detector geometry is discussed in Sec. 2. A comparison
between different muon propagation codes is done in Sec. 3.

##  2 Neutrino interaction

The previous version of gSeaGen [2] was limited to generating neutrinos up to
5 TeV due to GENIE’s range of validity. The newest version of GENIE includes
new modules (called HEDIS and GLRES) that compute neutrino differential cross
sections up to 1 EeV. This new version has been tested in gSeaGen and results
are presented in the following paragraphs.

In addition to the extension to generate neutrino interactions of any flavor
at high energies, the new version simulate tau leptons differently. In the
first version of gSeaGen, the outgoing tau from Charged Current interactions
was decayed on the spot and PYTHIA6 [7] was used to compute the kinematics of
the decay products. This method is a good approximation at energies below 1
TeV but is not valid at higher energies, hence several aspects for the
generation of tau neutrinos have been changed. First, the calculation of the
interaction volume is done accounting for the decay length and products of the
outgoing lepton. Second, the propagation of the tau includes energy losses
using TAUSIC [4]. Finally, the kinematics of the decay products are computed
using TAUOLA [8].

Fig. 2 shows the energy spectrum for the three neutrino species producing
detectable events in one building block of ARCA assuming a diffuse flux. Two
different DIS models implemented in HEDIS (CSMS and BGR) are compared. The
rate of down-going events is mainly dominated by the total neutrino-nucleus
cross sections, whereas up-going events are affected by Earth absorption at
high energies. The muon neutrino contribution dominates at these energies as
muons originated in CC interactions travel long distances and reach the active
volume of the detector. Above 10 TeV the tau and electron neutrino rates
differ since the tau leptons can propagate as muons before they decay.

Figure 2: Muon (left), tau (center), electron (right) neutrino and anti-
neutrino energy distributions generated with gSeaGen for up- and down-going
neutrinos. The diffuse flux shown in the legend is assumed. Two different
cross section models are compared: CSMS (dashed lines) and BGR (continuous
lines).

Another important development is related to the generation of Glashow resonant
events [9] since the latest version of GENIE allows to simulate this type of
interactions as well. Most of the studies of this channel have been focused on
trying to observe fully or partially contained showers in the detector coming
from the decay of a W boson into hadrons [10]. To simulate these events the
interaction volume can be restricted to the convex hull of the detector as
hadrons will travel only a few meters. However, when the W boson decays into a
muon or a tau the interaction volume must be significantly larger. gSeaGen can
compute different interaction volumes depending on the decay channel. Fig. 3
shows the expected rate of Glashow resonant events in one ARCA block. One can
observe that the main contribution comes from muons even though its branching
ratio is smaller than the hadronic decay.

Figure 3: Electron anti-neutrino energy distributions generated with gSeaGen
for up- and down-going neutrinos interacting via the Glashow resonance. Left:
comparison with the predictions for DIS interactions for electron neutrinos
and anti-neutrinos. Right: breakdown in the different decay channels of the W
boson.

##  3 Muon propagation

A fundamental step in the simulation of events around the detector volume is
the propagation of muons. On the one hand, muons produced in neutrino
interactions can travel long distances and produce light in the detector. On
the other hand, muon bundles from extensive air showers can also be detected
at large depths. Hence, it is crucial to properly simulate the energy losses
and scattering as muons propagate through matter. Currently, gSeaGen includes
two internal implementations to propagate muons in water/ice and rock called
PropaMuon and JPP. Additionally, the code can be linked to two external open
source programs: MUSIC [4] and PROPOSAL [5].

Fig. 4 shows the energy spectrum for the muon neutrinos producing detectable
events in one building block of ARCA assuming a diffuse flux. Two different
muon propagation codes (MUSIC and PROPOSAL) are compared for up-going and
down-going events (i.e. rock and water propagation respectively). The
agreement between the two codes is better than 5% and only subtle differences
are observed when propagating muons through rock. Recently, direct processing
of CORSIKA output files [11] has been implemented in gSeaGen. Consequently, it
is possible to propagate muons from the air showers until any custom depth of
ice or water as shown in Figure 4.

Figure 4: Left: Muon (anti-)neutrino energy distributions generated with
gSeaGen for up- and down-going neutrinos. The results obtained with MUSIC and
PROPOSAL are compared and found to be consistent between each other. Right:
Comparison of muon multiplicity distributions at the sea level and after
propagating to the ARCA-like detector geometry. The propagation is done using
PROPOSAL on atmospheric muons simulated with CORSIKA.

##  4 Conclusions

This work details the main developments of gSeaGen with respect to the
previous version. The two main new features are: possibility of generating
neutrino events up to PeV energies using the latest version of GENIE, and a
new tool that allows processing CORSIKA output files to propagate air showers
in water or ice. This upgrade is already publicly
available222https://zenodo.org/record/4766015#.YKOtspMzau4.

## Acknowledgments

A.G. acknowledges support from the European Union’s H2020-MSCA Grant Agreement
No. 101025085.

## References

  * [1] S. Adrian-Martinez et al. (KM3NeT collaboration), _Letter of intent for KM3NeT 2.0_ , J. Phys. G: Nucl. Part. Phys. 43 (2016) 084001.
  * [2] S. Aiello et al. (KM3NeT collaboration), _gSeaGen: the KM3NeT GENIE-based code for neutrino telescopes_ , Comput. Phys. Commun. 256 107477 (2020).
  * [3] C. Andreopoulos et al., _The GENIE Neutrino Monte Carlo Generator_ , Nucl. Instrum. Meth. A614 (2010) 87.
  * [4] P. Antonioli et al., _A Three-Dimensional Code for Muon Propagation through the Rock: MUSIC_ , Astropart.Phys. 7 (1997) 357-368.
  * [5] J.-H. Koehne et al., _PROPOSAL: A tool for propagation of charged leptons_ , Comput. Phys. Commun. 184 (2013) 2070-2090.
  * [6] L. Alvarez-Ruso et al., _Recent highlights from GENIE v3_ , arXiv:2106.09381 (2021).
  * [7] T. Sjöstrand et al., _PYTHIA 6.4 physics and manual_ , J. High Energy Phys. 05 (2006) 026.
  * [8] N. Davidson et al., _Universal interface of tauola: Technical and physics documentation_ , Comput. Phys. Commun. 183 (2012) 821.
  * [9] S. L. Glashow, _Resonant Scattering of Antineutrinos_ , Phys. Rev. 118 (1960) 316.
  * [10] M. G. Aartsen et al. (IceCube Collaboration), _Detection of a particle shower at the Glashow resonance with IceCube_ , Nature volume 591 (2021) 220–224.
  * [11] D. Heck et al., _CORSIKA: A Monte Carlo code to simulate extensive air showers_ , FZKA-6019 (1998).