# Pulses with Minimum Residual Intersymbol Interference for Faster than
Nyquist Signaling

Youssef Jaffal, Alex Alvarado, Senior Member, IEEE This work has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme under Grant 757791.
(Corresponding author: Youssef Jaffal.) The authors are with the Department of
Electrical Engineering, Eindhoven University of Technology, 5600 MB Eindhoven,
The Netherlands (e-mail<EMAIL_ADDRESS>a.alvarado@tue.nl).

###### Abstract

Faster than Nyquist signaling increases the spectral efficiency of pulse
amplitude modulation by accepting inter-symbol interference, where an
equalizer is needed at the receiver. Since the complexity of an optimal
equalizer increases exponentially with the number of the interfering symbols,
practical truncated equalizers assume shorter memory. The power of the
resulting residual interference depends on the transmit filter and limits the
performance of truncated equalizers. In this paper, we use numerical
optimizations and the prolate spheroidal wave functions to find optimal time-
limited pulses that achieve minimum residual interference. Compared to root
raised cosine pulses, the new pulses decrease the residual interference by an
order of magnitude, for example, a decrease by 32 dB is achieved for an
equalizer that considers four interfering symbols at $57\%$ faster
transmissions. As a proof of concept, for the $57\%$ faster transmissions of
binary symbols, we showed that using the new pulse with a $4$-state equalizer
has better bit error rate performance compared to using a root raised cosine
pulse with a $128$-state equalizer.

###### Index Terms:

Faster than Nyquist, residual intersymbol interference, prolate spheroidal
wave functions.

##  I Introduction

Faster than Nyquist (FTN) was introduced by Mazo in 1975 [1], where he found
that increasing the signaling rate of sinc pulses by $25\%$ does not affect
the minimum distance between binary coded codewords, which allows to increase
the data rates without affecting the probability of error when using an
optimum receiver. Almost 30 years later, the study of minimum distance in FTN
was extended to root raised cosine (RRC) filters in [2]. Constrained
theoretical information rates of FTN were investigated in [3] and [4], where
it was shown that FTN is beneficial from an information theoretic perspective.

As FTN signaling causes intersymbol interference (ISI), achieving the full
gains that are suggested by the analysis of the minimum distance and the
theoretical rates may require high-complexity (or impractical) receivers.
Optimal receivers with reduced complexity were proposed in [5] and [6] for
channels with ISI. Forney proposed to use a whitening filter followed by a
Viterbi algorithm in [5]. In [6], Ungerboeck proposed a modified Viterbi
algorithm that operates directly on the received symbols without using a
whitening filter.

As the complexity of the optimum equalizer increases exponentially with the
number of interfering symbols, the equalizer needs to be truncated in order to
be used in practical systems [7]. In [2], practical ways of coding and
equalization were proposed to achieve some gains of FTN, where a truncated
modified Viterbi algorithm (TMVA) was proposed to detect the transmitted
symbols. Truncated equalizers offer a trade-off between complexity and
probability of error, where the effect of the residual ISI (RISI) can be
decreased at the cost of increased complexity. The main issue of RISI is the
error floor at high signal to noise ratio, since the signal to RISI ratio
remains constant [8]. Some approaches to mitigate the effect of RISI were
proposed in the literature such as the decision feedback cancellation [8].
However, such methods increase the complexity of the system. The RISI depends
on the used pulse shaping filter, for example, [8, Table 1] shows the ISI
coefficients that are obtained from the RRC pulses.

In this paper, we take a different approach and we propose to mitigate RISI by
using pulse shaping filters that minimize the RISI power. Our proposal does
not increase the complexity of FTN systems as it requires only replacing the
RRC pulse by the pulses that are tailored for FTN regimes. We find optimal
time-limited pulses that minimize the RISI power while having a duration and
an out-of-band energy (OOBE) that are equal to the ones of the truncated RRC.
We propose a numerical method that makes use of the prolate spheroidal wave
functions (PSWFs). The PSWFs form a complete orthonormal set for time-limited
functions and they were used in [9] and [4] to design optimal time-limited
pulses. The main contribution of this paper is to introduce new pulses that
have the ability to improve the performance and the complexity of the receiver
in FTN systems.

The paper is organized as follows: in Sec. II we present the considered system
model. In Sec. III we formulate our optimization problem and we propose a
numerical solution to search for the optimal pulses. In Sec. IV we present the
results of the numerical optimizations and simulations. And we finally
conclude in Sec. V.

##  II System Model

$+$$n(t)$$p(t)$$p(-t)$TMVA
Equalizer$x(t)$TransmitterReceiver$-T$$T$$2T$$0$$A_{-1}$$A_{0}$$A_{1}$$A_{2}$$t$$y(t)\,$$\tilde{A}_{k}\quad$sample
at $t=kT$ Figure 1: System model under consideration. The transmit filter has
a duration $T_{s}$ (see (1)) and a bandwidth $W$ (see (3)), and the receiver
is based on the TMVA from [2].

We consider the baseband pulse amplitude modulation system illustrated in Fig.
1. The transmitted symbols $\left\\{A_{k}\right\\}_{k\in\mathbb{Z}}$ are
independent and identically distributed and are chosen from a constellation of
size $M$ that has a zero mean and a unit average energy. After pulse shaping,
the transmitted signal is

$x(t)=\sum_{l\in\mathbb{Z}}A_{l}\,p(t-lT)$

where $T$ is the modulation interval. The noise $n(t)$ is additive white
Gaussian with mean zero and power spectral density (PSD) $N_{0}/2$.

The transmit filter $p(t)$ is assumed to be time-limited to $T_{s}$ seconds,
i.e.,

$p(t)=0\text{ for }t\notin[-T_{s}/2,T_{s}/2],$ (1)

and to have a unit-norm (or unit energy)

$\int_{-T_{s}/2}^{T_{s}/2}p^{2}(t)dt=1.$ (2)

Being time-limited, $p(t)$ cannot be strictly band-limited. We define the
bandwidth of $p(t)$ by the value of $W$ that satisfies

$\int_{-W}^{W}\left|P(f)\right|^{2}df=1-\epsilon,$ (3)

where $\epsilon$ is the OOBE, and $P(f)=\int_{-T_{s}/2}^{T_{s}/2}p(t)e^{-j2\pi
ft}dt$ is the Fourier transform of $p(t)$.

At the receiver, the output of the matched filter $p(-t)$ is

$y(t)=\sum_{l\in\mathbb{Z}}A_{l}\,h(t-lT)+n(t)*p(-t)$

where $*$ is the convolution operation and

$h(t)=p(t)*p(-t)=\int_{-T_{s}/2}^{T_{s}/2}p(\tau)p(\tau-t)d\tau.$ (4)

Since $p(t)$ is unit norm, then $h(0)=1$. Sampling the signal $y(t)$ at $t=kT$
yields

$\tilde{A}_{k}=A_{k}+\underbrace{\sum_{l\neq
0}A_{k-l}\,h(lT)}_{\text{ISI}}+\underbrace{n(t)*p(-t)|_{t=kT}}_{\text{colored
noise}},$ (5)

where the noise term is a zero-mean Gaussian process with variance $N_{0}/2$,
and PSD

$\frac{N_{0}}{2T}\sum_{m\in\mathbb{Z}}{\left|P\left(\frac{f-m}{T}\right)\right|^{2}}.$
(6)

In this paper, we consider pulses and modulation intervals that are capable of
increasing the data rates by accepting ISI (see (5)). Note that if the pulse
$p(t)$ does not satisfy the Nyquist criterion, then the PSD in (6) is not
flat, and hence the noise at the input of the equalizer is colored. The role
of the equalizer in Fig. 1 is to detect the transmitted sequence of symbols in
the presence of ISI and the colored noise.

Here we consider the equalizer to be the TMVA proposed in [2]. In this case,
the equalizer considers $2L$ interfering symbols and treats the RISI as noise.
The complexity of the equalizer is proportional to $M^{L}$ and its performance
depends on the variance of the RISI [2]. The RISI is given by

$\text{RISI}=\sum_{|l|>L}A_{k-l}\,h(lT),$

which has zero mean and variance

$\sigma^{2}_{\text{RISI}}=\sum_{|l|>L}h^{2}(lT).$ (7)

In view of (4) and (7), the complexity and performance of the equalizer depend
on the used pulse $p(t)$. For a given $L$ (i.e., for a fixed complexity), the
performance is maximized when using a pulse that minimizes
$\sigma^{2}_{\text{RISI}}$.

##  III Problem Formulation and Proposed Solution

We consider the problem of finding the optimal pulse $p(t)$ that minimizes
$\sigma^{2}_{\text{RISI}}$ in (7). We formulate the optimization problem as
follows

$p^{*}(t)=\underset{p(t)}{\text{argmin}}\sum_{|l|>L}\left(p(t)*p(-t)|_{t=lT}\right)^{2},$
(8)

where $p(t)$ should satisfy (1), (2), and (3), and where minimizing the
objective function of (8) minimizes (7). The constraints are used to restrict
the search to unit norm time-limited pulses that have an OOBE equal to
$\epsilon$. In the following, we propose a numerical method to find the
optimal solution of (8).

Since $p(t)$ has unit energy and is $T_{s}$-seconds time-limited, it can be
written as a linear combination of normalized truncated PSWFs111Appendix -A
includes the definition and some properties of the PSWFs.

$p(t)=\sum\limits_{i=0}^{\infty}\alpha_{i}\frac{D\varphi_{c,i}(t)}{\sqrt{\lambda_{c,i}}},$
(9)

where $c=2T_{s}W$ is the time-bandwidth product of the PSWFs and
$\lambda_{c,i}$ is the eigenvalue of the $i^{th}$ PSWF. Note that the PSWFs
and their eigenvalues do not have closed form expressions. In this work we
generate them numerically using the software in [10].

In view of (9), solving (8) is equivalent to optimizing the coefficients
$\\{\alpha_{i}\\}_{i=0}^{\infty}$. Optimizing numerically over the infinite
set $\\{\alpha_{i}\\}_{i=0}^{\infty}$ is not feasible, so we approximate the
optimal solution using a finite subset of the truncated PSWFs and we optimize
over the set $\\{\alpha_{i}\\}_{i=0}^{N-1}$. We define the approximation error
to be the energy of
$\sum\limits_{i=N}^{\infty}\alpha_{i}\frac{D\varphi_{c,i}(t)}{\sqrt{\lambda_{c,i}}}$,
and by Property .2 in Appendix -A, the approximation error can be made smaller
than $\epsilon$ using finite values of $N$. The resulting optimization problem
at hand is

$\begin{aligned}
\\{\alpha^{*}_{i}\\}_{i=0}^{N-1}&=\,\,\underset{\\{\alpha_{i}\\}_{i=0}^{N-1}}{\text{argmin}}&&\sum_{|l|>L}\left(\sum\limits_{i=0}^{N-1}\alpha_{i}\frac{D\varphi_{c,i}(t)}{\sqrt{\lambda_{c,i}}}*\right.\\\
&&&\quad\left.\left.\sum\limits_{i=0}^{N-1}\alpha_{i}\frac{D\varphi_{c,i}(-t)}{\sqrt{\lambda_{c,i}}}\right|_{t=lT}\right)^{2}\\\
&\text{subject to}&&\sum_{i=0}^{N-1}\alpha_{i}^{2}=1\\\
&&&\sum\limits_{i=0}^{N-1}\alpha_{i}^{2}\lambda_{c,i}=1-\epsilon\end{aligned},$
(10)

where (1) is satisfied by writing $p(t)$ as a combination of truncated PSWFs,
and the equivalent constraints of (2) and (3) are obtained using (11) and (12)
respectively.

##  IV Numerical Results

We chose the RRC pulses as a baseline for our system model. RRC pulses are
characterized by their roll-off factor $\beta$, where the modulation interval
$T=\frac{1+\beta}{2W}$ avoids the ISI in (5). As RRC pulses have infinite time
duration, they are truncated to some finite time window $T_{s}$, which results
in non-zero ISI and some OOBE. We selected $\beta=0.1$ and we considered a
truncation window $T_{s}=\frac{15}{2W}$, where the obtained OOBE is
$\epsilon\approx 4.4\times 10^{-4}$. In Fig. 2 we show the obtained
$\sigma^{2}_{\text{RISI}}$ (which is defined in (7)), where $T=\frac{1.1}{2W}$
is the Nyquist signaling rate, and we refer to $T<\frac{1.1}{2W}$ as the FTN
signaling rates. For example, using $T=\frac{0.7}{2W}$ provides
$\left(\frac{1.1}{0.7}-1\right)\times 100\approx 57\%$ faster transmissions
which increases the data rates by $57\%$. Note that $L=0$ implies that no
equalizer is used at the receiver, and as $L$ increases, the power of the RISI
decreases since the equalizer considers $2L$ interfering symbols.

