of 16
CALT-TH-2022-007
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, CA 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, CA 91125, USA
4
Department of Physics, University of California, Santa Barbara, CA 93106, USA
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 (DFT) 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.
CONTENTS
I. Introduction
1
II. DM Interaction Rate Formalism
2
A. Absorption
2
B. Scattering
4
III. Detection Rates in ZrTe
5
5
A. Absorption
5
B. Scattering
6
IV. Conclusions
8
Acknowledgments
8
A. Numerical Details
8
1. Density Functional Theory (DFT)
8
2. DM Interaction Constraint Convergence and
Dielectric Function
9
B. Generalized Self-Energies
10
C. Analytic Approximations in Dirac Materials
13
References
15
I. INTRODUCTION
Detection of dark matter (DM) through non-
gravitational interactions remains one of the main goals
of particle physics. Electronic excitations have been iden-
tified as a promising path to lead the direct detection of
DM to sub-GeV masses, a region not kinematically acces-
sible 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], ex-
citations 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 semiconduc-
tors for which
O
(meV) band gaps (as opposed to typi-
cal
O
(eV) band gaps in semiconductors and insulators)
arise as a consequence of spin-orbit coupling (SOC) ef-
fects. Targets with such small band gaps can probe DM
masses down to
O
(keV) via scattering, and
O
(meV) via
absorption, while still suppressing thermal noise. More-
over, some of these SOC materials have tunable band
structures, a property which makes them interesting can-
didates 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
op-
erator, 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 numeri-
cally within
EXCEED-DM
[30, 31], which is publicly avail-
able 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 charac-
terized 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 ad-
ditional consequence is that they have weak electromag-
netic screening, even with a small energy gap between
the valence and conduction bands, making them a desir-
able target for sub-MeV dark matter coupled to electrons
via a dark photon mediator [25]. In this work, however,
arXiv:2202.11716v1 [hep-ph] 23 Feb 2022
2
we will focus only on ZrTe
5
properties which stem from
its SOC nature, and we will not exploit any of the ones
deriving from its Dirac nature (which is still debated [32–
43]).
To illustrate the variety of DM models that an SOC
target can probe, we consider several different DM mod-
els and processes. Specifically, we will study:
Standard
spin-independent
(SI)
and
spin-
dependent (SD) scattering via vector mediators.
The fundamental 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.
Scalar, pseudoscalar, and vector DM absorption. In
this case the fundamental interaction Lagrangians
take the form
L
int
=
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 gen-
eralize the DM interaction rate formalism to account for
spin-dependent wave functions in general (spin-orbit cou-
pled, 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 DFT calculation are presented
in App. A 1, and convergence tests for the results shown
in Sec. III can be found in App. A 2.
II. DM INTERACTION RATE FORMALISM
In this section we derive the rates for transitions be-
tween electronic energy levels induced by DM absorption
and scattering. For the targets of interest here, the elec-
tronic energy levels can be labelled by a band index
i
and a momentum
k
within the first Brillouin zone (1BZ),
which we collectively indicate with an index
I
=
{
i,
k
}
.
The wave functions of the electronic states can be written
in the Bloch form as:
Ψ
I
(
x
) =
1
V
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 non-relativistic (NR) ef-
fective 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 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Π
ˆ
φ
ˆ
φ
2
)
1
= 1 +
O
(
g
2
e
) is the wave function
renormalization which we will approximate as unity in
the following. The total absorption rate per unit target
mass,
R
, is given by
R
=
ρ
φ
ρ
T
m
φ
1
n
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 GeV cm
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
d
4
Q
(
A
λ
φ
η
)
(
Π
λλ
AA
Π
λη
Π
ηλ
φA
m
2
φ
δ
ηη
+ Π
ηη
φφ
)
(
A
λ
φ
η
)
,
(6)
where the implicit sum over
λ,λ
(
η,η
) runs over the
photon (DM) polarizations, and we have introduced the
self-energies polarization components defined as:
Π
λλ

