of 15
Dark matter direct detection in materials with spin-orbit coupling
Hsiao-Yi Chen,
1,2
Andrea Mitridate ,
3
Tanner Trickle ,
3
Zhengkang Zhang,
4
Marco Bernardi ,
1
and Kathryn M. Zurek
3
1
Department of Applied Physics and Materials Science, California Institute of Technology,
Pasadena, California 91125, USA
2
RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama, 351-0198, Japan
3
Walter Burke Institute for Theoretical Physics, California Institute of Technology,
Pasadena, California 91125, USA
4
Department of Physics, University of California, Santa Barbara, California 93106, USA
(Received 15 March 2022; accepted 13 July 2022; published 25 July 2022)
Semiconductors with
O
ð
meV
Þ
band gaps have been shown to be promising targets to search for sub-
MeV mass dark matter (DM). In this paper we focus on a class of materials where such narrow band gaps
arise naturally as a consequence of spin-orbit coupling (SOC). Specifically, we are interested in computing
DM-electron scattering and absorption rates in these materials using state-of-the-art density functional
theory techniques. To do this, we extend the DM interaction rate calculation to include SOC effects which
necessitates a generalization to spin-dependent wave functions. We apply our new formalism to calculate
limits for several DM benchmark models using an example ZrTe
5
target and show that the inclusion of
SOC can substantially alter projected constraints.
DOI:
10.1103/PhysRevD.106.015024
I. INTRODUCTION
The detection of dark matter (DM) through nongravita-
tional interactions remains one of the main goals of particle
physics. Electronic excitations have been identified as a
promising path to lead the direct detection of DM to
sub-GeV masses, a region not kinematically accessible
in experiments based on nuclear recoil. A variety of
avenues to search for DM-induced electronic excitations
have been proposed: ionization in noble gases
[1
8]
,
excitations across a band gap in crystal targets
[1
3,9
19]
, superconductors
[20,21]
, graphene
[22]
, Dirac materi-
als
[23
27]
, and transitions between molecular orbitals in
aromatic organics
[28,29]
.
In this work we focus on a specific class of semi-
conductors for which
O
ð
meV
Þ
band gaps [as opposed to
typical
O
ð
eV
Þ
band gaps in semiconductors and insulators]
arise as a consequence of spin-orbit coupling (SOC) effects.
Targets with such small band gaps can probe DM masses
down to
O
ð
keV
Þ
via scattering, and
O
ð
meV
Þ
via absorp-
tion, while still suppressing thermal noise. Moreover, some
of these SOC materials have tunable band structures, a
property which makes them interesting candidates for
direct-detection experiments
[27]
.
However, SOC effects introduce some intricacies in the
DM-electron interaction rate calculations since the Bloch
wave functions are no longer eigenstates of the
S
z
operator,
and therefore become two-component objects in spin space.
This implies that electron spin sums cannot be trivially
performed, and new transition form factors must be
computed. For example, spin-dependent vector-mediated
scattering can no longer be related to its spin-independent
counterpart, and must be computed from first principles.
We extend the framework in Sec.
II
and implement the new
spin-dependent form factors numerically within EXCEED-
DM
[30,31]
, which is publicly available on Github.
To showcase the formalism developed in this paper, we
apply it to a target with important SOC effects: ZrTe
5
.This
material has been extensively studied in the context of DM
direct detection
[24
26]
as a leading candidate for Dirac
material targets. Dirac materials are characterized by low-
energy excitations which behave like free electrons and
satisfy the Dirac equation. The properties of the electronic
excitations can then be understood by a simple extension of
the standard QED results. An additional consequence is that
they have weak electromagnetic screening, even with a small
energy gap between the valence and conduction bands,
making them a desirable target for sub-MeV dark matter
coupled to electrons via a dark photon mediator
[25]
. In this
work,however,wefocusonlyonZrTe
5
properties that stem
from its SOC nature, and we will not exploit any of those
derivingfromitsDirac nature (which is stilldebated
[32
43]
).
To illustrate the variety of DM models that an SOC target
can probe, we consider several different DM models and
processes. Specifically, we study the following.
Published by the American Physical Society under the terms of
the
Creative Commons Attribution 4.0 International
license.
Further distribution of this work must maintain attribution to
the author(s) and the published article
s title, journal citation,
and DOI. Funded by SCOAP
3
.
PHYSICAL REVIEW D
106,
015024 (2022)
2470-0010
=
2022
=
106(1)
=
015024(15)
015024-1
Published by the American Physical Society
(1) Standard spin-independent (SI) and spin-dependent
(SD) scattering via vector mediators: The funda-
mental interaction Lagrangians for these models take
the form
L
int
¼
(
φ
μ
ð
g
χ
̄
χγ
μ
χ
þ
g
e
̄
ψγ
μ
ψ
Þð
SI
Þ
;
φ
μ
ð
g
χ
̄
χγ
μ
γ
5
χ
þ
g
e
̄
ψγ
μ
γ
5
ψ
Þð
SD
Þ
;
ð
1
Þ
where
ψ
and
χ
are the electron and DM fermion
fields, respectively, and
φ
μ
is the dark mediator field.
(2) Scalar, pseudoscalar, and vector DM absorption: In
this case, the fundamental interaction Lagrangians
take the form
L
int
¼
8
>
>
<
>
>
:
g
e
φ
̄
ψψ
ð
scalar DM
Þ
;
g
e
φ
̄
ψ
i
γ
5
ψ
ð
pseudoscalar DM
Þ
;
g
e
φ
μ
̄
ψγ
μ
ψ
ð
vector DM
Þ
:
ð
2
Þ
The paper is organized as follows. In Sec.
II
we extend
the DM interaction rate formalism to account for spin-
dependent wave functions in general (spin-orbit coupled,
anisotropic) targets. Then, in Sec.
III
we apply these results
to the candidate material ZrTe
5
and compare the results
obtained with and without the inclusion of SOC effects.
Further details of the density functional theory (DFT)
calculation are presented in Appendix
A1
, and conver-
gence tests for the results shown in Sec.
III
can be found in
Appendix
A2
.
II. DM INTERACTION RATE FORMALISM
In this section we derive the rates for transitions between
electronic energy levels induced by DM absorption and
scattering. For the targets of interest here, the electronic
energy levels can be labeled by a band index
i
and a
momentum
k
within the first Brillouin zone (1BZ), which
we collectively indicate with an index
I
¼f
i;
k
g
. The wave
functions of the electronic states can be written in the Bloch
form as
Ψ
I
ð
x
Þ¼
1
ffiffiffiffi
V
p
e
i
k
·
x
u
I
ð
x
Þ
;
ð
3
Þ
where the periodic Bloch wave functions
u
I
are two-
component vectors in the spin basis, and
V
is the crystal
volume.
A. Absorption
In this subsection, we use the nonrelativistic (NR)
effective filed theory (EFT) developed in Ref.
[19]
, and
summarized in Appendix
B
, to compute DM absorption
rates in materials with sizable SOC.
The absorption rate of a state can be derived from the
imaginary part of its self-energy. In a medium, care must be
taken due to the possible mixing between the DM,
φ
, and
standard model (SM) states (in our case the SM photon,
A
).
In the presence of such mixing effects, the DM absorption
rate is related to the imaginary part of the self-energy of the
mostly DM
eigenstate,
Π
ˆ
φ
ˆ
φ
:
Γ
φ
abs
¼
Z
ˆ
φ
ω
Im
Π
ˆ
φ
ˆ
φ
;
ð
4
Þ
where
ω
m
φ
is the energy of the DM state, and
Z
ˆ
φ
¼
ð
1
d
Re
Π
ˆ
φ
ˆ
φ
d
ω
2
Þ
1
¼
1
þ
O
ð
g
2
e
Þ
is the wave-function renorm-
alization which we will approximate as unity in the
following. The total absorption rate per unit target mass,
R
,isgivenby
R
¼
ρ
φ
ρ
T
m
φ
1
n
X
n
η
¼
1
Γ
φ
η
abs
;
ð
5
Þ
where
n
is the number of degrees of freedom of the DM
particle (
n
¼
3
for vector DM and
n
¼
1
for scalar and
pseudoscalar DM) and we average over the incoming DM
polarizations. The DM density
ρ
φ
is taken to be
0
.
4
GeVcm
3
, and
ρ
T
is the target density.
To derive
Π
ˆ
φ
ˆ
φ
we need to diagonalize the in-medium
self-energy matrix, which in our case contains a mixing
between the DM and the SM photon:
S
eff
⊃−
1
2
Z
d
4
Q
ð
A
λ
φ
η
Þ
×
Π
λλ
0
AA
Π
λη
0
A
φ
Π
ηλ
0
φ
A
m
2
φ
δ
ηη
0
þ
Π
ηη
0
φφ
!

A
λ
0
φ
η
0

;
ð
6
Þ
where the implicit sum over
λ
;
λ
0
(
η
;
η
0
) runs over the photon
(DM) polarizations, and we have dropped the Lorentz
indices on
A
λ
(and
φ
λ
in the case of vector DM) for
simplicity. We have also introduced components of the self-
energy projected on the polarization vectors, defined as
Π
λλ
0
ε
λ
μ
Π
μν
ε
λ
0

ν
;
ð
7
Þ
where
ε
λ
μ
are polarization vectors. In general, the polari-
zation vectors that diagonalize this matrix are not the
typical longitudinal and transverse polarization vectors,
since mixing can occur (i.e.,
ε
λ
T;
μ
Π
μν
AA
ε
λ
0

L;
ν
0
, where
ε
T
,
ε
L
are the transverse and longitudinal polarization vectors,
respectively). However, one can always find an appropriate
basis to diagonalize the DM and photon self-energies. In
this basis, Eq.
(6)
becomes
HSIAO-YI CHEN
et al.
PHYS. REV. D
106,
015024 (2022)
015024-2
S
eff
⊃−
1
2
Z
d
4
Q
ð
A
λ
φ
η
Þ
×
Π
λ
AA
δ
λλ
0
Π
λη
0
A
φ
Π
ηλ
0
φ
A
ð
m
2
φ
þ
Π
η
φφ
Þ
δ
ηη
0
!

A
λ
0
φ
η
0

;
ð
8
Þ
where
Π
λ
AA
and
Π
η
φφ
are the eigenvalues of
Π
λλ
0
AA
and
Π
ηη
0
φφ
,
respectively.
The off-diagonal terms in Eq.
(8)
are perturbatively
suppressed by a factor of
g
e
with respect to the
Π
AA
terms.
Therefore, working at order
O
ð
g
2
e
Þ
, we find that the in-
medium self-energy for the
η
polarization of the mostly DM
eigenstate is given by
Π
η
ˆ
φ
ˆ
φ
¼
Π
η
φφ
þ
X
λ
Π
ηλ
φ
A
Π
λη
A
φ
m
2
φ
Π
λ
AA
:
ð
9
Þ
Since vector DM couples to electrons in the same way as
the photon, one can derive the relevant self-energies by
simply replacing the electromagnetic charge with
g
e
, e.g.,
Π
ηλ
φ
A
¼
ð
g
e
=e
Þ
Π
λ
AA
δ
ηλ
. Doing so allows us to write Eq.
(9)
in terms of the photon self-energy as
Π
η
ˆ
φ
ˆ
φ
¼

g
e
e

2
m
2
φ
Π
η
AA
m
2
φ
Π
η
AA
ð
vector DM
Þ
:
ð
10
Þ
Scalar and pseudoscalar DM only have one degree of
freedom, and therefore Eq.
(9)
takes the form
Π
ˆ
φ
ˆ
φ
¼
Π
φφ
þ
X
λ
Π
λ
A
φ
Π
λ
φ
A
m
2
φ
Π
λ
AA
½ð
pseudo
Þ
scalar DM

:
ð
11
Þ
As usual, the self-energies appearing in the previous
equations are computed from the sum of one-particle
irreducible diagrams. Working at one loop, there are two
graph topologies that can contribute:
ð
12
Þ
ð
13
Þ
where
O
ð
1
;
2
Þ
is any operator coupling the external field,
A
or
φ
(dashed lines), to the electron (solid lines). For vector
external states these operators carry Lorentz indices that are
inherited by
̄
Π
and
̄
Π
0
.
The full expressions for the self-energies involved in the
absorption calculation can be found in Appendix
B
.
However, as we discuss in the same appendix, due to
the absorption kinematics and the hierarchy between the
DM and electron velocities, a few diagrams dominate these
self-energies. Specifically, we find that the diagonalization
of the photon in-medium self-energy (and therefore the
derivation of
Π
λ
AA
) reduces to diagonalizing
̄
Π
v
i
;v
j
, where
the velocity operator is defined by
v
i
i
i
2
m
e
:
ð
14
Þ
From this it follows that the long-wavelength limit of the
dielectric function
ε
ð
0
;
ω
Þ
, which will enter explicitly in the
scattering rate calculation, can be derived from
̄
Π
v
i
;v
j
:
½
ε
ð
0
;
ω
Þ
ij
¼
1
þ
Π
ij
AA
ω
2
1
e
2
̄
Π
v
i
;v
j
ω
2
:
ð
15
Þ
The long-wavelength dielectric function
ε
ð
0
;
ω
Þ
, together
with details of its numeric calculation, are reported in
Appendix
A2
.
1
For scalar and pseudoscalar DM, the
leading-order terms in the self-energy of the mostly DM
eigenstate are found to be
Π
ˆ
φ
ˆ
φ

g
2
e
̄
Π
̄
v
2
;
̄
v
2
ð
scalar DM
Þ
;
g
2
e
ω
2
4
m
2
e
̄
Π
v
·
σ
;
v
·
σ
ð
pseudoscalar DM
Þ
;
ð
16
Þ
where we have introduced the operator
̄
v
2
2
8
m
2
e
:
ð
17
Þ
B. Scattering
In this subsection we proceed to derive the DM scatter-
ing rate with spin-dependent electronic wave functions.
Generalizing the formulas previously derived in
Refs.
[2,5,10,15,31,44]
, we can write the DM scattering
rate as
1
Strictly speaking, the dielectric function is a mixed-index
tensor, as evident from the defining equation,
J
i
¼
σ
i
j
A
j
¼
i
ω
ð
1
ε
i
j
Þ
A
j
, where
J
i
is the electronic current density,
σ
i
j
is
the conductivity tensor, and
ε
i
j
¼
δ
i
j
Π
i
j
=
ω
2
¼
1
þ
Π
ij
=
ω
2
(see
the discussion in Appendix A f Ref.
[25]
for more details). With a
slight abuse of notation, we define the matrix
ε
which has
components
½
ε

ij
¼
ε
i
j
.
DARK MATTER DIRECT DETECTION IN MATERIALS WITH
...
PHYS. REV. D
106,
015024 (2022)
015024-3
Γ
I
I
0
¼
π
8
Vm
2
e
m
2
χ
Z
d
3
q
δ
ð
E
I
0
E
I
ω
q
Þ
×




Z
d
3
k
ð
2
π
Þ
3
̃
Ψ

I
0
ð
k
þ
q
Þ
·
M
ð
q
Þ
·
̃
Ψ
I
ð
k
Þ




2
;
ð
18
Þ
where the bar indicates a spin average (sum) over the
incoming (outgoing) DM states,
̃
Ψ
I
ð
k
Þ
are the Fourier
transforms of the electronic wave functions defined in
Eq.
(3)
, and
ω
q
q
·
v
q
2
2
m
χ
:
ð
19
Þ
For the SI and SD models of interest here, we can write
the free electron scattering amplitude as
M
ss
0
;
σσ
0
ð
q
Þ¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
16
π
m
2
χ
m
2
e
̄
σ
e
μ
2
χ
e
s
f
e
f
0
e
F
med

q
q
0

S
ss
0
;
σσ
0
;
ð
20
Þ
where
F
med
ð
q
q
0
Þ
encodes the momentum dependence
induced by the mediator propagator,
f
e
ð
q
Þ
=f
0
e
is the
screening factor introduced by in-medium effects (defined
explicitly in Sec.
III B
), and
̄
σ
e
is a reference cross section
defined by
̄
σ
e
μ
2
χ
e
64
π
m
2
χ
m
2
e
X
ss
0
;
σσ
0
j
M
ss
0
;
σσ
0
ð
q
0
Þj
2
;
ð
21
Þ
with
q
0
¼
α
m
e
. The total rate per unit detector mass is then
R
¼
π
̄
σ
e
V
μ
2
χ
e
m
χ
ρ
χ
ρ
T
X
I;I
0
Z
d
3
q
ð
2
π
Þ
3

f
e
f
0
e

2
F
2
med
ð
q
Þ
g
ð
q
;
ω
Þ
F
II
0
ð
q
Þ
;
ð
22
Þ
where
g
ð
q
;
ω
Þ
is the velocity integral defined as
g
ð
q
;
ω
Þ
Z
d
3
v
f
χ
ð
v
Þ
2
πδ
ð
ω
ω
q
Þ
;
ð
23
Þ
with
f
χ
ð
v
Þ
being the DM velocity distribution in the
laboratory rest frame, which we take to be a boosted
Maxwell-Boltzmann distribution with parameters
v
0
¼
230
km s
1
,
v
esc
¼
600
km s
1
, and
v
e
¼
240
km s
1
.
The crystal form factor
F
II
0
is defined as
F
II
0
ð
q
Þ




Z
d
3
k
ð
2
π
Þ
3
̃
Ψ

I
0
ð
k
þ
q
Þ
·
S
·
̃
Ψ
I
ð
k
Þ




2
;
ð
24
Þ
where the spin operators for the models considered in this
work are given by
S
ss
0
;
σσ
0
¼
(
δ
ss
0
δ
σσ
0
ð
SI
Þ
;
1
ffiffi
3
p
P
i
σ
i
ss
0
σ
i
σσ
0
ð
SD
Þ
;
ð
25
Þ
and
σ
i
are the Pauli matrices. Given these expressions, the
form factors
F
II
0
for the SI and SD models take the form
2
F
II
0
¼

j
T
II
0
j
2
ð
SI
Þ
;
1
3
T

II
0
·
T
II
0
ð
SD
Þ
;
ð
26
Þ
where we have defined the DM model-independent tran-
sition form factors
T
II
0
and
T
II
0
as
T
II
0
¼
Z
d
3
k
ð
2
π
Þ
3
̃
Ψ

I
0
ð
k
þ
q
Þ
·
̃
Ψ
I
ð
k
Þ
;
ð
27
Þ
T
II
0
¼
Z
d
3
k
ð
2
π
Þ
3
̃
Ψ

I
0
ð
k
þ
q
Þ
·
σ
·
̃
Ψ
I
ð
k
Þ
:
ð
28
Þ
III. DETECTION RATES IN
ZrTe
5
We will now apply the formalism developed in the
previous section to our benchmark SOC target: ZrTe
5
. The
band structure of ZrTe
5
, with and without the inclusion of
SOC effects, is shown in Fig.
1
. The details of the DFT
calculation can be found in Appendix
A1
. The dominant
effect of SOC is to shift the valence and conduction bands
closer at the
Γ
point relative to the No SOC calculation.
While in theory the calculation of DM interaction rates is
identical for
O
ð
meV
Þ
and
O
ð
eV
Þ
gap semiconductors, in
practice one must be careful about sampling the 1BZ. This
is because these
O
ð
meV
Þ
energy differences generally only
occur in small volumes within the 1BZ. To account for this
we sample the 1BZ with a higher
k
-point density in regions
corresponding to the low-energy band structure. For ZrTe
5
this occurs near the
Γ
point, and we split the phase space
into two separate regions,
low
E
and
high
E
,
which we
describe now. The low-
E
region consists of the highest two
(one) valence bands and lowest two (one) conduction
bands, for the calculation with (without) SOC, sampled
on a
mini-BZ
grid. This mini-BZ grid is a rescaled
uniform Monkhorst-Pack grid
[46]
; each
k
is scaled by a
factor of
1
=
5
, giving a
125
×
k
-point sampling in that
region. This region will give the dominant contribution to
absorption of DM with mass
100
meV, as well as
DM scattering via a light mediator. The high-
E
region
includes all the bands outside the low-
E
region, sampled
with a standard Monkhorst-Pack uniform grid. The DM
2
The absence of the overall factor of 2, relative to the SI rate
formula given in Ref.
[31]
, canbe understood from the sum over the
states. If the wave functions are spin independent, then
P
IF
P
IF
P
ss
0
, where
s
(
s
0
) indexes the initial (final) electron
spin state. These spin sums contribute the extra factor of 2, bringing
Eq.
(22)
and the rate formula in Ref.
[31]
into agreement.
HSIAO-YI CHEN
et al.
PHYS. REV. D
106,
015024 (2022)
015024-4
absorption rates in Sec.
III A
are a combination of the low-
and high-
E
regions. The DM scattering rates in Sec.
III B
will be shown for both regions, and it will be clear when
one dominates the other.
We compute the Bloch wave functions
(3)
in both regions
within the framework of DFT with Quantum ESPRESSO
[47
49]
; details can be found in Appendix
A1
.TheDM
absorption and scattering rates are computed with an
extended version of EXCEED-DM
[30,31]
which includes
the formalism developed in Sec.
II
, and is publicly available
on Github.
For each of the models considered in the following
subsections, we show the projected constraints from three
different calculations. The curves labeled
SOC
are
computed with the inclusion of SOC effects, the curves
labeled
No SOC
do not include any SOC effects, and
those labeled
Partial SOC
are a combination of the SOC
and No SOC calculations, obtained using the energy levels
computed with SOC, and the wave functions without SOC.
While the Partial SOC results are not a consistent calcu-
lation, they aid in understanding how much of the differ-
ence between the SOC and No SOC results is due to the
changes in the band structure versus the inclusion of the
spin-dependent wave functions. Generally, we find that
the changes in the band structure are more influential than
the spin dependence in the wave functions, but the latter can
still be important.
Last, we note that previous works
[24
26]
derived exci-
tation rates analytically by exploiting the putative Dirac
nature of ZrTe
5
. While a direct comparison to assess the
validity ofthe analytic approximations is dubious sincewe do
not observe a conical band structure (see Appendix
A1
),
previous estimates from Ref.
[25]
are shown in Figs.
2
and
3
.
In Appendix
C
we discuss the validity of these analytic
approximations in more general Dirac materials.
A. Absorption
For the models considered [Eq.
(2)
] our results are shown
in Fig.
2
. For ease of comparison we map the constraints on
the
g
e
parameters in Eq.
(2)
to a more commonly used
notation,
FIG. 2. Comparison of projected 95% C.L. reach (three events, no background) assuming
1
kg-year exposure for scalar (left),
pseudoscalar (center), and vector (right) DM. We compare our results with (solid) and without (dotted) SOC for electronic absorption in
a ZrTe
5
target (red), with those for semiconductor silicon (Si, blue) and germanium (Ge, green) targets
[19]
, superconducting aluminum
(Al-SC, brown)
[19]
), phononic absorption in polar materials
[50,51]
(GaAs in orange and SiO
2
in purple), and previous estimates for
ZrTe
5
(teal)
[25]
. We also show the projected constraints combining the SOC energy levels with the No SOC wave functions (
Partial
SOC,
red dashed) to explicitly show the effect of the spin-dependent wave functions. Constraints are expressed in terms of the
commonly adopted parameters shown in Eq.
(29)
. Shaded red bands correspond to different parametrizations of the electron width
δ
½
10
1
.
5
;
10
0
.
5

ω
used in calculating the self-energies [see, e.g., Eq.
(B17)
], with the solid line corresponding to
δ
¼
10
1
ω
. Thin
lines indicate results obtained by rescaling the optical data. Also shown are the direct-detection limits from XENON
10
=
100
[13]
, fifth
force constraints
[52]
, and stellar cooling constraints from red giants (RG)
[53]
, and white dwarfs (WD)
[54]
. For the pseudoscalar
scenario we also report the couplings corresponding to the QCD axion in KSVZ and DFSZ models, for
0
.
28
tan
β
140
[55]
.
FIG. 1. ZrTe
5
band structure computed using DFT with SOC
(solid lines) and without SOC (dashed lines). The inset highlights
the low-energy (low
E
) band dispersion, whose details are
sampled using a denser
k
-point grid. The band gap for the
SOC band structure is set to the experimental value of 23.5 meV
[45]
, and the No SOC band structure is shifted accordingly, which
gives a larger band gap of 81.6 meV.
DARK MATTER DIRECT DETECTION IN MATERIALS WITH
...
PHYS. REV. D
106,
015024 (2022)
015024-5
g
e
¼
8
>
>
<
>
>
:
ffiffiffiffi
4
π
p
m
e
M
Pl
d
φ
ee
ð
scalar
Þ
;
g
aee
ð
pseudoscalar
Þ
;
e
κ
ð
vector
Þ
;
ð
29
Þ
where
M
Pl
¼
1
.
22
×
10
19
GeV is the Planck mass.
For all of the benchmark models, the inclusion of SOC
effects dominantly impacts the low-mass reach where the
SOC corrections to the band structure are most relevant.
Most notably, the lowest testable DM mass is shifted as a
consequence of the different band gaps: 23.5 meV with
SOC, and 81.6 meV without SOC. At higher masses the
SOC effects are milder and, as expected, the SOC reach
approaches the reach without SOC effects. The close
agreement between the
Partial SOC
and
SOC
curves
indicates that changes to the energy levels are what is
mainly driving the difference in the
SOC
and
No SOC
calculations.
For the scalar and vector DM models we find that ZrTe
5
is superior at low DM masses relative to a superconducting
aluminum target, another target material with an
O
ð
meV
Þ
gap (0.6 meV for the Al-SC curves shown here). However,
for pseudoscalar DM, for
m
φ
eV, Al-SC yields better
sensitivity than ZrTe
5
. This can be attributed to the large
amount of screening present in the vector DM case (but not
in the pseudoscalar DM case) for Al-SC.
Shaded bands correspond to different parametrizations
of the electron width,
δ
½
10
1
.
5
;
10
0
.
5

ω
, used when
computing the self-energies [see Eq.
(B17)
for more
details]. In theory, the absorption rate calculation is
independent of the choice of width; however, when
sampling the 1BZ discretely this is not the case, and
practically the goal is to find results that have a weak
dependence on this parameter. The discrepancy in the
shaded bands should be viewed as an uncertainty in the
calculation. The constraints turn up on the left-hand side
because of the band gap, and on the right-hand side because
of the finite number of bands used in the calculation. All
bands for which
E
E
F
<
4
eV, where
E
F
is the valence
band maximum, were included; see Appendix
A2
for more
details.
B. Scattering
We now consider DM-electron scattering in ZrTe
5
for
the two benchmark models
standard SI and standard SD
interactions
shown in Eq.
(1)
. Specifically, we consider a
light mediator for the SI model and a heavy mediator for the
SD model. For the SI model a light mediator was chosen
due to its high sensitivity to the lowest-energy excitations,
FIG. 3. Projected constraints on DM-electron scattering cross sections at 95% C.L. (three events, no background) assuming
1
kg-year
exposure for two benchmark models shown in Eq.
(1)
. Left: SI model with a light mediator [
F
med
¼ð
q
0
=q
Þ
2
], screened with the static
dielectric shown in Fig.
6
. The red solid (dashed) curve shows the constraints with (without) the inclusion of SOC effects. For
comparison we also show projected constraints from single phonon excitations in GaAs (orange) and SiO
2
(purple) computed with
PhonoDark
[56]
(assuming an energy threshold of
ω
min
¼
20
meV), electronic excitations in an aluminum superconductor
[57]
(brown), and previous estimates for ZrTe
5
(teal)
[25]
. We also show the projected constraints combining the SOC energy levels with the
No SOC wave functions (
Partial SOC
, red, dashed) to explicitly show the effect of the spin-dependent wave functions. Stellar
constraints (gray) are taken from Ref.
[58]
and the freeze-in benchmark (orange) is taken from Ref.
[59]
. Right: SD model with a heavy
mediator (
F
med
¼
1
). Curves labeled
low/high
E
include transitions restricted to the low-/high-
E
regions discussed in Sec.
III
.
HSIAO-YI CHEN
et al.
PHYS. REV. D
106,
015024 (2022)
015024-6
as well as for ease of comparison with other proposals
which commonly report constraints on this model. The SD
model was chosen to highlight the effect of spin-dependent
wave functions.
3
The results are shown in Fig.
3
and we discuss them in
detail here. Constraints computed in this work are shown
in red, with shaded bands corresponding to the uncertainty
in the calculation of the screening factor/dielectric function
from the electron width parameter, discussed previously in
Sec.
III A
.
When considering the SI model with a light mediator we
include anisotropic screening effects in the
f
e
ð
q
Þ
f
0
e
¼ð
ˆ
q
·
ε
ð
q
;
ω
Þ
·
ˆ
q
Þ
1
ð
30
Þ
factor.
ε
ð
q
;
ω
Þ
is the dielectric tensor, and this screening
factor is especially important for the sub-MeV DM masses
considered here. Since in this model the scattering rate is
dominated by events with small
q
, we approximate
ε
ð
q
;
ω
Þ
ε
ð
0
;
ω
Þ
, such that we replace the dielectric with
the anisotropic, long-wavelength dielectric function shown
in Fig.
6
. There is no such screening in the SD model:
f
e
¼
f
0
e
.
For the SI model, we find that the contributions from
transitions in the low-
E
region, discussed earlier in Sec.
III
,
dominate the scattering rate. Therefore, in the left panel of
Fig.
3
, we only show the results derived from transitions
within the low-
E
region. For the massive mediator SD
model, in the right panel of Fig.
3
, we see that the low-
E
contributions dominate at small DM masses. However,
for
m
χ
100
keV, when the high-
E
contributions at
O
ð
100
meV
Þ
become kinematically available, the high-
E
contributions are dominant. This is due to the fact that when
scattering via a heavy mediator the rate is no longer
dominated by the smallest momentum transfers. While
we did not explicitly include transitions between the low-
and high-
E
regions, we note that these are only expected to
be important for masses where the reach is comparable
between the regions, and will not affect the conclusions.
We find that, for the SI model with a light mediator, the
inclusion of SOC effects significantly alters the reach for
the whole DM mass range considered since the rate is
dominated by small energy/momentum depositions. For the
SD model with a massive mediator the SOC effects are
most prominent for low DM masses when the scattering is
probing the band structure near the band gap, which is the
most affected by SOC effects. We also see that at the lowest
masses the
Partial SOC
curve is closer to the
SOC
than
the
No SOC
lines. This shows that while the change to
the energy levels is the dominant effect when including
SOC, the spin dependence of the wave functions can give
O
ð
1
Þ
variations.
The left-hand side of all of the constraint curves are
determined by the band gap. The smallest kinematically
allowed DM mass is
m
χ
¼
6
keV for the SOC calculation
with
E
g
¼
23
.
5
meV, and
m
χ
¼
21
keV for the No SOC
with
E
g
¼
81
.
6
meV. As mentioned in Sec.
III A
, we only
consider bands up to 4 eV above the valence band
maximum. Kinematically, this means that we are only
including all contributions for
m
χ
<
MeV, and explains
why our projections stop there.
IV. CONCLUSIONS
Materials with strong spin-orbit coupling, such as ZrTe
5
,
are promising targets in which electronic excitations can be
utilized to search for sub-MeV DM. Their
O
ð
meV
Þ
band
gaps lead to sensitivity to new DM parameter space via
both absorption and scattering processes, without relying
on detecting single collective excitation modes.
However, due to the spin-orbit coupling, in these
materials the electron spin is no longer a good quantum
number, and the spin sums over electronic states cannot be
trivially reduced. This introduces interesting wrinkles in the
DM absorption and scattering rate calculations, which we
extended to account for these effects. In addition, we
updated the EXCEED-DM program
[30,31]
, which com-
putes DM-electron interaction rates from first principles, to
be compatible with this input for future study of general
targets with spin-orbit coupling.
We considered a wide range of DM models and
processes to which materials with SOC are sensitive:
absorption of vector, pseudoscalar, and scalar DM in
Sec.
III A
, and scattering via heavy and light mediators
via spin-independent and spin-dependent scattering poten-
tials in Sec.
III B
. We found that for sub-eV vector and
scalar DM absorption, ZrTe
5
is a far superior target relative
to an aluminum superconductor. We also found more
optimistic projections for SI scattering via a light mediator
than previous estimates, and computed, for the first time,
the projected constraints on an SD model with a heavy
mediator. Our projections for ZrTe
5
lay the foundation for
further first-principles studies of materials with strong spin-
orbit coupling as targets in direct-detection experiments.
ACKNOWLEDGMENTS
H.-Y. C. and M. B. were supported by the National
Science Foundation under Grant No. DMR-1750613.
A. M., T. T., K. Z., and Z. Z. were supported by the U.S.
Department of Energy, Office of Science, Office of High
Energy Physics, under Award No. DE-SC0021431, and the
Quantum Information Science Enabled Discovery
(QuantISED) for High Energy Physics (KA2401032).
3
In the SD model,
F
med
¼
1
also for a light mediator due to
thedominanceofthelongitudinalcomponent.Herewe focusonthe
heavy mediator case. To avoid perturbativity constraints on the
couplings,
g
χ
g
e
ð
4
π
Þ
2
, one needs
m
A
0
3
GeV
ð
10
37
cm
2
̄
σ
SD
e
Þ
1
=
4
for
keV
<m
χ
<
MeV.
DARK MATTER DIRECT DETECTION IN MATERIALS WITH
...
PHYS. REV. D
106,
015024 (2022)
015024-7
K. Z. was also supported by a Simons Investigator Award. Z.
Z. was also supported by the U.S. Department of Energy
under the Grant No. DE-SC0011702. The computations
presented here were conducted in the Resnick High
Performance Computing Center, a facility supported by
Resnick Sustainability Institute at the California Institute
of Technology.
APPENDIX A: NUMERICAL DETAILS
1. Density functional theory
The DFT calculations are carried out within the gener-
alized gradient approximation
[60]
using the Quantum
ESPRESSO code
[61]
with and without SOC included. We
use the ZrTe
5
experimental lattice constants
a
¼
3
.
9797
Å,
b
¼
14
.
470
Å, and
c
¼
13
.
676
Å of the orthorhombic
crystal structure
[62]
. We employ fully relativistic pseu-
dopotentials for calculations including SOC, and scalar
relativistic pseudopotentials for calculations without SOC,
in both cases generated with Pseudo Dojo
[63
65]
. In each
case, we use a 3265 eV kinetic energy cutoff on a uniform
4
×
4
×
2
Brillouin zone grid to compute the electron
density. To systematically converge the absorption and
scattering rates, for the high-
E
region we compute the
electronic wave functions with 200, 300, and 400 eV
cutoffs on
10
×
10
×
10
,
12
×
12
×
12
, and
14
×
14
×
14
k
-grids. For the low-
E
region, we compute the wave
functions with 650, 750, and 850 eV cutoffs on
8
×
8
×
8
,
9
×
9
×
9
, and
10
×
10
×
10
uniform
k
-grids
in a small reciprocal-space volume that includes the low-
energy band dispersion. The convergence of these calcu-
lations is discussed in Appendix
A2
.
The computed band structure of ZrTe
5
is presented in
Fig.
1
, where we correct the band gap with a scissor shift to
match the experimental band gap for the calculation with
SOC. The inset shows in detail the dispersion near the band
edges, highlighting the linear dispersion along the intra-
layer directions
Γ
-Y and
Γ
-Z. Note that in interlayer
directions (not shown in Fig.
1
) the dispersion is not linear
or conical. This band structure obtained by combining the
experiment lattice constant and the Perdew-Burke-
Ernzerhof exchange correlation functional is consistent
with a previous study
[66]
. While the presence of a
Dirac cone in ZrTe
5
is still under debate
[32
43]
, pursuing
FIG. 4. Convergence of the constraints on DM absorption, for the models discussed in Sec.
III A
, with respect to the
k
-point sampling
(
k
-grid) and plane-wave energy cutoff,
E
cut
. The first row includes SOC effects, while the second row does not. Absorption rates were
computed by adding the contributions from the low-
E
and high-
E
regions, and the first (second) value in the legends corresponds to the
parameter used in the low (high)
E
calculation. For example, the red dotted line corresponds to a calculation in which the low (high)
E
region was sampled on an
8
×
8
×
8
ð
10
×
10
×
10
Þ
Monkhorst-Pack grid in the 1BZ, with
E
cut
¼
650
ð
200
Þ
eV. All curves assume a
width parameter of
δ
¼
10
1
ω
.
HSIAO-YI CHEN
et al.
PHYS. REV. D
106,
015024 (2022)
015024-8