$\frac{0.5}{2W}$$\frac{0.6}{2W}$$\frac{0.7}{2W}$$\frac{0.8}{2W}$$\frac{0.9}{2W}$$\frac{1}{2W}$$\frac{1.1}{2W}$$-80$$-60$$-40$$-20$$0$$L=1$$L=2$$L=4$$L=0$$L=4$$L=2$$L=1$$\left(\frac{0.61}{2W},-14\right)$$32$
dB$T$ (s)$\sigma^{2}_{\text{RISI}}$ (dB)RRC$p^{*}(t)$ Figure 2: Performance of
obtained $p^{*}(t)$ and truncated RRC with $\beta=0.1$, for
$T_{s}=\frac{15}{2W}$ and $\epsilon=4.4\times 10^{-4}$.

To have a fair comparison with the truncated RRC pulses, we use the parameters
$T_{s}=\frac{15}{2W}$ and $\epsilon=4.4\times 10^{-4}$ in our optimization
problem (10), where the time-bandwidth product of the PSWFs is $c=2T_{s}W=15$.
We use the MATLAB optimization function fmincon to search for the optimal
solution of (10). We performed the optimization using different values of $N$
(the number of used truncated PSWFs). We found that the solutions converge for
$N=22$ where using $N=30$ did not improve the achieved
$\sigma^{2}_{\text{RISI}}$, which is consistent with the analysis in Sec. III
where we showed that the optimal solution can be approximated using a finite
value of $N$. In Fig. 2 we present the obtained $\sigma^{2}_{\text{RISI}}$ of
the optimal pulses for $\epsilon=4.4\times 10^{-4}$,
$T\in\left[\frac{0.5}{2W},\frac{1+\beta}{2W}\right]$ with a step of
$\frac{0.01}{2W}$, and $L\in\\{0,1,2,4\\}$. For $L=0$, the truncated RRC pulse
is very close to the optimal one. For $L>0$, the obtained pulses show a large
decrease in the RISI variance. This allows for example to use the optimal
pulses with $L=1$ instead of the truncated RRC with $L=4$ for $T\geq
0.61/(2W)$ and $\sigma^{2}_{\text{RISI}}\leq-14$ dB. We also considered
comparing with the RRC pulse with $\beta=0.2$, where the obtained OOBE is
$\approx 9.5\times 10^{-5}$ for $T_{s}=\frac{15}{2W}$. We do not include the
obtained results here since they have the same pattern as the ones in Fig. 2.

In Fig. 3a we show the obtained optimal pulse and the RRC pulse at ‘$\square$’
in Fig. 2, i.e. for $L=2$ and $T=0.7/(2W)$. We also show their frequency
spectrum in Fig. 2b. In Fig. 4 we show the corresponding $|h(t)|^{2}$ of the
pulses of Fig. 3, where $h(t)$ is defined in (4). In table I we show the
coefficients of the PSWFs for the optimal pulse $p^{*}(t)$ of Fig. 3, where
the presented $\\{\alpha_{i}^{*}\\}_{i=0}^{N-1}$ are the solution of (10) for
the design parameters $\epsilon=4.4\times 10^{-4}$, $T=0.7/(2W)$, and $L=2$.
We found that the optimal pulse changes (i.e.
$\\{\alpha_{i}^{*}\\}_{i=0}^{N-1}$ change) when changing the value of one of
the design parameters ($\epsilon$, $T$, and $L$). We also noticed that the
obtained optimal pulses are symmetric, which means that we can write them as a
linear combination of only symmetric truncated PSWFs, i.e., the ones with even
index $i$ (see Appendix -A). Therefore we can reduce the complexity of
numerical optimization by optimizing (10) over the $\alpha_{i}$ with even
index $i$, while fixing $\alpha_{i}=0$ for odd indices.

(a)(b)$\frac{-7.5}{2W}$$\frac{-5}{2W}$$\frac{-2.5}{2W}$$0$$\frac{2.5}{2W}$$\frac{5}{2W}$$\frac{7.5}{2W}$$0$$\sqrt{\frac{W}{2}}$$\sqrt{2W}$$t$ (s)$p(t)$RRC$p^{*}(t)$$-3W$$-2W$$-W$$0$$W$$2W$$3W$$\frac{10^{-8}}{2W}$$\frac{10^{-6}}{2W}$$\frac{10^{-4}}{2W}$$\frac{10^{-2}}{2W}$$\frac{10^{0}}{2W}$$f$ (Hz)$|P(f)|^{2}$RRC$p^{*}(t)$ Figure 3: (a) Optimal pulse versus RRC at ‘$\square$’ in Fig. 2, and (b) their corresponding $|P(f)|^{2}$. (a)(b)$\frac{-15}{2W}$$\frac{-10}{2W}$$\frac{-5}{2W}$$0$$\frac{5}{2W}$$\frac{10}{2W}$$\frac{15}{2W}$$10^{-10}$$10^{-5}$$10^{0}$$t$ (s)$|h(t)|^{2}$$p^{*}(t)$Considered by the equalizerIgnored by the equalizer$\frac{-15}{2W}$$\frac{-10}{2W}$$\frac{-5}{2W}$$0$$\frac{5}{2W}$$\frac{10}{2W}$$\frac{15}{2W}$$10^{-10}$$10^{-5}$$10^{0}$$t$ (s)$|h(t)|^{2}$RRCConsidered by the equalizerIgnored by the equalizer Figure 4: (a) $|h(t)|^{2}$ corresponding to $p^{*}(t)$ in Fig. 3, and (b) $|h(t)|^{2}$ corresponding to RRC of Fig. 3. The equalizer considers $2L=4$ interfering symbols. TABLE I: The coefficients $\\{\alpha_{i}^{*}\\}_{i=0}^{21}$ for the optimal pulse $p^{*}(t)$ of Fig. 3. $\alpha^{*}_{0}$ | $\alpha^{*}_{2}$ | $\alpha^{*}_{4}$ | $\alpha^{*}_{6}$ | $\alpha^{*}_{8}$ | $\alpha^{*}_{10}$ | $\alpha^{*}_{12}$ | $\alpha^{*}_{14}$ | $\alpha^{*}_{16}$ | $\alpha^{*}_{18}$ | $\alpha^{*}_{20}$ |  $\alpha^{*}_{i}$ for odd $i$
---|---|---|---|---|---|---|---|---|---|---|---
$0.8053$ | $-0.442$ | $0.2923$ | $-0.1996$ | $0.136$ | $-0.0905$ | $0.0562$ | $-0.03$ | $0.0107$ | $0.00014$ | $0.0011$ | $0$

To conclude, we performed numerical simulations to evaluate the bit error rate
(BER) of our system. We considered a binary input constellation (i.e.,
$A_{k}\in\\{-1,1\\}$), the modulation interval $0.7/(2W)$, and the pulses of
Fig. 3. For the optimal pulse we used a TMVA equalizer with $2^{L}=4$ states,
and we used $4$-state, $16$-state, and $128$-state TMVA equalizers for the
truncated RRC pulse (i.e., using $L=2$, $L=4$, and $L=7$). Fig. 5 shows that
the obtained optimal pulse with $4$-state TMVA outperforms the truncated RRC
with $128$-state TMVA in terms of BER, where $E_{b}$ is the average energy per
bit ($E_{b}=1$ in our case) and $N_{0}$ determines the PSD of the noise.
Hence, the optimal pulse can simultaneously improve the BER performance and
reduce the complexity of the receiver. Note that the truncated RRC suffers
from an error floor where the RISI limits the BER for very low noise power (or
equivalently high $E_{b}/N_{0}$).

$6$$10$$14$$18$$22$$26$$30$$34$$38$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$Error
floor caused by RISI$E_{b}/N_{0}$ (dB)BERRRC, $4$-state TMVARRC, $16$-state
TMVARRC, $128$-state TMVA$p^{*}(t)$, $4$-state TMVA Figure 5: BER simulations
of the pulses of Fig. 3(a).

##  V Conclusions

The main conclusion of this paper is that optimizing the pulse shaping filters
is essential for the practical implementations of FTN systems, where we showed
that the optimal pulses can simultaneously improve the performance and
decrease the complexity of truncated equalizers. We also found that there is
no universal optimal filter, where we obtained different optimal pulses for
the different design parameters: the modulation interval, the out-of-band
energy, and the number of interfering terms considered by the equalizer.

###  -A The Prolate Spheroidal Wave Functions

For any positive time-bandwidth product $c=2T_{s}W$, the PSWFs
$\\{\varphi_{c,i}(t)\\}_{i=0}^{\infty}$ are the normalized solutions of the
integral equation [11]222Note that our definition of the time-bandwidth
product is different than the one in [11]. For example, $\varphi_{1,i}(t)$
here is the same as $\varphi_{\frac{\pi}{2},i}(t)$ in [11].

$\lambda_{c,i}\,\varphi_{c,i}(t)=\int_{-T_{s}/2}^{T_{s}/2}\frac{\sin{2\pi
W(t-s)}}{\pi(t-s)}\varphi_{c,i}(s)\,ds,$

where $\lambda_{c,i}$ is the eigenvalue of $\varphi_{c,i}(t)$. The PSWFs are
real band-limited continuous functions that are symmetric when $i$ is even and
antisymmetric when $i$ is odd.

We denote by $D\varphi_{c,i}(t)$ the truncated PSWF in the time window
$[-T_{s}/2,T_{s}/2]$. The energy of $D\varphi_{c,i}(t)$ is equal to
$\lambda_{c,i}$ [11]. Hence the normalized truncated PSWFs are

$\left\\{\frac{D\varphi_{c,i}(t)}{\sqrt{\lambda_{c,i}}}\right\\}_{i=0}^{\infty},$

and they form a complete orthonormal set for $T_{s}$-seconds time-limited
functions that have finite energies.

Consider a unit norm $T_{s}$-seconds time-limited signal $p(t)$, then $p(t)$
can be uniquely written as a linear combination of the normalized truncated
PSWFs as given in (9), i.e.
$p(t)=\sum_{i=0}^{\infty}\alpha_{i}\frac{D\varphi_{c,i}(t)}{\sqrt{\lambda_{c,i}}}$.
By Parseval, the energy of $p(t)$ depends only on the coefficients
$\left\\{\alpha_{i}\right\\}_{i=0}^{\infty}$ and it is equal to

$\int_{-T_{s}/2}^{T_{s}/2}|p(t)|^{2}dt=\sum_{i=0}^{\infty}\alpha_{i}^{2}.$
(11)

The eigenvalue $\lambda_{c,i}$ is equal to the energy of
$\frac{D\varphi_{c,i}(t)}{\sqrt{\lambda_{c,i}}}$ in the bandwidth $[-W,\,W]$,
and in the following property we derive the energy of $p(t)$ in the bandwidth
$[-W,\,W]$.

######  Property .1.

The energy of $p(t)$ in the frequency band $[-W,\,W]$ is equal to

$\int_{-W}^{W}|P(f)|^{2}df=\sum_{i=0}^{\infty}\alpha_{i}^{2}\lambda_{c,i}$
(12)

where $P(f)$ is the Fourier transform of $p(t)$.

###### Proof.

The Fourier transform of a truncated PSWF is equal to a scaled version of the
PSWF [12, eq. (4)]. Therefore the Fourier transforms of the truncated PSWFs
are orthogonal over $[-W,\,W]$, and hence the left-hand side of (12) is equal
to the sum of the energies of
$\alpha_{i}\frac{D\varphi_{c,i}(t)}{\sqrt{\lambda_{c,i}}}$ in the band
$[-W,\,W]$. ∎

The eigenvalue $\lambda_{c,i}$ decreases as $i$ increases and
$1>\lambda_{c,i}>\lambda_{c,i+1}>0$ for all $i\geq 0$. For relatively low
index $i$ (i.e., $i<<c$), the eigenvalues are very close to $1$. And as $i$
increases beyond $c$, $\lambda_{c,i}$ decreases exponentially towards zero. As
an example, Fig. 6 shows the eigenvalues for $c=15$ where a logarithmic scale
is used on the y-axis. The exponential decay of $\lambda_{c,i}$ allows to
approximate the sum in (12) using a finite number of terms, and as shown in
the following property, it also allows to approximate $p(t)$ using a finite
number of truncated PSWFs if the OOBE of $p(t)$ is very small.

$0$$5$$10$$15$$20$$25$$30$$35$$10^{-20}$$10^{-10}$$10^{0}$$i$$\lambda_{15,i}$
Figure 6: Eigenvalues of the PSWFs for $c=15$ and $i=0,1,\dots,35$

######  Property .2.

If $\int_{-W}^{W}|P(f)|^{2}df=1-\epsilon$, then there exists a finite integer
$N$ where an upper bound on the energy of
$\sum_{i=N}^{\infty}\alpha_{i}\frac{D\varphi_{c,i}(t)}{\sqrt{\lambda_{c,i}}}$
is well approximated by $\epsilon$.

###### Proof.

Using (11) and (12), the OOBE is given by

$\sum\limits_{i=0}^{\infty}\alpha_{i}^{2}(1-\lambda_{c,i})=\epsilon.$