λ
μ
Π
μν

λ
ν
(7)
3
where

λ
μ
are polarization vectors. In general, the polar-
ization vectors which diagonalize this matrix are not the
typical longitudinal and transverse polarization vectors,
since mixing can occur (i.e., Π
L,T
AA
6
= 0). However one can
always find an appropriate basis to diagonalize the DM
and photon self-energies. In this basis Eq. (6) becomes
S
eff
⊃−
1
2
d
4
Q
(
A
λ
φ
η
)
(
Π
λ
AA
δ
λλ
Π
λη
Π
ηλ
φA
(
m
2
φ
+ Π
η
φφ
)
δ
ηη
)
(
A
λ
φ
η
)
,
(8)
where Π
λ
AA
and Π
η
φφ
are the eigenvalues of Π
λλ
AA
and Π
ηη
φφ
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:
Π
η
ˆ
φ
ˆ
φ
= Π
η
φφ
+
λ
Π
ηλ
φ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 free-
dom, and therefore Eq. (9) takes the form
Π
ˆ
φ
ˆ
φ
= Π
φφ
+
λ
Π
λ
Π
λ
φA
m
2
φ
Π
λ
AA
((pseudo)scalar DM)
.
(11)
As usual, the self-energies appearing in the previous
equations are computed from the sum of 1PI diagrams.
Working at one loop, there are two graph topologies that
can contribute
Q
−→
O
1
O
2
≡ −
i
̄
Π
O
1
,
O
2
(
Q
)
,
(12)
Q
−→
O
≡ −
i
̄
Π
O
(
Q
)
,
(13)
where
O
(1
,
2)
is any operator to which the external field,
A
or
φ
(dashed lines), couples. For vector external states
these operators carry Lorentz indices that are inherited
by
̄
Π and
̄
Π
.
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 diagonaliza-
tion 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, is reported in Ap-
pendix A 2.
1
For scalar and pseudoscalar DM, the lead-
ing order terms in the self-energy of the mostly DM eigen-
state 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)
1
Strictly speaking, the dielectric function is a mixed index
tensor, as evident from the defining equation,
J
i
=
σ
i
j
A
j
=
(
1
ε
i
j
)
A
j
, where
ε
i
j
=
δ
i
j
Π
i
j
2
= 1 + Π
ij
2
(see the
discussion in Appendix A of Ref. [25] for more details). With a
slight abuse of notation we define the matrix
ε
which has compo-
nents [
ε
]
ij
=
ε
i
j
.
4
2
1
0
1
2
Energy [eV]
SOC
No SOC
0
.
2
0
.
0
0
.
2
Γ
low
E
Y
Γ
Z
T
S
R
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 [44], and the No SOC band structure is shifted accordingly,
which gives a larger band gap of 81
.
6 meV.
B. Scattering
In this subsection we proceed to derive the DM scat-
tering rate with spin-dependent electronic wave func-
tions. Generalizing the formulas previously derived in
Refs. [2, 5, 10, 15, 31, 45], we can write the DM scatter-
ing rate as
Γ
I
I
=
π
8
V m
2
e
m
2
χ
d
3
q δ
(
E
I
E
I
ω
q
)
×
d
3
k
(2
π
)
3
̃
Ψ
I
(
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
transform 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
,σσ
(
q
) =
16
πm
2
χ
m
2
e
σ
e
μ
2
χe
f
e
f
0
e
F
med
(
q
q
0
)
S
ss
,σσ
,
(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, and
σ
e
is a reference cross section defined by
σ
e
μ
2
χe
64
πm
2
χ
m
2
e
ss
,σσ
|M
ss
,σσ
(
q
0
)
|
2
,
(21)
with
q
0
=
αm
e
. The total rate per unit detector mass is
then
R
=
π
σ
e
V μ
2
χe
m
χ
ρ
χ
ρ
T
I,I
d
3
q
(2
π
)
3
(
f
e
f
0
e
)
2
F
2
med
(
q
)
g
(
q
)
F
II
(
q
)
(22)
where
g
(
q
) is the velocity integral defined as
g
(
q
)
d
3
v
f
χ
(
v
) 2
πδ
(
ω
ω
q
)
,
(23)
with
f
χ
(
v
) being the DM velocity distribution in the
5
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
is defined as
F
II
(
q
)
d
3
k
(2
π
)
3
̃
Ψ
I
(
k
+
q
)
·
S
·
̃
Ψ
I
(
k
)
2
(24)
where the spin operators for the models considered in
this work are given by
S
ss
,σσ
=
δ
ss
δ
σσ
(SI)
,
1
3
i
σ
i
ss
σ
i
σσ
(SD)
,
(25)
and
σ
i
are the Pauli matrices. Given these expressions
the form factors,
F
II
, for the SI and SD models take the
form,
2
F
II
=
|T
II
|
2
(SI)
,
1
3
T
II
·
T
II
(SD)
,
(26)
where we have defined the DM model independent tran-
sition form factors,
T
II
and
T
II
, as
T
II
=
d
3
k
(2
π
)
3
̃
Ψ
I
(
k
+
q
)
·
̃
Ψ
I
(
k
)
,
(27)
T
II
=
d
3
k
(2
π
)
3
̃
Ψ
I
(
k
+
q
)
·
σ
·
̃
Ψ
I
(
k
)
.
(28)
III. DETECTION RATES IN
ZrTe
5
We will now apply the formalism developed in the pre-
vious 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 A 1. 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 gener-
ally only occur in small volumes within the 1BZ. To ac-
count for this we sample the 1BZ with a higher
k
-point
density in regions corresponding to the low energy band
2
The absence of the overall factor of two, relative to the SI
rate formula given in Ref. [31], can be understood from the sum
over the states. If the wave functions are spin independent then
IF
IF
ss
, where
s
(
s
) indexes the initial (final) electron
spin state. These spin sums contribute the extra factor of two,
bringing Eq. (22) and the rate formula in Ref. [31] into agreement.
structure. For ZrTe
5
this occurs near the Γ point, and we
split the phase space in to 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 low-
est 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 media-
tor. The high
E
region includes all the bands outside the
low
E
region, sampled with a standard Monkhorst-Pack
uniform grid. The DM 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 re-
gions, and it will be clear when one dominates the other.
We compute the Bloch wave functions, Eq. (3), in
both regions within the framework of DFT with
Quan-
tum ESPRESSO
[47–49]; details can be found in Ap-
pendix A 1. The DM 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 sub-
sections, we will show the projected constraints from
three different calculations. The curves labelled “SOC”
are computed with the inclusion of SOC effects, the
curves labelled “No SOC” do not include any SOC ef-
fects, and those labelled “Partial SOC” are a combina-
tion of the SOC and No SOC calculations, obtained using
the energy levels computed with SOC, and the wave func-
tions without SOC. While the Partial SOC results are not
a consistent calculation, they aid in understanding how
much of the difference 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. Gen-
erally 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.
Lastly, we note that previous works [24–26] have de-
rived excitation rates analytically by exploiting the pu-
tative Dirac nature of ZrTe
5
. While a direct comparison
to assess the validity of the analytic approximations is
dubious since we do not observe a conical band structure
(see Appendix A 1), previous estimates from Ref. [25] are
shown in Figs. (2, 3). In App. 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
6
10
3
10
2
10
1
1
10
10
2
m
φ
[eV]
10
6
10
7
10
8
10
9
10
10
d
φee
RG
Fifth Force
Al
SC
Si
Ge
ZrTe
5
Scalar DM
SOC
Partial SOC
No SOC
10
3
10
2
10
1
1
10
10
2
m
φ
[eV]
10
13
10
12
10
11
10
10
10
9
g
aee
WD
KSVZ
DFSZ
Al
SC
Ge
Si
ZrTe
5
Pseudoscalar DM
10
3
10
2
10
1
1
10
10
2
m
φ
[eV]
10
16
10
15
10
14
10
13
κ
XENON10
/
100
Ge
Si
Al
SC
ZrTe
5
ZrTe
5
[25]
SiO
2
GaAs
Vector DM
FIG. 2. Comparison of projected 95% C.L. reach (3 events, no background) assuming one 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 the ones 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 parameterizations 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 XENON10/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].
commonly used notation,
g
e
=
4
πm
e
M
Pl
d
φee
(scalar)
g
aee
(pseudoscalar)
(vector)
,
(29)
where
M
Pl
= 1
.
22
×
10
19
GeV is the Planck mass.
For all 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 su-
perconducting aluminum target, another target mate-
rial 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 width parame-
terizations,
δ
[10
1
.
5
,
10
0
.
5
]
ω
. In theory, the absorp-
tion 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 calcula-
tion. All bands for which
E
E
F
<
4 eV, where
E
F
is the
valence band maximum were included; see Appendix A 2
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 en-
ergy excitations, 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
3
In the SD model,
F
med
= 1 also for a light mediator due to
the dominance of longitudinal component. Here we focus on the
heavy mediator case. To avoid perturbativity constraints on the
couplings,
g
χ
g
e
.
(4
π
)
2
, one needs
m
A
.
3 GeV
(
10
37
cm
2
σ
SD
e
)
1
/
4
for keV
< m
χ
<
MeV.
7
1
10
10
2
10
3
m
χ
[keV]
10
43
10
42
10
41
10
40
10
39
10
38
10
37
10
36
σ
SI
e
[cm
2
]
Freeze
In
Stellar
ZrTe
5
[25]
Al
SC
ZrTe
5
SiO
2
GaAs
1
10
10
2
10
3
m
χ
[keV]
10
39
10
38
10
37
10
36
10
35
10
34
10
33
σ
SD
e
[cm
2
]
low
E
high
E
SOC
Partial SOC
No SOC
FIG. 3. Projected constraints on DM-electron scattering cross sections at the 95% C.L. (three events, no background) assuming
one 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 labelled “low/high
E
” include transitions
restricted to the low/high
E
regions discussed in Sec. III.
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 uncer-
tainty 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
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 di-
electric function shown in Fig. 6.
For the SI model, we find that the contribution 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 me-
diator 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
contri-
butions 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 scat-
tering 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 depen-
dence of the wave functions can give
O
(1) variations.
The left hand side of all the constraint curves are de-
termined 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
8
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 excita-
tions can be utilized to search for sub-MeV DM. Their
O
(meV) band gaps lead to sensitivity to new DM param-
eter space via both absorption and scattering processes,
without relying on detecting single collective excitation
modes.
However, due to the spin-orbit coupling, in these ma-
terials 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 ad-
dition, we updated the
EXCEED-DM
program [30, 31],
which computes 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 pro-
cesses to which materials with SOC are sensitive: absorp-
tion of vector, pseudoscalar, and scalar DM in Sec. III A,
and scattering via heavy and light mediators via spin-
independent and spin-dependent scattering potentials 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 media-
tor 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 foun-
dation 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 Sci-
ence Foundation under Grant No. DMR-1750613. A.M.,
T.T., K.Z. and Z.Z. were supported by the U.S. Depart-
ment of Energy, Office of Science, Office of High En-
ergy Physics, under Award No. DE-SC0021431, and the
Quantum Information Science Enabled Discovery (Quan-
tISED) for High Energy Physics (KA2401032). 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 DE-SC0011702. The computations presented
here were conducted in the Resnick High Performance
Computing Center, a facility supported by Resnick Sus-
tainability Institute at the California Institute of Tech-
nology.
Appendix A: Numerical Details
1. Density Functional Theory (DFT)
The DFT calculations are carried out within the generalized gradient approximation (GGA) [60] using the
Quantum
Espresso
code [61] with and without spin-orbit coupling (SOC) included. We use the ZrTe
5
experimental lattice
constants,
a
= 3
.
9797
̊
A,
b
= 14
.
470
̊
A, and
c
= 13
.
676
̊
A of the orthorhombic crystal structure [62]. We employ fully
relativistic pseudopotentials 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 (BZ) 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 calculations is
discussed in Appendix A 2.
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 intralayer 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 (PBE) 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 more
extensive tests of crystal structure and DFT functionals, or carrying out beyond-DFT band structure calculations, is
beyond the scope of this work.
9
10
1
1
10
m
φ
[eV]
10
5
10
6
10
7
10
8
10
9
d
φee
Scalar DM
k
grid
(8
×
8
×
8
,
10
×
10
×
10)
(9
×
9
×
9
,
12
×
12
×
12)
(10
×
10
×
10
,
14
×
14
×
14)
10
1
1
10
m
φ
[eV]
10
13
10
12
10
11
10
10
10
9
g
aee
Pseudoscalar DM
E
cut
(650 eV
,
200 eV)
(750 eV
,
300 eV)
(850 eV
,
400 eV)
10
1
1
10
m
φ
[eV]
10
17
10
16
10
15
κ
SOC
Vector DM
10
1
1
10
m
φ
[eV]
10
5
10
6
10
7
10
8
10
9
d
φee
10
1
1
10
m
φ
[eV]
10
13
10
12
10
11
10
10
10
9
g
aee
10
1
1
10
m
φ
[eV]
10
17
10
16
10
15
κ
No SOC
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
ω
.
2. DM Interaction Constraint Convergence and Dielectric Function
In this appendix we will discuss some details of the DM scattering and absorption rate calculations, as well as
the long wavelength, anisotropic dielectric function,
ε
(0
). Since the main focus of this paper is the effect of SOC,
only the electronic wave functions near the Fermi surface are needed. This is because, in ZrTe
5
, SOC effects are
approximately
O
(10 meV), and therefore a very small perturbation for states
>
eV away from the Fermi surface. We
are therefore safely within the “valence to conduction” regime, discussed in more detail in Ref. [31], and do not need to
study deeper, core electronic levels, or larger energy states where the electrons are close to free. DFT is the preferred
tool for studying these transitions, and the two main convergence parameters are the number of
k
-points in the 1BZ
sampling, and the plane wave expansion cutoff,
E
cut
. In both the low
E
and high
E
regions we sample
k
points
uniformly with a Monkhorst-Pack grid. The only difference is that the low
E
points are scaled by 1
/
5 relative to the
high
E
region. Convergence of the DM absorption and scattering constraints with respect to the
k
point sampling
and
E
cut
parameters are shown in Fig. 4 and Fig. 5 respectively. The constraints in the main text are identical to the
most converged constraints shown in Figs. (4, 5). Generally we see faster convergence with respect to
E
cut
than the
k
point density, and slightly faster convergence for the DFT calculation which omits SOC effects than those which
include them. We also note that all-electron reconstruction effects were omitted here since we are focusing on very
small DM masses, and therefore kinematically limited to small
q
transitions. However these effects could be important
for studies of DM scattering in ZrTe
5
at higher masses, or for larger experimental thresholds.
The dielectric function in the long wavelength limit is shown in Fig. 6, was used as an intermediate to compute a
few different constraints. Specifically it was used to screen the SI scattering rate, and it can be shown that the vector
DM absorption rate, as well as the pseudoscalar DM absorption rate when wave functions are spin independent, can
be related to the dielectric function. Moreover this calculation serves as a useful benchmark to compare future DFT
calculations.