For any integer $N>c$, $(1-\lambda_{c,i})>(1-\lambda_{c,N-1})>0$ $\forall
i\geq N$. For sufficiently large $N$, $(1-\lambda_{c,N-1})\approx 1$ since
$\lambda_{c,N-1}$ decays exponentially with $N$ (see the example in Fig. 6).
Then

$\displaystyle\epsilon$
$\displaystyle=\sum\limits_{i=0}^{N-1}\alpha_{i}^{2}(1-\lambda_{c,i})+\sum\limits_{i=N}^{\infty}\alpha_{i}^{2}(1-\lambda_{c,i})$
$\displaystyle\epsilon$
$\displaystyle>\sum\limits_{i=N}^{\infty}\alpha_{i}^{2}(1-\lambda_{c,N-1})$
$\displaystyle\epsilon$
$\displaystyle\approx\epsilon/(1-\lambda_{c,N-1})>\sum\limits_{i=N}^{\infty}\alpha_{i}^{2},$

where $\sum\limits_{i=N}^{\infty}\alpha_{i}^{2}$ is the energy of
$\sum_{i=N}^{\infty}\alpha_{i}\frac{D\varphi_{c,i}(t)}{\sqrt{\lambda_{c,i}}}$.
∎

## References

  * [1] J. E. Mazo, “Faster-than-Nyquist signaling,” The Bell System Technical Journal, vol. 54, no. 8, pp. 1451–1462, Oct.1975.
  * [2] A. D. Liveris and C. N. Georghiades, “Exploiting faster-than-Nyquist signaling,” IEEE Transactions on Communications, vol. 51, no. 9, pp. 1502–1511, Sept. 2003.
  * [3] F. Rusek and J. B. Anderson, “Constrained capacities for faster-than-Nyquist signaling,” IEEE Transactions on Information Theory, vol. 55, no. 2, pp. 764–775, Feb. 2009.
  * [4] Y. Jaffal and I. Abou-Faycal, “Achievable rates using PAM time-limited pulses over band-limited channels: From Nyquist to FTN,” in 2019 IEEE Wireless Communications and Networking Conference (WCNC), pp. 1–6, IEEE, Oct. 2019.
  * [5] G. Forney, “Maximum-likelihood sequence estimation of digital sequences in the presence of intersymbol interference,” IEEE Transactions on Information theory, vol. 18, no. 3, pp. 363–378, May 1972.
  * [6] G. Ungerboeck, “Adaptive maximum-likelihood receiver for carrier-modulated data-transmission systems,” IEEE transactions on communications, vol. 22, no. 5, pp. 624–636, May 1974.
  * [7] P. McLane, “A residual intersymbol interference error bound for truncated-state Viterbi detectors,” IEEE Transactions on Information Theory, vol. 26, no. 5, pp. 548–553, Sept. 1980.
  * [8] A. D. Liveris and C. N. Georghiades, “Fundamentals of faster-than-nyquist signaling,” in 5G Mobile Communications, pp. 203–222, Springer, 2017.
  * [9] P. Halpern, “Optimum finite duration Nyquist signals,” IEEE Transactions on Communications, vol. 27, no. 6, pp. 884–888, Jun. 1979.
  * [10] R. Adelman, N. A. Gumerov, and R. Duraiswami, “Software for Computing the Spheroidal Wave Functions Using Arbitrary Precision Arithmetic,” arXiv preprint arXiv:1408.0074, Aug. 2014.
  * [11] D. Slepian and H. O. Pollak, “Prolate Spheroidal Wave Functions, Fourier Analysis and Uncertainty—I,” Bell System Technical Journal, vol. 40, no. 1, pp. 43–63, Jan. 1961.
  * [12] Y. Jaffal and I. Abou-Faycal, “Time-Limited Codewords over Band-Limited Channels: Data Rates and the Dimension of the WT Space,” Entropy, vol. 22, no. 9, p. 924, Aug. 2020.


# Sub-MeV Dark Matter Detection with Bilayer Graphene

Anirban Das  0000-0002-7880-9454 Jiho Jang  0000-0002-0687-891X and  Hongki
Min  0000-0001-5043-2432

###### Abstract

After decades of effort in a number of dark matter direct detection
experiments, we still do not have any conclusive signal yet. Therefore it is
high time we broaden our horizon by building experiments with increased
sensitivity for light dark matter, specifically in the sub-MeV mass regime. In
this paper, we propose a new detector material, a bilayer stack of graphene
for this purpose. Its voltage-tunable low energy sub-eV electronic band gap
makes it an excellent choice for the detector material of a light dark matter
search experiment. We compute its dielectric function using the random phase
approximation and estimate the projected sensitivity for sub-MeV dark matter-
electron scattering and sub-eV dark matter absorption. We show that a bilayer
graphene dark matter detector can have competitive sensitivity as other
candidate target materials, like a superconductor, in this mass regime. The
dark matter scattering rate in bilayer graphene is also characterized by a
daily modulation from the rotation of the Earth which may help us mitigate the
backgrounds in a future experiment.

##  1 Introduction

The search for sub-MeV particle dark matter (DM) has become the next frontier
in the direct detection experiments. Most of the direct detection experiments
in the last few decades have concentrated in heavier WIMP-like DM regime above
a GeV, but have failed to yield any conclusive signal yet. On the other hand,
the lack of experiments in the lighter mass regime have left a vast parameter
space unexplored till date. However, particle DM with sub-MeV mass can be
produced with the observed relic abundance in several theoretical scenarios,
such as strongly interacting DM [1], freeze-in or freeze-out DM [2, 3, 4, 5,
6], and other nonstandard DM scenarios [7]. The direct detection experiments
looking for heavier DM use nuclear and electronic recoils or excitations in
the detector from DM scattering, and have larger energy threshold typically
above $\mathcal{O}(10)\,\mathrm{eV}$. However, DM being nonrelativistic with a
typical velocity $240\,\mathrm{km\,s^{-1}}$ in the Milky Way halo, its kinetic
energy is about six orders of magnitude smaller than its mass. As a result,
those experiments lose sensitivity for lighter sub-MeV DM scattering as the
energy transfer $\lesssim\mathcal{O}(\mathrm{few})\,\mathrm{eV}$ falls short
of their threshold.

Therefore new ideas are needed for light DM search. In the past few years,
several target materials and experimental techniques have been suggested for
light DM search using electronic excitations of lower energy some of which
have already been realized in experiments. A few examples are athermal phonon
sensors made of semiconductors [8, 9, 10, 11, 12, 13], single charge detectors
[14, 15, 16], quantum defects in diamond and sillicon carbide [17, 18],
superfluid helium [19, 20, 21], magnetic microcalorimeters [22],
superconducting devices [23, 24, 25, 26, 27, 28], Migdal effect in
semiconductor [29], three-dimensional Dirac materials [30, 31, 32], quantum
dots [33] etc. We encourage the reader to refer to Ref. [34] for a detailed
compilation of such ideas. While these studies have broadened the future
outlook of light DM search, each detection technology comes with its unique
signal readout and background challenges that are needed to be overcome. Also,
some of the experiments have been seeing a background signals in the lower
energy channels below a few $100\,\mathrm{eV}$ that is limiting the progress
toward achieving lower energy threshold [35]. Hence, it is important to
explore various novel materials to optimize the future sensitivity in an
experiment.

Apart from the halo DM, low threshold devices are also needed to probe Earth-
bound thermalized DM population near the surface of the Earth. Such a
population of DM is predicted in a class of models where DM interacts
relatively strongly with the ordinary matter [36, 37, 38, 39, 40, 41]. In such
a scenario, a fraction of the DM coming from the Milky Way halo will scatter
multiple times in the Earth’s crust, lose energy, and get gravitationally
captured. These DM particles will thermalize with the local temperature which
is about $300\,\mathrm{K}$ near the Earth’s surface. This translates to about
$26\,\mathrm{meV}$ for the typical kinetic energy of a DM particle. Therefore,
detecting this thermalized DM population asks for either detectors with energy
threshold below $\lesssim 100\,\mathrm{meV}$, or ways to boosting them to
higher energy, or annihilation into Standard Model particles [42, 43, 44, 27,
45, 46, 47].

In this paper, we explore the possibility of using bilayer graphene (BLG) as
the target material to detect light DM. Graphene is a two-dimensional one-atom
thick hexagonal lattice of carbon atoms with unique electronic properties
[48]. In monolayer graphene, the electrons have a linear dispersion relation
in the low energy regime. The valence and conduction bands of an intrinsic
monolayer graphene touch each other at a single point in the Brillouin zone,
giving rise to its semimetallic properties with point-like Fermi surface.
Previously, monolayer graphene has been suggested as a detector material for
sub-GeV mass DM search [49, 50, 51, 52, 53]. The two-dimensional nature of
monolayer graphene was used in some of these studies to achieve directional
sensitivity of the DM scattering in the detector which may be helpful in
isolating the DM signal from the background in a future experiment. Most of
them depend on the removal of an electron from the graphene sheet, and a
minimum energy of at least a few eV is needed for a $\pi$-band electron
ejection. However, detecting sub-MeV DM requires sub-eV energy threshold
target making monolayer graphene unsuitable for this purpose. Moreover, the
electronic bands in a monolayer graphene do not have any energy gap which can
be set as a threshold. Having an energy gap $E_{g}$ is crucial to suppress the
thermal fluctuation noise in an experiment which can be done by cooling down
the setup below the energy gap ($k_{B}T<E_{g}$). Therefore, achieving a sub-eV
threshold for sub-MeV DM requires us to venture beyond monolayer graphene.

Bernal stacked (AB stacked) BLG consists of two coupled monolayer graphene
sheets, where half of the carbon atoms in the first layer are positioned
directly above the hexagonal center of the second layer, while the remaining
carbon atoms align directly above the carbon atoms in the second layer. The
interlayer coupling between the electrons from the two layers gives BLG
distinct electronic properties that are absent in monolayer graphene [54, 55,
56]. However, most importantly, a _tunable_ energy gap can be opened up to
$\sim 300\,\mathrm{meV}$ between the valence and conduction bands by
externally applying a gate voltage across the BLG surface [48, 57, 58, 59,
60]. Then electronic excitation between the bands is possible for DM
scattering with sufficient energy deposition. We compute the dielectric
function of BLG from the model Hamiltonian in the low-energy continuum limit,
and use it to compute the DM scattering and absorption rates for two example
values of sub-eV energy gap: $50$ and $100\,\mathrm{meV}$. The excited
electrons can be detected using a charge read-out technique as in other
experiments with a solid state target. We show that BLG can have good
sensitivity for DM-electron scattering cross section and for dark photon DM
mixing parameter, even with a relatively small exposure. Moreover, the
intrinsic 2D nature of BLG makes it a suitable material to measure the daily
modulation of DM signal expected due to Earth’s rotation about its axis. This
modulation can also potentially help us reject the background sources that do
not share the same time variation.

This paper is organized as follows. In section 2, we give an introduction to
the electronic properties of BLG and the effects of applying a gate voltage.
Then in section 3, we compute its dielectric function in the isotropic limit
and compute the sensitivity of a BLG detector for DM-electron scattering cross
section and vector DM absorption. Finally, we show and discuss the results in
section 4, and conclude in section 5. We use the natural units
$\hbar=c=\epsilon_{0}=1$ throughout this paper.

##  2 Electronic properties of bilayer graphene

BLG is a two-dimensional material consisting of two stacked layers of
graphene, each being a single layer of carbon atoms arranged in a hexagonal
lattice. BLG can be stacked in various ways, including AA stacked graphene
[61, 62], twisted BLG [63, 64, 65, 66], and so on. Among these, we investigate
a prominent type known as Bernal stacked (AB stacked) BLG [67, 55] for a dark
matter detector. From now on, unless otherwise specified, referring to BLG
will mean Bernal stacked (AB stacked) bilayer graphene. Figure 1 shows the
structure of monolayer and bilayer graphene in both real space and reciprocal
lattice space. In both cases, the lattice structure is described by the
primitive lattice vectors

$\displaystyle\boldsymbol{a}_{1}=a\left(\frac{1}{2},\frac{\sqrt{3}}{2}\right),\quad\boldsymbol{a}_{2}=a\left(\frac{1}{2},-\frac{\sqrt{3}}{2}\right),$
(2.1)

where $a=2.46\ \mathrm{\AA}$ is the lattice constant. The distance between
adjacent carbon atoms is $a_{\mathrm{CC}}=a/\sqrt{3}=1.42\ \mathrm{\AA}$. Each
layer of monolayer graphene has two inequivalent sublattices, labeled by A
(red) and B (blue) in figure 1. Therefore, in the case of bilayer graphene,
there are four sublattices in the unit cell. The corresponding reciprocal
lattice vectors $\boldsymbol{b}_{1}$ and $\boldsymbol{b}_{2}$ satisfy the
condition $\boldsymbol{a}_{i}\cdot\boldsymbol{b}_{j}=2\pi\delta_{ij}$ and are
given as

$\displaystyle\boldsymbol{b}_{1}=\frac{2\pi}{a}\left(1,\frac{1}{\sqrt{3}}\right),\quad\boldsymbol{b}_{2}=\frac{2\pi}{a}\left(1,-\frac{1}{\sqrt{3}}\right).$
(2.2)

Figure 1: (_Left_) Crystal structure of bilayer graphene. The unit cell is
highlighted in gray and contains four sublattices: A1, B1, A2, and B2. The
numbers 1 and 2 represent the layer indices for the top (empty circles) and
bottom (filled circles) layers, respectively. (_Center_) A schematic
illustration of the stacking sequence for bilayer graphene with the hoppings
represented by $\gamma_{0}$ and $\gamma_{1}$. The B sites of the top layer
(B1) are aligned with the A sites of the bottom layer (A2). (_Right_) The
reciprocal lattice with the reciprocal lattice vectors, and high-symmetry
points $\mathrm{\Gamma},\mathrm{K},\mathrm{K^{\prime}}$, and $\mathrm{M}$ are
indicated with black dots.

The low-energy band structure of the monolayer graphene is characterized by a
two-dimensional massless Dirac equation with linear dispersion at two
inequivalent hexagonal corners (often referred to as the valley) in the
Brillouin zone, namely the $\mathrm{K}$ and $\mathrm{K}^{\prime}$ points. For
a quantitative understanding of the low-energy band structure of bilayer
graphene, we introduce the continuum model Hamiltonian that describes the low-
energy band structure near the valley ($\mathrm{K}$ or $\mathrm{K}^{\prime}$)
with the basis of four atomic sites in the unit cell (A1, B1, A2, B2):

$\displaystyle
H=\begin{pmatrix}\varepsilon_{\mathrm{A1}}&v\pi^{\dagger}&0&0\\\
v\pi&\varepsilon_{\mathrm{B1}}&\gamma_{1}&0\\\
0&\gamma_{1}&\varepsilon_{\mathrm{A2}}&v\pi^{\dagger}\\\
0&0&v\pi&\varepsilon_{\mathrm{B2}}\end{pmatrix},$ (2.3)

where $\pi=\xi k_{x}+ik_{y}$ is the momentum measured from the $\mathrm{K}$
($\xi=1$) or $\mathrm{K}^{\prime}$ ($\xi=-1$) points, and
$v=\sqrt{3}a\gamma_{0}/2~{}(\approx 3\times 10^{-3}$ for large $k$ within the
continuum model) is the effective velocity, and $\gamma_{0}$ and $\gamma_{1}$
are the intralayer and interlayer coupling between the nearest neighbors,
respectively, as described in figure 1. The diagonal terms
$\varepsilon_{\mathrm{A1}}$, $\varepsilon_{\mathrm{B1}}$,
$\varepsilon_{\mathrm{A2}}$, $\varepsilon_{\mathrm{B2}}$ are the on-site
energies of the four atomic sites in the unit cell. Bilayer graphene consists
of two stacked monolayer graphene. Thus, at the $\mathrm{K}$ and
$\mathrm{K}^{\prime}$ points, there exist two copies of the Dirac Hamiltonian,
originating from the top and bottom monolayer graphene sheets, which
correspond to the upper-left and lower-right 2$\times$2 blocks of $H$ in eq.
(2.3). Due to the interlayer coupling $\gamma_{1}$, these two Dirac
Hamiltonians are coupled. This results in a pair of bands with quadratic
dispersion at the Fermi energy and another pair of bands that is split away of
the order of $\gamma_{1}$ from the Fermi energy, as shown in figure 2. For the
hopping parameters, we used the values $\gamma_{0}=3.16\,\mathrm{eV}$ and
$\gamma_{1}=0.381\,\mathrm{eV}$ [55]. More distant hoppings have a minor
effect on the electronic structure, and thus exerting only a minor influence
on the physics of interest in this study. Hence, we only consider the two
dominant hoppings, for simplicity. As $k$ moves away from the $\mathrm{K}$
($\mathrm{K}^{\prime}$) point ($k\gtrsim 100\,\mathrm{eV}$), the impact of
$\gamma_{1}$ on the band structure diminishes. The bands start to show a
linear dispersion with a band velocity $v$, similar to the monolayer case. We
emphasize that the eigenenergies of the continuum model Hamiltonian in eq.
(2.3) are isotropic with respect to the direction of the momentum
$\boldsymbol{k}$. It is also important to note that the continuum model is
valid roughly within the region where $k\lesssim 1\,\mathrm{keV}$. Beyond this
limit, the system becomes anisotropic, and a tight-binding model, which is
valid across the entire Brillouin zone, is required to describe the accurate
behavior of the band structure, rather than a continuum model [68, 69].

Figure 2: Band structure of bilayer graphene near the $\mathrm{K}$ or
$\mathrm{K}^{\prime}$ points: without any external electric field (_left_) and
with an external electric field (_right_) inducing on-site energy differences
of $U=100\,\mathrm{meV}$. The dashed line represents the Fermi energy. The
band structure near the Fermi energy is magnified on the right. Here, $|U|$
indicates the energy difference between the two bands labeled by c1 and v1 at
the $\mathrm{K}$ or $\mathrm{K}^{\prime}$ point, and $E_{g}$ represents the
band gap.

Intrinsic bilayer graphene does not have an energy gap like monolayer
graphene. However, an external potential applied to bilayer graphene induces a
gap [48, 57, 58]. Experimental studies have shown that, through doping or
electric gating, a gap of up to approximately $300\,\mathrm{meV}$ can be
achieved [59, 60]. This effect can be directly modeled from the Hamiltonian in
eq. (2.3). In the presence of an external electric field perpendicular to the
graphene ($\boldsymbol{E}_{\mathrm{ext}}\parallel\boldsymbol{z}$), the same
potential is applied within a single layer, but the two layers experience
different potentials. Denoting the potential difference between the two layers
as $U$, the on-site energies become
$\varepsilon_{\mathrm{A1}}=\varepsilon_{\mathrm{B1}}=U/2$ and
$\varepsilon_{\mathrm{A2}}=\varepsilon_{\mathrm{B2}}=-U/2$. The band
structures of BLG without and with an external potential are depicted in
figure 2. The electric field opens a gap in the low-energy band structure of
bilayer graphene, creating a band with a Mexican hat structure, where the band
gap is given by $E_{g}=|U|\gamma_{1}/\sqrt{\gamma_{1}^{2}+U^{2}}$. For small
values of $U$, the difference between $U$ and $E_{g}$ is small. For example,
$E_{g}=49.6$ meV when $U=50$ meV, and $E_{g}=96.7$ meV when $U=100$ meV.
However, we will not distinguish between the values of $U$ and $E_{g}$ in the
rest of the paper and often use both interchangeably.

##  3 Dark matter-electron scattering in bilayer graphene

In this section, we describe the methods to calculate the DM scattering and
absorption rates using the energy loss function of BLG.

###  3.1 Kinematics of dark matter

In any detector, DM scattering rate can be written as [70, 71],

$\displaystyle R(t)=\frac{1}{\rho_{T}}\frac{\rho_{\chi}}{m_{\chi}}\int
d\omega\int\frac{d^{3}q}{(2\pi)^{3}}\frac{\pi\sigma(q)}{\mu_{\chi
e}^{2}}g(\boldsymbol{q},\omega,t)S(\boldsymbol{q},\omega).$ (3.1)

Here, $\rho_{T}$ is the detector material density, $\rho_{\chi}$ is the local
DM density that we take to be $0.4\,\mathrm{GeVcm^{-3}}$ [72], $m_{\chi}$ is
the DM mass, $\sigma(q)$ is the DM-electron scattering cross section, and
$\mu_{\chi e}$ is the DM-electron reduced mass. The dynamic structure factor
$S(\boldsymbol{q},\omega)$ characterizes the response of the detector material
for a momentum $\boldsymbol{q}$ and energy deposition $\omega$.

The scattering cross section can be split into a constant factor $\sigma_{\chi
e}$ and a mediator form factor $F_{\mathrm{med}}(q)$ which is momentum-
dependent, such that $\sigma(q)=\sigma_{\chi e}F_{\mathrm{med}}(q)^{2}$, where

$\sigma_{\chi e}\equiv\dfrac{\mu_{\chi
e}^{2}}{\pi}\left(\dfrac{g_{e}g_{\chi}}{q_{0}^{2}+m_{\phi}^{2}}\right)^{2}\,,\qquad
F_{\mathrm{med}}(q)=\dfrac{q_{\mathrm{0}}^{2}+m_{\phi}^{2}}{q^{2}+m_{\phi}^{2}}\,.$
(3.2)

Here, $q_{0}=\alpha m_{e}$ is a reference momentum, $m_{\phi}$ is the mediator
mass, and $g_{e},g_{\chi}$ are the couplings of the mediator with electron and
DM, respectively. As mentioned in section 2, we will compute the BLG response
function assuming the isotropic approximation which is not valid for large $q$
and $\omega$. Hence in this paper, we will always assume a massless mediator,
i.e. $q\gg m_{\phi}$.

The function $g(\boldsymbol{q},\omega,t)$ is the only time-varying quantity in
eq.(3.1) and contains the directional information of the scattering rate. It
is obtained by integrating the DM velocity distribution, which is assumed to
be Maxwellian with a mean velocity $v_{0}=230\,\mathrm{km\,s^{-1}}$ and cut
off at the local galactic escape velocity $v_{\rm esc}$. A minimum velocity
$v_{-}$ is taken for a certain $q$ and $\omega$ [70, 71],

$g(\boldsymbol{q},\omega,t)=\dfrac{2\pi^{2}v_{0}^{2}}{qN_{0}}\left[\exp\left(-\dfrac{v_{-}(\boldsymbol{q},t)^{2}}{v_{0}^{2}}\right)-\exp\left(-\dfrac{v_{\rm
esc}^{2}}{v_{0}^{2}}\right)\right]\,,\quad
v_{-}(\boldsymbol{q},t)=\text{min}\\{v_{\mathrm{esc}},\frac{\omega}{q}+\frac{q}{2m_{\chi}}+\hat{q}\cdot\boldsymbol{v}_{\mathrm{lab}}(t)\\}\,.$
(3.3)

Here $N_{0}$ is the normalization constant of the truncated Maxwell
distribution

$\displaystyle
N_{0}=\pi^{3/2}v_{0}^{2}\left[\mathrm{erf}\left(\frac{v_{\mathrm{esc}}}{v_{0}}\right)-\frac{2v_{\mathrm{esc}}}{\pi^{1/2}}\exp\left(-\frac{v_{\mathrm{esc}}^{2}}{v_{0}^{2}}\right)\right]\,,$
(3.4)

and $\boldsymbol{v}_{\mathrm{lab}}(t)$ is the time-dependent velocity of the
laboratory on the Earth in the galactic frame,

$\boldsymbol{v}_{\mathrm{lab}}(t)=|\boldsymbol{v}_{E}|\begin{pmatrix}\sin\theta_{e}\sin\phi\\\
\sin\theta_{e}\cos\theta_{e}(\cos\phi-1)\\\
\cos^{2}\theta_{e}+\sin^{2}\theta_{e}\cos\phi\end{pmatrix}\,,$ (3.5)

with $\phi=2\pi\times(t/24\,\mathrm{h}),\theta_{e}=42^{\circ}$, and
$|\boldsymbol{v}_{E}|=240\,\mathrm{km\,s^{-1}}$. All directional information
of the DM flux is included in $v_{-}(\boldsymbol{q},t)$.

The dynamic structure factor $S(\boldsymbol{q},\omega)$ can be computed in
terms of the dielectric function $\varepsilon(\boldsymbol{q},\omega)$ using
the fluctuation-dissipation theorem as long as the DM-electron interaction is
weak [73, 74]. We describe this method in the next section.

###  3.2 Calculation of the dynamic structure factor

Equation (3.1) above needs to be modified appropriately as we are going to use
BLG in this study, which is a 2D material. The dynamic structure factor
encodes the information of the correlation of the density fluctuations in a
system at a momentum $\boldsymbol{q}$ and energy $\omega$. It can be expressed
using the electron density operator as the Fourier transform of the density-
density correlation function [75, 76], i.e.

$\displaystyle S(\boldsymbol{q},\omega)=\frac{1}{V}\int dte^{i\omega
t}\langle\hat{n}(\boldsymbol{q},t)\hat{n}(-\boldsymbol{q},0)\rangle\,.$ (3.6)

Here, $\hat{n}(\boldsymbol{q},t)$ is the electron density operator at a
momentum $\boldsymbol{q}$ and time $t$:

$\displaystyle\hat{n}(\boldsymbol{q},t)=\int
d^{3}r\,e^{-i\boldsymbol{q}\cdot\boldsymbol{r}}\hat{n}(\boldsymbol{r},t)\,,$
(3.7)

where $\hat{n}(\boldsymbol{r},t)$ is the electron density operator in real
space. For 2D materials like BLG, the electron density is confined to the 2D
space (near $r_{z}=0$). We can assume that $\hat{n}(\boldsymbol{r},t)$ can be
separated into the material’s in-plane ($r_{x}$, $r_{y}$) and out-of-plane
($r_{z}$) components. Then the electron density operator can be written as

$\displaystyle\hat{n}(\boldsymbol{r},t)=\hat{n}_{\mathrm{2D}}(\boldsymbol{r}_{\parallel},t)\left[\frac{1}{d}\Theta(d/2-|r_{z}|)\right].$
(3.8)

Here, we make the assumption that the electron density remains constant
between the two layers of BLG, defined within the range $-d/2<r_{z}<d/2$,
where $d=3.35$ $\mathrm{\AA}$ represents the interlayer distance. The spatial
coordinates in the in-plane direction are denoted as
$\boldsymbol{r}_{\parallel}=(r_{x},r_{y})$, and $\Theta(x)$ represents the
Heaviside step function.

Note that we have assumed a step function-like profile for the real-space
electron density, as such an assumption does not significantly impact the
results. Given that the DM mass is on the sub-MeV regime, the order of
magnitude for $q$ is roughly in the sub-keV range, implying $qd\ll 1$. This
implies that the length scale of the interlayer separation in bilayer graphene
is much smaller than the inverse of the typical momentum transfer $\sim 1/q$,
rendering the detailed structure of the real-space density profile less
crucial for our considerations. Plugging eq. (3.8) into eq. (3.7), we get

$\displaystyle\hat{n}(\boldsymbol{q},t)$ $\displaystyle=\int
d^{3}re^{-i\boldsymbol{q}\cdot\boldsymbol{r}}\hat{n}(\boldsymbol{r},t)$
$\displaystyle=\left[\int
d^{2}r_{\parallel}e^{-i\boldsymbol{q}_{\parallel}\cdot\boldsymbol{r}_{\parallel}}n_{\mathrm{2D}}(\boldsymbol{r}_{\parallel},t)\right]\left[\int
dr_{z}e^{-iq_{z}r_{z}}\frac{1}{d}\Theta(d/2-|r_{z}|)\right]$
$\displaystyle=\hat{n}_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},t)\left(\frac{\sin{q_{z}d/2}}{q_{z}d/2}\right).$
(3.9)

Using this, we can rewrite eq. (3.6) in terms of the electron density operator
in 2D as

$\displaystyle S(\boldsymbol{q},\omega)$ $\displaystyle=\frac{1}{V}\int
dte^{i\omega
t}\langle\hat{n}(\boldsymbol{q},t)\hat{n}(-\boldsymbol{q},0)\rangle$
$\displaystyle=\frac{1}{d}\left(\frac{1}{A}\int dte^{i\omega t}\langle
n_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},t)n_{\mathrm{2D}}(-\boldsymbol{q}_{\parallel},0)\rangle\right)\left(\frac{\sin{q_{z}d/2}}{q_{z}d/2}\right)^{2}$
$\displaystyle=\frac{1}{d}S_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)\left(\frac{\sin{q_{z}d/2}}{q_{z}d/2}\right)^{2},$
(3.10)

where $V=Ad$ is the 3D volume, with $A$ being the 2D volume or surface area.
Here, $S_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)$ denotes the dynamic
structure factor of the 2D material, and the following square of the sinc
function contains the variation in the $q_{z}$-direction. Combining the
aforementioned information, the detector component
$\frac{1}{\rho_{T}}S(\boldsymbol{q},\omega)$ in the DM-electron scattering
rate in eq. (3.1) transforms as

$\displaystyle\frac{1}{\rho_{T}}S(\boldsymbol{q},\omega)\ \rightarrow\
\frac{1}{\sigma_{T}}S_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)\left(\frac{\sin{q_{z}d/2}}{q_{z}d/2}\right)^{2},$
(3.11)

where $\sigma_{T}=1.53\times 10^{-7}\mathrm{g/cm^{2}}$ is the surface density
of BLG. From here on, unless otherwise specified, the dynamic structure factor
will refer to $S_{\mathrm{2D}}$. Specifically, the dynamic structure factor
can be written as the imaginary (dissipative) part of the 2D density-density
response function $\chi_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)$
derived from the fluctuation-dissipation theorem [75, 76]:

$\displaystyle
S_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)=-\frac{2}{1-e^{-\beta\omega}}\mathrm{Im}\left[\chi_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)\right],$
(3.12)

where $\beta=1/(k_{B}T)$. The 2D density-density response function in the
frequency-momentum space is the Fourier transform of the 2D retarded density
correlation function given by

$\displaystyle\chi_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)=-\frac{i}{A}\int
dte^{i\omega
t}\left<\left[\hat{n}_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},t),\hat{n}_{\mathrm{2D}}(-\boldsymbol{q}_{\parallel},0)\right]\right>.$
(3.13)

Note that the expectation value $\langle\cdots\rangle$ refers to the thermal
average over many-body states, taking into account the Coulomb interaction
between the electrons. It is often practical to introduce the dielectric
function of the detector. The dielectric function, in the frequency-momentum
space, is defined as the ratio of the screened potential to the bare potential
in a material. The 2D version of the dielectric function is given by

$\displaystyle\frac{1}{\varepsilon(\boldsymbol{q}_{\parallel},\omega)}=1+v_{\mathrm{C}}(\boldsymbol{q}_{\parallel})\chi_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)\,,$
(3.14)

where $v_{\mathrm{C}}(q)=\dfrac{e^{2}}{2\varepsilon_{r}q}$ is the 2D Coulomb
interaction. Here, $\varepsilon_{r}$ is the dielectric constant of the medium
which we assume to be a vacuum. At zero temperature, the dynamic structure
factor in eq. (3.12) becomes

$\displaystyle
S_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)=-2\Theta(\omega)\mathrm{Im}\chi_{\mathrm{2D}}(\omega)=\frac{2\Theta(\omega)}{v_{\mathrm{C}}(\boldsymbol{q}_{\parallel})}\mathcal{W}(\boldsymbol{q}_{\parallel},\omega),\quad\mathcal{W}(\boldsymbol{q}_{\parallel},\omega)=\mathrm{Im}\left[-\frac{1}{\varepsilon(\boldsymbol{q}_{\parallel},\omega)}\right].$
(3.15)

This implies a direct connection between the dynamic structure factor and the
imaginary part of the inverse dielectric function, termed the energy loss
function (ELF) $\mathcal{W}(\boldsymbol{q},\omega)$. The ELF represents
electronic energy loss in materials and is thus related to an experimentally
measurable quantity. The information about the response to a small density
perturbation is contained in the dielectric function.

To calculate the ELF, we can use the random phase approximation (RPA) to
compute the expectation value in eq. (3.13). Diagrammatically, the RPA
corresponds to an infinite series of the non-interacting polarization function
connected by the Coulomb interaction lines, neglecting other higher-order
terms [77, 78, 75, 76]:

$\displaystyle\chi_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)=\frac{\chi_{\mathrm{2D}}^{(0)}(\boldsymbol{q}_{\parallel},\omega)}{\varepsilon(\boldsymbol{q}_{\parallel},\omega)}\,,\quad\varepsilon(\boldsymbol{q}_{\parallel},\omega)=1-v_{\mathrm{C}}(\boldsymbol{q}_{\parallel})\chi_{\mathrm{2D}}^{(0)}(\boldsymbol{q}_{\parallel},\omega)\,,$
(3.16)

where the response function in eq. (3.13) and the dielectric function in eq.
(3.14) are expressed in terms of the non-interacting polarization function
$\chi^{(0)}_{\mathrm{2D}}(\boldsymbol{q},\omega)$ using the RPA. This approach
is known to be valid in the weak coupling (typically high density) limit or in
the large fermion flavor limit.
$\chi^{(0)}_{\mathrm{2D}}(\boldsymbol{q},\omega)$ is given by [75]

$\displaystyle\chi^{(0)}_{\mathrm{2D}}(\boldsymbol{q}_{\parallel},\omega)=g\sum_{n,m}\int\frac{d^{2}k}{(2\pi)^{2}}\frac{f_{\boldsymbol{k},n}-f_{\boldsymbol{k}+\boldsymbol{q}_{\parallel},m}}{\omega+\varepsilon_{\boldsymbol{k},n}-\varepsilon_{\boldsymbol{k}+\boldsymbol{q}_{\parallel},m}+i\eta}|\langle\boldsymbol{k},n|\boldsymbol{k}+\boldsymbol{q}_{\parallel},m\rangle|^{2},$
(3.17)

where $\varepsilon_{\boldsymbol{k},n}$ and
$f_{\boldsymbol{k},n}=1/[1+e^{\beta(\varepsilon_{\boldsymbol{k},n}-\mu)}]$ are
the eigenenergy and the Fermi-Dirac distribution for the wave vector
$\boldsymbol{k}$ and band index $n$, respectively, $g=4$ is the spin-valley
degeneracy factor for graphene, and $\mu$ is the chemical potential, $\eta$ is
a phenomenological broadening parameter proportional to the inverse lifetime
(or decay width) of quasiparticles which takes $\eta\rightarrow 0^{+}$ for the
clean limit. For numerical calculations, we set $\eta=1$ meV, which is the
typical order of magnitude for BLG [79, 80, 81]. We note that variations in
$\eta$ have a negligible impact on the overall results.

In this work, we will also consider DM absorption in BLG. To this end, we will
show the results for a dark photon model where the DM mixes with the Standard
Model (SM) photon through a small mixing parameter $\kappa$,

$\displaystyle\mathcal{L}\supset\kappa F_{\mu\nu}F^{\prime\mu\nu}\,,$ (3.18)

where $F^{\mu\nu}$ and $F^{\prime\mu\nu}$ are the field strengths of the SM
photon and the dark photon, respectively. Because of the mixing between dark
photon DM and the Standard Model photon, DM absorption rate is again related
to the ELF of the material rescaled by the mixing parameter $\kappa$ as
follows [82],

$\displaystyle R_{\rm
abs}=\kappa^{2}\rho_{\chi}\dfrac{d}{\sigma_{T}}\,\mathrm{Im}\left[-\frac{1}{\varepsilon(q,\omega)}\right]\,,$
(3.19)

with $\omega=m_{\chi}$ and $q=m_{\chi}v$. Because $q\ll\omega$ for DM
absorption, it is more convenient to calculate the dielectric function in 2D
using the following expression [83, 78]

$\displaystyle\varepsilon(q\rightarrow
0,\omega)=1+\dfrac{i}{2\varepsilon_{r}}\dfrac{q}{\omega}\sigma(\omega).$
(3.20)

Here $\varepsilon_{r}$ denotes the dielectric constant of the medium which we
assume to be vacuum, and the optical conductivity $\sigma(\omega)$ of BLG is
calculated using the Kubo formula [77]. We will show projection of constraint
on the mixing parameter $\kappa$.

##  4 Results

We now have all the tools needed to calculate the ELF of BLG, and the DM
scattering and absorption rates. In this section, we will show and explain our
main results.

###  4.1 Dielectric properties of bilayer graphene

Figure 3: The color plot of the energy loss function (ELF) for BLG with
$U=50\,\mathrm{meV}$ (_left_) and $100\,\mathrm{meV}$ (_right_), represented
in the $q$ \- $\omega$ plane. The black dashed line indicates the boundary of
the electron-hole continuum, with the more intensely colored region above
signifying the continuum itself. The kinematically-allowed boundaries for DM
with mass $m_{\chi}=20,100$, and $1000\,\mathrm{keV}$ are depicted by the red
solid lines. The inset in the right panel figure shows the different
electronic transitions that give rise to the features of the ELF. See section
4.1 for more details.

The ELF of BLG in the presence of a potential, calculated through eqs.
(3.15)-(3.17), is shown as a color plot in the $q$ \- $\omega$ plane in figure
3. Since the low-energy band structure described by the continuum model
Hamiltonian in eq. (2.3) is isotropic, the ELF also depends only on the
magnitude of $\boldsymbol{q}$ and not on its direction.

The magnitude of the ELF is related to the electronic energy loss at the
momentum transfer $q$ and energy transfer $\omega$. Especially when there is a
gap and the Fermi energy lies within it, the origin of the ELF is the
interband transitions of electrons. This represents cases where momentum
transfer by $q$ and energy transfer by $\omega$ are possible from an occupied
band to an empty band, and this region is referred to as the electron-hole
continuum [75]. For example, in the region where $\omega<E_{g}$, the ELF is
close to zero, which means the absence of available states for excitation
below the gap. The boundary of the electron-hole continuum is indicated by a
dashed line in figure 3. Additionally, for $q\gtrsim 100\,\mathrm{keV}$, the
boundary of the electron-hole continuum in $q$ \- $\omega$ plane shows a
linear dispersion with the velocity $v$ in eq. (2.3). We note that since the
Fermi energy lies within the gap, implying no states exist at that energy,
contributions from intraband transitions or plasmons to the ELF are absent in
this case.

The red curves indicate the kinematic phase space of DM of mass $m_{\chi}$.
Here, the mean DM velocity is assumed to be $v_{0}=230\,\mathrm{km\,s^{-1}}$.
In the phase space where there is an overlap between this kinematically-
allowed region and the dark areas of the BLG’s ELF, DM-electron scattering
occurs. If there is no overlap, the scattering rate vanishes. For example, in
the case where $m_{\chi}=20\,\mathrm{keV}$, the kinematically-allowed region
exists outside the electron-hole continuum, leading to no DM-electron
scattering. This means that the $U$ applied to the BLG determines the minimum
detectable mass threshold for DM. The continuum model we employed is valid for
ranges approximately $q\lesssim 1\,\mathrm{keV}$ and $\omega\lesssim
2\,\mathrm{eV}$. For DM with a mass $m_{\chi}=1\,\mathrm{MeV}$, the
kinematically-allowed boundary extends into regions with $q>1\,\mathrm{keV}$.
However, the region contributing to the DM-electron scattering rate, that is,
the area overlapping with the electron-hole continuum, is approximately
$q\lesssim 0.6\,\mathrm{keV}$. Therefore, the description using the continuum
model Hamiltonian remains valid. Consequently, our results are relevant for DM
mass range up to approximately $m_{\chi}\lesssim 1\,\mathrm{MeV}$.

For larger values of $m_{\chi}$, a tight-binding model that is valid across
the entire Brillouin zone is required, rather than a continuum model. Then the
band structure no longer remains isotropic in the momentum space, leading to
the loss function developing a directional dependence on $\boldsymbol{q}$.
Additionally, when the magnitude of the $q$ exceeds the size of the Brillouin
zone, considerations like the local field effects [75, 84] become necessary in
the determination of the dielectric function. On the low mass end, we are
limited by the finite value of $U$. Energy deposition below $U$ is not
possible which will help reduce the thermal fluctuations.

Figure 4: Energy spectrum of DM scattering rate in bilayer graphene with two
band gaps $U=50\,\mathrm{meV}$ (_left_) and $100\,\mathrm{meV}$ (_right_) for
$m_{\chi}=100$ (violet), 200 (pink), and $300\,\mathrm{keV}$ (sky blue). A
massless mediator and $\sigma_{\chi e}=10^{-38}\,\mathrm{cm^{2}}$ are assumed
in all cases. The energy thresholds are marked with a vertical gray line. The
gray-shaded below-threshold regions are not taken into account in our results.

Having discussed the properties of the ELF of bilayer graphene, we now compute
the DM scattering rate following eq. (3.1). We show the scattering rate as a
function of energy in figure 4 for $m_{\chi}=100,200,$ and $300\,\mathrm{keV}$
for two values of $U=50,100\,\mathrm{meV}$ and cross section $\sigma_{\chi
e}=10^{-38}\,\mathrm{cm^{2}}$. These thresholds are marked with a vertical
gray line. As mentioned before, we will only show results for the massless
mediator case. From figure 4, we can see that the rates are maximum near the
thresholds and decreases toward higher energies. This maximum rate comes from
the electronic transitions from the valence to the conduction bands $v_{1}\to
c_{1}$ as shown in the inset in the right panel of figure 3 (green arrows).
Moreover, the overall rate increases for heavier DM. This is understandable
from the overlap between the DM kinematic phase space and the ELF of BLG in
the $q$ \- $\omega$ plane in figure 3. For low mass DM, the overlap region is
small and limited toward the low energy corner of the plane, hence the
scattering rate is maximum around $E_{g}$. For heavier mass, the overlap
region is extended to higher energies and has support from the additional
$v_{1}\to c_{2}$ and $v_{2}\to c_{1}$ electronic transitions marked with the
orange arrows in the inset of figure 3. This broadens the spectrum toward
higher energy for heavier mass DM as can be seen in figure 4. Finally, we note
that the rates are higher for the case $E_{g}=50\,\mathrm{meV}$ because the
ELF in this case peaks at a relatively lower energy ($v_{1}\to c_{1}$
transition) yielding a greater overlap with the DM kinematic phase space. This
is also evident from figure 3, indicating that the sensitivity could further
improve for smaller $E_{g}$. However, a smaller value of $E_{g}$ will also be
subject to the technical challenge of maintaining the same gate voltage over a
large volume of the target. The optimal value of $E_{g}$ will thus depend on
the details of the experimental design.

Before ending this section, we want to comment on the additional below-
threshold peaks that appear in the gray-shaded region in figure 4. These are
arising from the small but nonzero value of the ELF in the $\omega<E_{g}$
energy regime which, in turn, is a result of our assumption of a nonzero decay
width $\eta$ of the excited electron-hole pairs in eq. (3.17). This small
value of the ELF, combined with the diverging nature of the form factor
$F_{\mathrm{med}}(q)\propto 1/q^{2}$ at small $q$, gives rise to these
additional peaks. However, we note that if $\eta$ were zero, the ELF would
completely vanish below $\omega<E_{g}$, resulting in a zero value of
$dR/d\omega$ in this energy range. In our calculation, we do not include the
energy range below $\omega<E_{g}=E_{\mathrm{th}}$, and assume the thresholds
to be the same as the band gap in each case.

###  4.2 Projected sensitivity & daily modulation

In this section, we will explore the future projections with bilayer graphene
as the detector to search for DM. We envisage an experimental setup similar to
the ones developed for photon detectors with monolayer and bilayer graphene
[85, 86, 87, 88, 89, 90]. Graphene being an atomic-thick material, one needs a
relatively large surface area to achieve a reasonable amount of detector mass.
One way is to prepare multiple sheets of BLG and build a single detector
module. A single experiment then can deploy multiple such modules to achieve
larger target mass. In this regard, we want to mention about an upcoming
experiment called PTOLEMY that is planning to use tritium adsorbed onto
graphene substrate to detect low energy cosmic neutrinos [91, 92]. The
technical design of a BLG DM detector can get help from PTOLEMY’s prototype
detector design and background study as their goal is also to achieve a large
exposure ($\sim$ 100 g-yr). To reach the desirable energy threshold, the BLG
needs to be biased with the appropriate voltage, and the setup needs to be
cooled down below that temperature to mitigate the background noise from
thermal fluctuations. Moreover, the excited electron-hole pairs created from
DM scattering and absorption need to be detected. This can be done with
sensitive charge detectors similar to the charge-coupled devices [93, 94]. In
this paper, we do not go into any further technical detail about the
experimental setup and leave it for future work.

Figure 5: Projected sensitivties on DM-electron scattering cross section
$\sigma_{\chi e}$ (_left_) and DM mixing parameter $\kappa$ for absorption
(_right_) in bilayer graphene with two energy thresholds $E_{\mathrm{th}}=50$
and $100\,\mathrm{meV}$ assuming 3 events in 10 mg-yr exposure with no
background. We use scattering rate averaged over a 24 hour period. The blue
band marked with ‘Freeze-in’ in the left panel is the region of theoretical
interest. For comparison, we also show the projection for superconducting
aluminum (SC-Al) as a gray, dotted line [95]. Bilayer graphene is superior to
superconducting aluminum for DM scattering in the
$20\,\mathrm{keV}-\,\mathrm{MeV}$ mass range. The dips appearing for the DM
absorption sensitivity curve in the right panel are due to opening up of
different band transitions of electrons. See the text for more details. The
gray-shaded regions are excluded by other laboratory and astrophysical
constraints [14, 96, 10, 97, 98, 99, 100].

We show the projected sensitivity for DM-electron scattering in the left panel
for the case of a massless mediator, and DM absorption in the right panel of
figure 5. As mentioned before, we show the results for two values of
$U=50\,\mathrm{meV}$ (orange) and $100\,\mathrm{meV}$ (teal) assuming a cutoff
of 3 events for an exposure of 10 mg-yr for both cases with no background. We
lose sensitivity below $m_{\chi}\simeq 20$ and $35\,\mathrm{keV}$,
respectively, for the two band gaps which can be understood from their
respective energy thresholds. For comparison, we also show the projection for
superconducting aluminum as a gray, dotted line [95]. The left panel of figure
5 shows that BLG performs better than superconducting aluminum for DM
scattering in the $10\,\mathrm{keV}-\,\mathrm{MeV}$ mass range for the same
amount of exposure. With the assumed 3 events per 10 mg-yr exposure limit, BLG
sensitivity can reach $\sigma_{\chi e}\simeq 10^{-38}\,\mathrm{cm^{2}}$ for
$m_{\chi}\simeq 60\,\mathrm{keV}$, and hence, the band of freeze-in benchmark
model marked with blue in figure 5.

In the case of DM absorption, we show our sensitivity projections on the
mixing parameter $\kappa$ of dark photon DM in the mass range
$0.01<m_{\chi}<1\,\mathrm{eV}$ as the whole DM mass energy is absorbed
($\omega=m_{\chi}$). It is computed using eq. (3.19) and the same limit 3
events per 10 mg-yr exposure with no background. We see that the sensitivity
can reach $\kappa\simeq 10^{-12}$ as shown in the right panel of figure 5. The
dips in the DM absorption sensitivity curve are coming from the opening up of
electronic transitions between different bands, such as $v_{1}\to c_{1}$,
$v_{1}\to c_{2}$, or $v_{2}\to c_{1}$, in the ELF as was discussed in section
4.1. The gray-shaded regions are excluded by existing laboratory and
astrophysical constraints [14, 96, 10, 97, 98, 99, 100]. We computed the
sensitivity of BLG for pseudoscalar DM absorption too which is also related to
the conductivity [82]. However, we do not show it here as it does not reach a
coupling strength that is not already constrained by other astrophysical
observations.

Another important characteristic of a BLG detector is the daily modulation of
the DM scattering rate. This arises from the change in the direction of the
incoming DM wind due to the rotation of the Earth about its axis. The two-
dimensional nature of BLG facilitates the daily modulation even more as it
makes the detector response $S(\boldsymbol{q},\omega)$ inherently anisotropic.
The time modulation information in the scattering rate is included in the
$g(\boldsymbol{q},\omega,t)$ factor in eq.(3.1). In figure 6, we show the
daily modulation of the scattering rate normalized by the average
$R(t)/\langle R\rangle$ over a 24 hr period for a few DM masses in BLG with
$U=100\,\mathrm{meV}$. The rate varies by upto $R(t)/\langle R\rangle\simeq
80\%$ for $m_{\chi}=120\,\mathrm{keV}$ and by smaller amount for heavier DM
mass. The amount of modulation is more for lighter DM mass. This is
understandable as the typical energy transfer is closer to the threshold for
lighter DM and a small change in the relative velocity of DM can affect the
scattering rate significantly. Such daily modulation of DM scattering rate can
be leveraged to isolate any potential DM signal in a future experiment from
various background scattering events from cosmic rays, radioactivity, or any
internal systematic noise that do not share such time variation. We note that
on top of the daily modulation, there is also an annual modulation expected
from the rotation of the Earth around the sun resulting in a slightly
different relative velocities of the DM wind throughout the year. We do not
calculate it in this work.

Figure 6: Daily modulation of the scattering rate of a few DM masses in
bilayer graphene with $E_{g}=50$ (_left_) and $100\,\mathrm{meV}$ (_right_).
The modulation is more pronounced for lighter DM as the typical energy
transfer is closer to the threshold, and it is milder for lower threshold
$50\,\mathrm{meV}$ compared to $100\,\mathrm{meV}$ due to the same reason.

In passing, we want to point out that the DM in the Milky Way halo may have a
kinematic substructure in the solar neighborhood with a different density and
velocity distribution. In such a case, the sensitivity projections and the
daily modulation may get affected because of the substructure [101, 102].

##  5 Conclusions & Outlook

In this paper, we proposed bilayer graphene (BLG) as a target material for DM
direct detection experiment and studied its sensitivity to probe DM-ordinary
matter interaction. BLG is a 2D material consisting of two honeycomb lattices
of carbon atoms stacked together. The electrons in a BLG possess a tunable
band gap that can be adjusted with a gate voltage applied across the material.
This band gap can be used as the energy threshold in the experiment giving us
a convenient control on the thermal noise background. We computed the
dielectric function of BLG using the random phase approximation and used it to
compute the DM scattering and absorption rates. We then showed projections for
limits on the DM-electron scattering cross section $\sigma_{\chi e}$ for a
massless mediator, and on the kinetic mixing parameter $\kappa$ between dark
photon DM and ordinary photon. Our results show that a 10 mg-yr experimental
exposure with no background will allow us to probe $\sigma_{\chi e}\gtrsim
10^{-38}\,\mathrm{cm^{2}}$ around $m_{\chi}\simeq 50\,\mathrm{keV}$, and reach
the freeze-in benchmark models in the DM-electron scattering parameter space.
Thanks to the small value of the band gap, BLG can detect $\sim
20\,\mathrm{keV}$ mass DM, and potentially even reach below $10\,\mathrm{keV}$
if lower energy threshold is achievable. Such small energy threshold may also
allow us to probe the very low energy Earth-bound thermalized DM population
predicted by some models. In the context of a dark photon model, DM absorption
in BLG can probe $\kappa\gtrsim 10^{-12}$ in the sub-eV mass range. We also
compute the expected daily modulation in the DM scattering rate due to the
rotation of the Earth about its axis. We found that the scattering rate can
vary by upto $\sim 80\%$ from the average value for lighter DM. This
modulation can help us distinguish the DM signal from other constant
background sources.

The promising results from this study motivates us to ask next - how to build
a DM detector using BLG? A detailed and more dedicated study is required to
answer this question. We plan to pursue that in a future work. In this regard,
we pointed out the future experiment called PTOLEMY that is planning to detect
low energy cosmic neutrinos using tritium adsorbed on graphene layers.
Perhaps, a future experiment using BLG can get technical help from PTOLEMY as
their aim is to achieve an ambitious 100 g-yr exposure [91, 103, 92]. In the
projection calculation of this paper, we did not assume any background.
However, in reality one should expect various sources of backgrounds, such as
cosmic rays, radioactive sources, thermal noise etc. The muons and gamma rays
from the cosmic rays can be mitigated by doing the experiment in an
underground laboratory. The environmental radioactive sources are more
difficult to get rid of. Appropriate shielding of the experimental setup and
careful modeling of the remaining radioactive sources is one solution.
Additionally, there could be hardware-specific systematics that one might need
to pay attention to [35].

Both, direct detection of sub-MeV DM and graphene are relatively new areas of
research with many novel ideas being developed today. Hence, there are plenty
of scopes to improve the idea presented in this paper. For example, figure 3
clearly tells us that the sensitivity of bilayer graphene as DM target can be
further improved by matching its band velocity $v$ with the DM velocity, or by
extending the overlap between the detector response and DM phase space using
flatter bands. Twisted bilayer [63, 64, 65, 66] or multilayer ($>2$) graphene
[104, 105] could be useful for this purpose. Recent research has seen
significant progress in studying flat band materials, where the bands near the
Fermi level exhibit flat characteristics [106]. Various materials, such as the
kagome lattice, have been discovered [107, 108]. Such flat band materials can
offer improvement in the sensitivity of DM detectors. We want to pursue these
possibilities in future.

## Acknowledgment

We thank Yoni Kahn and Noah Kurinsky for helpful discussions during the work.
AD was supported by Grant Korea NRF-2019R1C1C1010050. J.J. and H.M.
acknowledge support from the National Research Foundation of Korea (NRF)
grants funded by the Korea government (MSIT) (Grants No. 2023R1A2C1005996) and
the Creative-Pioneering Researchers Program through Seoul National University
(SNU).

## References

  * [1] Y. Hochberg, E. Kuflik, T. Volansky and J.G. Wacker, _Mechanism for Thermal Relic Dark Matter of Strongly Interacting Massive Particles_ , _Phys. Rev. Lett._ 113 (2014) 171301 [1402.5143].
  * [2] C. Boehm, P. Fayet and J. Silk, _Light and heavy dark matter particles_ , _Phys. Rev. D_ 69 (2004) 101302 [hep-ph/0311143].
  * [3] C. Boehm and P. Fayet, _Scalar dark matter candidates_ , _Nucl. Phys. B_ 683 (2004) 219 [hep-ph/0305261].
  * [4] P. Fayet, _Light spin 1/2 or spin 0 dark matter particles_ , _Phys. Rev. D_ 70 (2004) 023514 [hep-ph/0403226].
  * [5] E. Izaguirre, G. Krnjaic, P. Schuster and N. Toro, _Analyzing the Discovery Potential for Light Dark Matter_ , _Phys. Rev. Lett._ 115 (2015) 251301 [1505.00011].
  * [6] J.H. Chang, R. Essig and A. Reinert, _Light(ly)-coupled Dark Matter in the keV Range: Freeze-In and Constraints_ , _JHEP_ 03 (2021) 141 [1911.03389].
  * [7] R.T. D’Agnolo, C. Mondino, J.T. Ruderman and P.-J. Wang, _Exponentially Light Dark Matter from Coannihilation_ , _JHEP_ 08 (2018) 079 [1803.02901].
  * [8] CRESST collaboration, _Results on MeV-scale dark matter from a gram-scale cryogenic calorimeter operated above ground_ , _Eur. Phys. J. C_ 77 (2017) 637 [1707.06749].
  * [9] CRESST collaboration, _First results from the CRESST-III low-mass dark matter program_ , _Phys. Rev. D_ 100 (2019) 102002 [1904.00498].
  * [10] SuperCDMS collaboration, _Light Dark Matter Search with a High-Resolution Athermal Phonon Detector Operated Above Ground_ , _Phys. Rev. Lett._ 127 (2021) 061801 [2007.14289].
  * [11] R. Ren et al., _Design and characterization of a phonon-mediated cryogenic particle detector with an eV-scale threshold and 100 keV-scale dynamic range_ , _Phys. Rev. D_ 104 (2021) 032010 [2012.12430].
  * [12] CPD collaboration, _Performance of a large area photon detector for rare event search applications_ , _Appl. Phys. Lett._ 118 (2021) 022601 [2009.14302].
  * [13] “The TESSERACT dark matter project, SNOWMASS LOI.” www.snowmass21.org/docs/files/summaries/CF/SNOWMASS21-CF1_CF2-IF1_IF8-120.pdf.
  * [14] SENSEI collaboration, _SENSEI: Direct-Detection Results on sub-GeV Dark Matter from a New Skipper-CCD_ , _Phys. Rev. Lett._ 125 (2020) 171802 [2004.11378].
  * [15] DAMIC-M collaboration, _First Constraints from DAMIC-M on Sub-GeV Dark-Matter Particles Interacting with Electrons_ , _Phys. Rev. Lett._ 130 (2023) 171003 [2302.02372].
  * [16] Oscura collaboration, _The Oscura Experiment_ , 2202.10518.
  * [17] S. Rajendran, N. Zobrist, A.O. Sushkov, R. Walsworth and M. Lukin, _A method for directional detection of dark matter using spectroscopy of crystal defects_ , _Phys. Rev. D_ 96 (2017) 035009 [1705.09760].
  * [18] M.C. Marshall, M.J. Turner, M.J.H. Ku, D.F. Phillips and R.L. Walsworth, _Directional detection of dark matter with diamond_ , _Quantum Sci. Technol._ 6 (2021) 024011 [2009.01028].
  * [19] S.A. Hertel, A. Biekert, J. Lin, V. Velan and D.N. McKinsey, _Direct detection of sub-GeV dark matter using a superfluid 4He target_, _Phys. Rev. D_ 100 (2019) 092007 [1810.06283].
  * [20] H.J. Maris, G.M. Seidel and D. Stein, _Dark Matter Detection Using Helium Evaporation and Field Ionization_ , _Phys. Rev. Lett._ 119 (2017) 181303 [1706.00117].
  * [21] R. Anthony-Petersen et al., _Applying Superfluid Helium to Light Dark Matter Searches: Demonstration of the HeRALD Detector Concept_ , 2307.11877.
  * [22] G.-B. Kim, _Self-absorption and Phonon Pulse Shape Discrimination in Scintillating Bolometers_ , _J. Low Temp. Phys._ 199 (2020) 1004.
  * [23] Y. Hochberg, Y. Zhao and K.M. Zurek, _Superconducting Detectors for Superlight Dark Matter_ , _Phys. Rev. Lett._ 116 (2016) 011301 [1504.07237].
  * [24] Y. Hochberg, I. Charaev, S.-W. Nam, V. Verma, M. Colangelo and K.K. Berggren, _Detecting Sub-GeV Dark Matter with Superconducting Nanowires_ , _Phys. Rev. Lett._ 123 (2019) 151802 [1903.05101].
  * [25] J. Chiles et al., _New Constraints on Dark Photon Dark Matter with Superconducting Nanowire Detectors in an Optical Haloscope_ , _Phys. Rev. Lett._ 128 (2022) 231802 [2110.01582].
  * [26] Y.-H. Kim, S.-J. Lee and B. Yang, _Superconducting detectors for rare event searches in experimental astroparticle physics_ , _Superconductor Science and Technology_ 35 (2022) 063001.
  * [27] A. Das, N. Kurinsky and R.K. Leane, _Dark Matter Induced Power in Quantum Devices_ , 2210.09313.
  * [28] C.W. Fink, C.P. Salemi, B.A. Young, D.I. Schuster and N.A. Kurinsky, _The Superconducting Quasiparticle-Amplifying Transmon: A Qubit-Based Sensor for meV Scale Phonons and Single THz Photons_ , 2310.01345.
  * [29] Z.-L. Liang, C. Mo, F. Zheng and P. Zhang, _Phonon-mediated Migdal effect in semiconductor detectors_ , _Phys. Rev. D_ 106 (2022) 043004 [2205.03395].
  * [30] Y. Hochberg, Y. Kahn, M. Lisanti, K.M. Zurek, A.G. Grushin, R. Ilan et al., _Detection of sub-MeV Dark Matter with Three-Dimensional Dirac Materials_ , _Phys. Rev. D_ 97 (2018) 015004 [1708.08929].
  * [31] R.M. Geilhufe, F. Kahlhoefer and M.W. Winkler, _Dirac Materials for Sub-MeV Dark Matter Detection: New Targets and Improved Formalism_ , _Phys. Rev. D_ 101 (2020) 055005 [1910.02091].
  * [32] Z. Huang, C. Lane, S.E. Grefe, S. Nandy, B. Fauseweh, S. Paschen et al., _Dark Matter Detection with Strongly Correlated Topological Materials: Flatband Effect_ , 2305.19967.
  * [33] C. Blanco, R. Essig, M. Fernandez-Serra, H. Ramani and O. Slone, _Dark matter direct detection with quantum dots_ , _Phys. Rev. D_ 107 (2023) 095035 [2208.05967].
  * [34] R. Essig, G.K. Giovanetti, N. Kurinsky, D. McKinsey, K. Ramanathan, K. Stifter et al., _Snowmass2021 Cosmic Frontier: The landscape of low-threshold dark matter direct detection in the next decade_ , in _Snowmass 2021_ , 3, 2022 [2203.08297].
  * [35] P. Adari et al., _EXCESS workshop: Descriptions of rising low-energy spectra_ , _SciPost Phys. Proc._ 9 (2022) 001 [2202.05097].
  * [36] A. Gould and G. Raffelt, _Thermal Conduction by Massive Particles_ , _Astrophys. J._ 352 (1990) 654.
  * [37] H. Banks, S. Ansari, A.C. Vincent and P. Scott, _Simulation of energy transport by dark matter scattering in stars_ , _JCAP_ 04 (2022) 002 [2111.06895].
  * [38] V. De Luca, A. Mitridate, M. Redi, J. Smirnov and A. Strumia, _Colored Dark Matter_ , _Phys. Rev. D_ 97 (2018) 115024 [1801.01135].
  * [39] B. Dasgupta, A. Gupta and A. Ray, _Dark matter capture in celestial objects: light mediators, self-interactions, and complementarity with direct detection_ , _JCAP_ 10 (2020) 023 [2006.10773].
  * [40] M. Pospelov and H. Ramani, _Earth-bound millicharge relics_ , _Phys. Rev. D_ 103 (2021) 115031 [2012.03957].
  * [41] R.K. Leane and J. Smirnov, _Floating dark matter in celestial bodies_ , _JCAP_ 10 (2023) 057 [2209.09834].
  * [42] S. Baum, L. Visinelli, K. Freese and P. Stengel, _Dark matter capture, subdominant WIMPs, and neutrino observatories_ , _Phys. Rev. D_ 95 (2017) 043007 [1611.09665].
  * [43] S. Rajendran and H. Ramani, _Composite solution to the neutron lifetime anomaly_ , _Phys. Rev. D_ 103 (2021) 035014 [2008.06061].
  * [44] D. Budker, P.W. Graham, H. Ramani, F. Schmidt-Kaler, C. Smorra and S. Ulmer, _Millicharged Dark Matter Detection with Ion Traps_ , _PRX Quantum_ 3 (2022) 010330 [2108.05283].
  * [45] J. Billard, M. Pyle, S. Rajendran and H. Ramani, _Calorimetric Detection of Dark Matter_ , 2208.05485.
  * [46] D. McKeen, M. Moore, D.E. Morrissey, M. Pospelov and H. Ramani, _Accelerating Earth-Bound Dark Matter_ , 2202.08840.
  * [47] D. McKeen, D.E. Morrissey, M. Pospelov, H. Ramani and A. Ray, _Dark Matter Annihilation inside Large-Volume Neutrino Detectors_ , _Phys. Rev. Lett._ 131 (2023) 011005 [2303.03416].
  * [48] A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov and A.K. Geim, _The electronic properties of graphene_ , _Reviews of Modern Physics_ 81 (2009) 109 [0709.1163].
  * [49] S.-Y. Wang, _Graphene-based detectors for directional dark matter detection_ , _Eur. Phys. J. C_ 79 (2019) 561 [1509.08801].
  * [50] Y. Hochberg, Y. Kahn, M. Lisanti, C.G. Tully and K.M. Zurek, _Directional detection of dark matter with two-dimensional targets_ , _Phys. Lett. B_ 772 (2017) 239 [1606.08849].
  * [51] D. Kim, J.-C. Park, K.C. Fong and G.-H. Lee, _Detection of Super-light Dark Matter Using Graphene Sensor_ , 2002.07821.
  * [52] R. Catena, T. Emken, M. Matas, N.A. Spaldin and E. Urdshals, _Direct searches for general dark matter-electron interactions with graphene detectors: Part II. Sensitivity studies_ , 2303.15509.
  * [53] R. Catena, T. Emken, M. Matas, N.A. Spaldin and E. Urdshals, _Direct searches for general dark matter-electron interactions with graphene detectors: Part I. Electronic structure calculations_ , 2303.15497.
  * [54] E.V. Castro, N.M.R. Peres, J.M.B. Lopes dos Santos, F. Guinea and A.H. Castro Neto, _Bilayer graphene: gap tunability and edge properties_ , in _Journal of Physics Conference Series_ , vol. 129 of _Journal of Physics Conference Series_ , p. 012002, Oct., 2008, DOI [1004.5079].
  * [55] E. McCann and M. Koshino, _The electronic properties of bilayer graphene_ , _Reports on Progress in Physics_ 76 (2013) 056503 [1205.6953].
  * [56] A.V. Rozhkov, A.O. Sboychakov, A.L. Rakhmanov and F. Nori, _Electronic properties of graphene-based bilayer systems_ , _Phys. Rep._ 648 (2016) 1 [1511.06706].
  * [57] E. McCann, _Asymmetry gap in the electronic band structure of bilayer graphene_ , _Phys. Rev. B_ 74 (2006) 161403.
  * [58] H. Min, B. Sahu, S.K. Banerjee and A.H. MacDonald, _Ab initio theory of gate induced gaps in graphene bilayers_ , _Phys. Rev. B_ 75 (2007) 155115.
  * [59] J.B. Oostinga, H.B. Heersche, X. Liu, A.F. Morpurgo and L.M.K. Vandersypen, _Gate-induced insulating state in bilayer graphene devices_ , _Nature Materials_ 7 (2008) 151.
  * [60] E.V. Castro, K.S. Novoselov, S.V. Morozov, N.M.R. Peres, J.M.B.L. dos Santos, J. Nilsson et al., _Biased bilayer graphene: Semiconductor with a gap tunable by the electric field effect_ , _Phys. Rev. Lett._ 99 (2007) 216802.
  * [61] Z. Liu, K. Suenaga, P.J.F. Harris and S. Iijima, _Open and closed edges of graphene layers_ , _Phys. Rev. Lett._ 102 (2009) 015501.
  * [62] P.L. de Andres, R. Ramírez and J.A. Vergés, _Strong covalent bonding between two graphene layers_ , _Phys. Rev. B_ 77 (2008) 045403.
  * [63] R. Bistritzer and A.H. MacDonald, _Moiré bands in twisted double-layer graphene_ , _Proceedings of the National Academy of Sciences_ 108 (2011) 12233.
  * [64] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras et al., _Unconventional superconductivity in magic-angle graphene superlattices_ , _Nature_ 556 (2018) 43 [1803.02342].
  * [65] Y. Cao, V. Fatemi, A. Demir, S. Fang, S.L. Tomarken, J.Y. Luo et al., _Correlated insulator behaviour at half-filling in magic-angle graphene superlattices_ , _Nature_ 556 (2018) 80.
  * [66] E.Y. Andrei and A.H. MacDonald, _Graphene bilayers with a twist_ , _Nature Materials_ 19 (2020) 1265.
  * [67] T. Ohta, A. Bostwick, T. Seyller, K. Horn and E. Rotenberg, _Controlling the electronic structure of bilayer graphene_ , _Science_ 313 (2006) 951 [https://www.science.org/doi/pdf/10.1126/science.1130681].
  * [68] E. McCann and M. Koshino, _The electronic properties of bilayer graphene_ , _Reports on Progress in Physics_ 76 (2013) 056503.
  * [69] J. Jung and A.H. MacDonald, _Accurate tight-binding models for the $\pi$ bands of bilayer graphene_, _Phys. Rev. B_ 89 (2014) 035405.
  * [70] T. Trickle, Z. Zhang, K.M. Zurek, K. Inzani and S.M. Griffin, _Multi-Channel Direct Detection of Light Dark Matter: Theoretical Framework_ , _JHEP_ 03 (2020) 036 [1910.08092].
  * [71] Y. Kahn and T. Lin, _Searches for light dark matter using condensed matter systems_ , _Rept. Prog. Phys._ 85 (2022) 066901 [2108.03239].
  * [72] J.I. Read, _The Local Dark Matter Density_ , _J. Phys. G_ 41 (2014) 063101 [1404.1938].
  * [73] S. Knapen, J. Kozaczuk and T. Lin, _Migdal Effect in Semiconductors_ , _Phys. Rev. Lett._ 127 (2021) 081805 [2011.09496].
  * [74] Y. Hochberg, Y. Kahn, N. Kurinsky, B.V. Lehmann, T.C. Yu and K.K. Berggren, _Determining Dark-Matter–Electron Scattering Rates from the Dielectric Function_ , _Phys. Rev. Lett._ 127 (2021) 151802 [2101.08263].
  * [75] G.F. Giuliani and G. Vignale, _Quantum Theory of the Electron Liquid_ , Cambridge University Press, Cambridge U.K. (2005).
  * [76] M. Girvin and K. Yang, _Modern Condensed Matter Physics_ , Cambridge University Press, Cambridge U.K. (2019).
  * [77] G.D. Mahan, _Many-Particle Physics_ , Springer, New York USA (2000).
  * [78] H. Bruus and K. Flensber, _Many-Body Quantum Theory in Condensed Matter Physics: An Introduction_ , Oxford University Press, Oxford U.K. (2004).
  * [79] N. Prasad, G.W. Burg, K. Watanabe, T. Taniguchi, L.F. Register and E. Tutuc, _Quantum lifetime spectroscopy and magnetotunneling in double bilayer graphene heterostructures_ , _Phys. Rev. Lett._ 127 (2021) 117701.
  * [80] G. Li, A. Luican and E.Y. Andrei, _Scanning tunneling spectroscopy of graphene on graphite_ , _Phys. Rev. Lett._ 102 (2009) 176804.
  * [81] G. Moos, C. Gahl, R. Fasel, M. Wolf and T. Hertel, _Anisotropy of quasiparticle lifetimes and the role of disorder in graphite from ultrafast time-resolved photoemission spectroscopy_ , _Phys. Rev. Lett._ 87 (2001) 267402.
  * [82] A. Mitridate, T. Trickle, Z. Zhang and K.M. Zurek, _Dark matter absorption via electronic excitations_ , _JHEP_ 09 (2021) 123 [2106.12586].
  * [83] J.D. Jackson, _Classical Electrodynamics_ , Wiley, New York (1999).
  * [84] S. Knapen, J. Kozaczuk and T. Lin, _Dark matter-electron scattering in dielectrics_ , _Phys. Rev. D_ 104 (2021) 015031 [2101.08275].
  * [85] D. Brida, A. Tomadin, C. Manzoni, Y.J. Kim, A. Lombardo, S. Milana et al., _Ultrafast collinear scattering and carrier multiplication in graphene_ , _Nature Communications_ 4 (2013) 1987 [1209.5729].
  * [86] K.J. Tielrooij, J.C.W. Song, S.A. Jensen, A. Centeno, A. Pesquera, A. Zurutuza Elorza et al., _Photoexcitation cascade and multiple hot-carrier generation in graphene_ , _Nature Physics_ 9 (2013) 248 [1210.1205].
  * [87] F.H.L. Koppens, T. Mueller, P. Avouris, A.C. Ferrari, M.S. Vitiello and M. Polini, _Photodetectors based on graphene, other two-dimensional materials and hybrid systems_ , _Nature Nanotechnology_ 9 (2014) 780.
  * [88] J.O.D. Williams, J.A. Alexander-Webber, J.S. Lapington, M. Roy, I.B. Hutchinson, A.A. Sagade et al., _Towards a graphene-based low intensity photon counting photodetector_ , _Sensors_ 16 (2016) .
  * [89] A. De Sanctis, J.D. Mehew, M.F. Craciun and S. Russo, _Graphene-based light sensing: fabrication, characterisation, physical properties and performance_ , _Materials_ 11 (2018) 1762 [1809.04862].
  * [90] P. Seifert, X. Lu, P. Stepanov, J.R. DurÃ¡n Retamal, J.N. Moore, K.-C. Fong et al., _Magic-angle bilayer graphene nanocalorimeters: Toward broadband, energy-resolving single photon detection_ , _Nano Letters_ 20 (2020) 3459 [https://doi.org/10.1021/acs.nanolett.0c00373].
  * [91] S. Betts et al., _Development of a Relic Neutrino Detection Experiment at PTOLEMY: Princeton Tritium Observatory for Light, Early-Universe, Massive-Neutrino Yield_ , in _Snowmass 2013: Snowmass on the Mississippi_ , 7, 2013 [1307.4738].
  * [92] “Ptolemy: Relic neutrino direct detection.” https://indico.cern.ch/event/1199289/contributions/5447612/.
  * [93] SENSEI collaboration, _SENSEI: Characterization of Single-Electron Events Using a Skipper Charge-Coupled Device_ , _Phys. Rev. Applied_ 17 (2022) 014022 [2106.08347].
  * [94] G. Papadopoulos, _Using scientific-grade ccds for the direct detection of dark matter with the damic-m experiment_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 1040 (2022) 167184.
  * [95] Y. Hochberg, B.V. Lehmann, I. Charaev, J. Chiles, M. Colangelo, S.W. Nam et al., _New constraints on dark matter from superconducting nanowires_ , _Phys. Rev. D_ 106 (2022) 112005 [2110.01586].
  * [96] DAMIC collaboration, _Constraints on Light Dark Matter Particles Interacting with Electrons from DAMIC at SNOLAB_ , _Phys. Rev. Lett._ 123 (2019) 181802 [1907.12628].
  * [97] FUNK Experiment collaboration, _Limits from the Funk Experiment on the Mixing Strength of Hidden-Photon Dark Matter in the Visible and Near-Ultraviolet Wavelength Range_ , _Phys. Rev. D_ 102 (2020) 042001 [2003.13144].
  * [98] XENON collaboration, _Light Dark Matter Search with Ionization Signals in XENON1T_ , _Phys. Rev. Lett._ 123 (2019) 251801 [1907.11485].
  * [99] H. An, M. Pospelov, J. Pradler and A. Ritz, _Direct Detection Constraints on Dark Photon Dark Matter_ , _Phys. Lett. B_ 747 (2015) 331 [1412.8378].
  * [100] H. An, M. Pospelov, J. Pradler and A. Ritz, _New limits on dark photons from solar emission and keV scale dark matter_ , _Phys. Rev. D_ 102 (2020) 115022 [2006.13929].
  * [101] J. Buch, M.A. Buen-Abad, J. Fan and J.S.C. Leung, _Dark Matter Substructure under the Electron Scattering Lamppost_ , _Phys. Rev. D_ 102 (2020) 083010 [2007.13750].
  * [102] T.N. Maity and R. Laha, _Dark matter substructures affect dark matter-electron scattering in xenon-based direct detection experiments_ , _JHEP_ 02 (2023) 200 [2208.14471].
  * [103] PTOLEMY collaboration, _Neutrino physics with the PTOLEMY project: active neutrino properties and the light sterile case_ , _JCAP_ 07 (2019) 047 [1902.05508].
  * [104] H. Min and A.H. MacDonald, _Electronic Structure of Multilayer Graphene_ , _Progress of Theoretical Physics Supplement_ 176 (2008) 227.
  * [105] J. Nilsson, A.H. Castro Neto, F. Guinea and N.M.R. Peres, _Electronic properties of bilayer and multilayer graphene_ , _Phys. Rev. B_ 78 (2008) 045405.
  * [106] N. Regnault, Y. Xu, M.-R. Li, D.-S. Ma, M. Jovanovic, A. Yazdani et al., _Catalogue of flat-band stoichiometric materials_ , _Nature_ 603 (2022) 824.
  * [107] M. Kang, S. Fang, L. Ye, H.C. Po, J. Denlinger, C. Jozwiak et al., _Topological flat bands in frustrated kagome lattice CoSn_ , _Nature Communications_ 11 (2020) 4004.
  * [108] M. Kang, L. Ye, S. Fang, J.-S. You, A. Levitan, M. Han et al., _Dirac fermions and flat bands in the ideal kagome metal FeSn_ , _Nature Materials_ 19 (2020) 163.