# Gravitational Focusing of Wave Dark Matter

Hyungjin Kim<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany  Alessandro Lenoci
<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron DESY, Notkestr. 85,
22607 Hamburg, Germany

###### Abstract

A massive astrophysical object deforms the local distribution of dark matter,
resulting in a local overdensity of dark matter. This phenomenon is often
referred to as gravitational focusing. In the solar system, the gravitational
focusing due to the Sun induces modulations of dark matter signals on
terrestrial experiments. We consider the gravitational focusing of light
bosonic dark matter with a mass of less than about 10 eV. The wave nature of
such dark matter candidates leads to unique signatures in the local
overdensity and in the spectrum, both of which can be experimentally relevant.
We provide a formalism that captures both the gravitational focusing and the
stochasticity of wave dark matter, paying particular attention to the
similarity and difference to particle dark matter. Distinctive patterns in the
density contrast and spectrum are observed when the de Broglie wavelength of
dark matter becomes comparable or less than the size of the system and/or when
the velocity dispersion of dark matter is sufficiently small. While
gravitational focusing effects generally remain at a few percent level for a
relaxed halo dark matter component, they could be much larger for dark matter
substructures. With a few well-motivated dark matter substructures, we
investigate how each substructure responds to the gravitational potential of
the Sun. The limit at which wave dark matter behaves similar to particle dark
matter is also discussed.

††preprint: DESY-21-217

##  I Introduction

The evidence for dark matter (DM) has been accumulating for decades. Since the
evidence is based on the gravitational interaction, it still remains a mystery
how dark matter communicates with the standard model. While there exists a
plethora of possibilities, a particular class of DM candidates, namely weakly
interacting massive particle, has given particular attention due to its
predictive power and theoretical motivations. However, it has avoided intense
scrutiny of terrestrial DM searches as well as cosmological and astrophysical
probes.

A light bosonic particle with a mass smaller than $10\,{\rm eV}$ provides a
compelling alternative dark matter candidate. Such a DM candidate arises from
various beyond the standard models with different theoretical and
phenomenological motivations. One specific example is the QCD axion [1, 2, 3],
which is originally proposed as a solution to the strong CP problem in the
standard model [4, 5, 6, 7, 8, 9, 10]. Not only from a solution to the strong
CP problem, it also arises from dynamical and anthropic solutions to the
electroweak hierarchy problem [11, 12, 13, 14, 15, 16]. Moreover, a large
number of light bosonic states may arise from higher form fields on compact
extra dimensions in some ultraviolet theories [17, 18], and one of them can be
dark matter in the present universe.

A light bosonic DM candidate behaves like a classical wave rather than an
individual particle. With a typical value of DM mass density in a collapsed
halo, the number of DM particles within a volume of de Broglie wavelength
takes a gigantic value. This allows us to treat DM as a collection of waves,
and hence, we call such candidates wave dark matter candidates. The wave
nature makes them qualitatively different from usual particle DM candidates,
leaving unique cosmological and astrophysical signatures throughout the entire
history of the universe. See a recent review [19] for an overview.

As it appears ubiquitously in many beyond the standard models, many efforts
have been put forward to find such wave dark matter in Milky Way for past
decades (see e.g. [20] for a review on axion and axion-like particle
searches). When it comes to terrestrial searches of dark matter, one essential
input is the dark matter local distribution. A common starting point is to
assume that the dark matter forms an isothermal halo and its velocity
distribution is given by the Maxwell-Boltzmann distribution [21, 22]. Such a
model is referred to as the standard halo model, and as its name implies, it
has been the standard assumption for various analyses of terrestrial dark
matter searches, not only for the wave dark matter but also for particle dark
matter. While the standard halo model is simple and is reproduced reasonably
well from N-body studies [23, 24], it ignores several important aspects of our
current understanding of the local dark matter distribution, which might be
important to determine a precise map of local dark matter.

Firstly, it does not account for the potential existence of DM substructures
near the solar neighborhood. With precision astrometric data, recent analyses
have identified a stellar population near the solar neighborhood with a
distinct kinematic structure [25, 26]. This stellar substructure is called
Gaia-Enceladus or Gaia-Sausage due to its elongated structure in the velocity
space. It is inferred that this stellar substructure was created through a
recent merger with a satellite galaxy of mass $\sim 10^{10}M_{\odot}$ about
$8\textrm{--}10$ Gyr ago. The same merger event is expected to deposit DM into
the Milky Way. Numerical studies showed that accreted dark matter from the
same merger event has a kinematic structure very similar to its stellar
counterpart [27, 28]. A nontrivial fraction of dark matter in the inner halo
is therefore expected to be in form of substructures [29, 30, 31]. In addition
to Gaia-Enceladus, there could be other substructures with completely
different kinematic properties. Examples include streams with
characteristically small velocity dispersion, and a dark disk which may arise
from the merger events that are responsible for the formation of the thick
stellar disk in the Milky Way [32].

Secondly, it does not account for the deformation of local dark matter
distribution due to the gravitational potential of the Sun. As we are bound to
the solar system, all the obervations are eventually affected by the presence
of the Sun. Specifically, the gravitational potential of the Sun deflects the
trajectory of DM, gravitationally focusing DM particles in the solar system.
This leaves a tail of local overdensity in the opposite direction to the
direction of the solar system in the Milky Way. It is therefore important to
examine how the dark matter profile near the solar system is modified by the
Sun, especially along the orbital trajectory of the Earth. This gravitational
focusing of DM has been studied, especially in the context of weakly
interacting massive particle DM searches [33, 34, 35, 36, 37, 38]. While it
was found that the event rate due to the local overdensity remains only a few
percent level in the standard halo model, it can lead to an interesting
pattern of DM flux in the sky map [36, 37] as well as modifying the DM annual
modulation signal [38].

The purpose of this work is to extend the discussion to the gravitational
focusing of wave DM and to examine how wave dark matter substructures are
focused in the solar system. We consider a spinless, light, bosonic dark
matter candidate, minimally coupled to the gravity, with a mass smaller than
$10\,\,{\rm eV}$. We first develop a formalism which allows us to study the
gravitational focusing and the stochasticity of wave dark matter
simultaneously. Unlike particle DM, systems of wave DM include an additional
spatial scale characterized by de Broglie wavelength, and much of wave
properties of the problem is controlled by the de Broglie wavelength; for
instance, the overdensity becomes constant when the radial distance from an
astrophysical source becomes smaller than de Broglie wavelength and the width
of density wake field is controlled by de Broglie wavelength. All of these
features are distinct from particle dark matter and they are commonly observed
in the other context of wave dark matter studies at different scales; for
instance the same wave nature flattens the inner profile of DM halo, and
hence, may resolve small-scale issues in the cold dark matter paradigm [39].
In addition, this wave nature reduces the dynamical friction, friction induced
by the density wake field on astrophysical objects [40, 41, 42]. We discuss
below how this wave nature manifests itself in the solar system and also
discuss possible implications on terrestrial experiments.

The paper is organized as follows. In Section II, we briefly review the
gravitational focusing of particle DM for the purpose of later comparison with
the wave DM. We compute the local density contrast for the particle DM and
discuss its behavior as a function of parameters of the system, such as DM
mean velocity and dispersion. In Section III, we investigate the gravitational
response of wave DM. The investigation is done in several steps. We first
begin with an action of a minimally coupled spinless bosonic particle, and
expand the field in terms of creation and annihilation operator with a mode
function. Then, we solve the Schrödinger equation for the mode function in
Section III.1, which encodes an information on the gravitational potential and
spatial distribution of the wave DM. We then determine statistical properties
of the DM field by specifying the density operator of the system, which is
discussed in Section III.2. Specification of density operator allows us to
compute the ensemble average of arbitrary operators constructed from the
field. This procedure together with the mode function obtained from the
Schrödinger equation enables us to consider the gravitational focusing and the
stochasticity of the field at the same time. After discussing ensemble
averages of the field, we compare the wave gravitational focusing with
particle dark matter focusing, contrasting the similarities and differences
between them. In Section IV, we apply the formalism discussed in the previous
section to the DM substructures in the solar system as well as relaxed halo
dark matter to see how each substructure responds to the gravitational
potential of the Sun. We discuss the effect of coordinate transformation as
well as the classical limit of wave DM in Section V and conclude in Section
VI. Throughout the work, we use natural unit $c=\hbar=1$ except for the
section we discuss the classical limit of wave DM.

##  II Particle focusing

We review the gravitational focusing of particle dark matter. We follow [36,
37]. The results summarized in this section are to compare the particle DM
result with that of wave DM. Readers who are already familiar with the
gravitational focusing of particle DM may continue to Section III.

Figure 1: We show the density contrast of the particle dark matter. We choose
${\bf v}_{\rm dm}=(-240,0,0)$ km/s, $v_{\rm dm}/\sigma=\sqrt{2}$ and the
source mass to be the solar mass.

Let us consider the dark matter phase space distribution, $f(t,{\bf x},{\bf
p})$. In the collisionless system, the phase space density is conserved, i.e.
$df/dt=0$. According to Liouville’s theorem, the phase space volume should be
conserved along the trajectory given by the Hamiltonian of the system.

The dark matter local overdensity can be easily estimated by Liouville’s
theorem. Consider particles confined in an infinitesimal phase space volume,
$d^{3}v_{0}d^{3}x_{0}$, asymptotically far away from a source of mass $M$. For
simplicity, let us assume a constant phase space density over this
infinitesimal volume. Due to the gravitational potential, particles are
accelerated and their speed at a radial distance $r$ becomes
$v^{2}(r)=v_{0}^{2}+2GM/r$, where $v_{0}$ is the initial speed far away from
the source $M$. According to Liouville’s theorem, the phase space volume
remains constant, $d^{3}x_{0}\,d^{3}v_{0}=d^{3}x\,d^{3}v$, where
$d^{3}x\,d^{3}v$ is the phase space volume of particles at $r$. This leads to
$d^{3}x=[v_{0}/v(r)]d^{3}x_{0}$. We find that the particle mass density at $r$
is

$\displaystyle\rho(r)\approx\rho_{0}[1+v_{e}^{2}(r)/v_{0}^{2}]^{1/2},$ (1)

where $v_{e}^{2}(r)=2GM/r$ is the escape velocity and $\rho_{0}$ is the mass
density far away from the source $M$.

More quantitative analysis is possible. Suppose that $f({\bf v}_{0})$ is the
phase space distribution at $r\to\infty$, normalized as $\int
d^{3}v_{0}\,f({\bf v}_{0})=1$. The differential particle mass in an
infinitesimal phase space volume is $mdN=\rho_{0}f({\bf
v}_{0})d^{3}x_{0}d^{3}v_{0}=\rho_{0}f({\bf v}_{0})d^{3}xd^{3}v$. The mass
density at the position ${\bf x}$ relative to its asymptotic value is given by

$\displaystyle\frac{\rho({\bf x})}{\rho_{0}}=1+\delta$ $\displaystyle=$
$\displaystyle\int d^{3}v\,f({\bf v}_{0})$ (2) $\displaystyle=$
$\displaystyle\int d^{3}v_{0}\,d({\bf v}_{0},{\bf x})f({\bf v}_{0}).$

where $\delta$ is the density contrast. The Jacobian $d({\bf v}_{0},{\bf x})$
is [36]

$d({\bf v}_{0},{\bf
x})=\frac{1}{2}\left[\left(1+\frac{2v_{e}^{2}/v_{0}^{2}}{1-\hat{v}_{0}\cdot\hat{x}}\right)^{\frac{1}{2}}\\!\\!+\left(1+\frac{2v_{e}^{2}/v_{0}^{2}}{1-\hat{v}_{0}\cdot\hat{x}}\right)^{-\frac{1}{2}}\right].$

For a given velocity distribution, the above integral can be numerically
performed. One can directly compute the first integral by using an explicit
expression of ${\bf v}_{0}={\bf v}_{0}({\bf v},{\bf x})$ [37].

For the further discussion, we consider the Maxwell-Boltzmann distribution
function $f({\bf v})$,

$\displaystyle f({\bf
v})=\frac{1}{(2\pi\sigma^{2})^{3/2}}\exp\left[-\frac{({\bf v}-{\bf v}_{\rm
dm})^{2}}{2\sigma^{2}}\right].$ (3)

The velocity distribution is characterized by the mean velocity is
$\langle{\bf v}\rangle={\bf v}_{\rm dm}$ and the velocity dispersion $\sigma$.
With the Maxwell-Boltzmann distribution, the density contrast is completely
described by three parameters,

$v_{\rm dm}/\sigma,\qquad r/\bar{r},\qquad\hat{x}\cdot\hat{v}_{\rm
dm}(\equiv\mu),$

where

$\bar{r}=\frac{GM}{\sigma^{2}}=0.03\;{\rm
AU}\;\left(\frac{M}{M_{\odot}}\right)\left(\frac{240\;{\rm
km/sec}}{\sqrt{2}\sigma}\right)^{2}$

is the radius of gravitational influence. One may choose $v_{e}(r)/\sigma$
instead of $r/\bar{r}$. We present the density contrast of particle dark
matter in Figure 1 with ${\bf v}_{\rm dm}=(-240,\,0,\,0)$ km/sec and $v_{\rm
dm}/\sigma=\sqrt{2}$ for the purpose of demonstration. The density contrast is
invariant under the rotation with respect to $\hat{v}_{\rm dm}$-axis due to an
azimuthal symmetry.

A characteristic density wake field is formed along $\hat{v}_{\rm dm}$. Two
variables might be useful to inspect the density wake field,

$\displaystyle\delta_{\rm ds}=\delta|_{\mu=1},\qquad\delta_{\rm
avg}=\int\frac{d\Omega_{\hat{x}}}{4\pi}\delta=\int\frac{d\mu}{2}\delta,$ (4)

where $\delta_{\rm ds}$ is the density contrast along the downstream
($\hat{x}\cdot\hat{v}_{\rm dm}=1$) and $\delta_{\rm avg}$ is the density
contrast averaged over a sphere of radius $r$. While above variables are a
function of two parameters, $v_{\rm dm}/\sigma$ and $r/\bar{r}$, they
practically depend only on one specific combination of these two. For the
downstream density contrast, we find

$\displaystyle 1+\delta_{\rm
ds}\approx\Big{(}1+\frac{v_{e}^{2}}{\sigma^{2}}\Big{)}^{\frac{1}{2}}=\Big{(}1+\frac{2\bar{r}}{r}\Big{)}^{\frac{1}{2}}.$
(5)

while, for the averaged density contrast, we find

$\displaystyle 1+\delta_{\rm avg}$ $\displaystyle=$ $\displaystyle\int
d^{3}v_{0}\,f({\bf
v}_{0})\big{(}1+\frac{v_{e}^{2}}{v_{0}^{2}}\big{)}^{\frac{1}{2}}$ (6)
$\displaystyle\approx$
$\displaystyle\Big{(}1+\frac{v_{e}^{2}}{\sigma^{2}+v_{\rm
dm}^{2}}\Big{)}^{\frac{1}{2}}.$

We check these expressions numerically in Appendix A. The numerical results
are within a factor of two for wide range of parameters. Note also that the
ratio $\delta_{\rm avg}/\delta_{\rm ds}$ represents the angular scale over
which the density contrast takes nonvanishing value. For $\sigma>v_{e}$ and
$v_{\rm dm}>\sigma$, this angular scale is approximately given by
$\Delta\mu\sim(\sigma/v_{\rm dm})^{2}$.

##  III Wave focusing

We consider the wave dark matter in this section. For concreteness, we
consider a spinless bosonic dark matter with a mass $m$, minimally coupled to
the gravity. We approach the problem in several steps. We begin with the
action and expand the field in terms of creation and annihilation operators
with a mode function. A mode function for the system is obtained by solving
the equation of motion, which is discussed in Section III.1. The mode function
obtained in this way contains information about spatial distribution of the
wave near the source $M$. Statistical properties of the field are determined
by specifying the underlying density operator of the system, which is
discussed in Section III.2.

The action of the system is

$\displaystyle S=\int
d^{4}x\,\sqrt{-g}\left[\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-\frac{1}{2}m^{2}\phi^{2}\right],$
(7)

where $g=\det g_{\mu\nu}$ is the determinant of the metric, $g_{\mu\nu}$ is
the metric tensor, and $m$ is the dark matter mass. The metric is given as

$\displaystyle ds^{2}=(1+2\Phi)dt^{2}-(1-2\Phi)d\vec{x}\,^{2}.$ (8)

where $d\vec{x}\,^{2}=dr^{2}+r^{2}d\Omega$ and $\Phi$ is the gravitational
potential. Since we are interested in non-relativistic dark matter, we expand
the field as

$\displaystyle\hat{\phi}(x)=\sum_{i}\frac{1}{\sqrt{2mV}}\Big{[}\hat{a}_{i}\psi_{i}(x)e^{-imt}+\hat{a}_{i}^{\dagger}\psi_{i}^{*}(x)e^{imt}\Big{]},$
(9)

where $\hat{a}_{i}$ and $\hat{a}_{i}^{\dagger}$ are an annihilation and
creation operator, satisfying the canonical commutation relation
$[\hat{a}_{i},\hat{a}_{j}^{\dagger}]=\delta_{ij}$. We have expanded the field
in a finite box of volume $V$. The index $i$ denotes all possible quantum
numbers. We solve the equation of motion to determine the mode function
$\psi_{i}(x)$ and specify the density operator in the following sections as
advertised.

###  III.1 Wave function

We first solve the equation of motion for the mode function. The wave function
is given by the solution to the Klein-Gordon equation,
$(\square+m^{2})\phi=0$. In the nonrelativistic liimt, we expand the field as
in Eq. (9). In the limit $|\ddot{\psi}/\dot{\psi}|,\,|\dot{\psi}/\psi|\ll m$
and $|\dot{\Phi}/\Phi|\ll m$, the Klein-Gordon equation is reduced to the
Schrödinger equation,

$\displaystyle
i\partial_{t}\psi=\Big{(}-\frac{\nabla^{2}}{2m}+m\Phi\Big{)}\psi.$ (10)

The gravitational potential $\Phi$ is determined by the Poisson equation,
$\nabla^{2}\Phi=4\pi G(\rho_{M}+\rho_{\rm dm})$, where $\rho_{M}$ and
$\rho_{\rm dm}$ are the mass density of the source and dark matter,
respectively. Ignoring the self-gravity, the gravitational potential is
approximated as $\Phi(r)=-GM/r$. In this case, solving the Schrödinger
equation Eq. (10) is identical to solving the Coulomb scattering problem of an
electron under the central $1/r$ potential generated by the proton. With
$\psi_{i}(x)\to e^{-i(k_{i}^{2}/2m)t}\psi_{i}({\bf x})$, the wave function is
given as [43, 40]

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\psi_{{\bf k}}({\bf x})\\!=\\!e^{i{\bf
k}\cdot{\bf
x}}\Gamma(1-i\beta)e^{\pi\beta/2}M\big{[}i\beta,1,ikr(1-\hat{k}\cdot\hat{x})\big{]}.$
(11)

where $\beta=GMm^{2}/k$ and $M(a,b,c)$ is the confluent hypergeometric
function. Each mode is characterized by the wavenumber ${\bf k}$. For $kr\gg
1$ and $\mu=\hat{k}\cdot\hat{x}<0$, the wave function is approximated a plane
wave $\psi_{{\bf k}}({\bf x})\sim e^{i{\bf k}\cdot{\bf x}}$.

###  III.2 Statistical properties

The discussion in the previous section shows how a plane wave of wavenumber
${\bf k}$ responds to the gravitational potential of the source $M$. A
realistic dark matter halo is a superposition of different modes, and so is
the field $\hat{\phi}$. To determine the statistical properties of
$\hat{\phi}$, we must specify the density operator of the system.

####  III.2.1 Density operator

Before we discuss the density operator for the field, let us begin with a
simple harmonic oscillator in quantum mechanics as a toy model. The
statistical properties of simple harmonic oscillator are described by the
density operator $\hat{\rho}$. If the system is in the thermal equilibrium,
the density operator is given by $\hat{\rho}=e^{-\beta\hat{H}}/\,{\rm
tr}(e^{-\beta\hat{H}})$ with a Hamiltonian of the system $\hat{H}=\omega
a^{\dagger}a$. Alternatively, the density matrix can be written as

$\displaystyle\hat{\rho}=\frac{1}{1+\langle
n\rangle}\sum_{n}\left(\frac{\langle n\rangle}{\langle
n\rangle+1}\right)^{n}|n\rangle\langle n|,$ (12)

where $|n\rangle$ is an eigenstate of the number operator
$\hat{N}=a^{\dagger}a$. The mean occupation number is $\langle
n\rangle=(e^{\beta\omega}-1)^{-1}$. While the above form is obtained for the
thermal system, it applies for much broader systems as long as the statistical
properties of the excitation are chaotic [44]. More specifically, the above
density matrix maximizes the entropy of the system for a fixed mean occupation
number $\langle n\rangle$. Without assuming thermal equilibrium, we use the
above form of density matrix for the following discussion.

Since we are interested in a system with a large occupation number, it is
convenient to express the above density matrix in the coherent state
representation. The coherent state is defined as
$|\alpha\rangle=e^{-|\alpha|^{2}/2}\sum\frac{\alpha^{n}}{\sqrt{n!}}|n\rangle$
with a complex number $\alpha$. It is an eigenstate of the annihilation
operator, $\hat{a}|\alpha\rangle=\alpha|\alpha\rangle$. With the coherent
state, the density operator can be expressed [45, 46],

$\displaystyle\hat{\rho}$ $\displaystyle=$ $\displaystyle\int
d^{2}\alpha\,P(\alpha)|\alpha\rangle\langle\alpha|.$ (13)

where $P(\alpha)$ is a quasi-probability distribution, defined as

$\displaystyle P(\alpha)=\frac{1}{\pi\langle
n\rangle}\exp\left[-\frac{|\alpha|^{2}}{\langle n\rangle}\right].$ (14)

For the density matrix Eq. (12), the quasi-probability distribution has
desired properties to be interpreted as a probability distribution for a
complex number $\alpha$; it is positive, $P(\alpha)\geq 0$, and it is
integrated up to unity, $\int d^{2}\alpha\,P(\alpha)=1$. Moreover, the quasi-
probability distribution can be factorized into two probability distributions,
one for the modulus $|\alpha|$ and one for the phase $\theta={\rm
arg}(\alpha)$,

$P(\alpha)=P(|\alpha|)P(\theta)$

where the modulus $|\alpha|$ follows the Rayleigh distribution,
$P(|\alpha|)=2(|\alpha|/\langle n\rangle)\exp[-|\alpha|^{2}/\langle
n\rangle]$, while the phase $\theta$ follows the uniform distribution,
$P(\theta)=1/2\pi$. The expectation value of any operator $\hat{\cal
O}(\hat{a},\hat{a}^{\dagger})$ is given by $\langle\hat{\cal O}\rangle=\,{\rm
tr}(\hat{\cal O}\hat{\rho})=\int(d^{2}\alpha/\pi)\langle\alpha|\hat{\cal
O}\hat{\rho}|\alpha\rangle$. Since the coherent state is the eigenstate of
annihilation operator, we may replace $\hat{\cal
O}(\hat{a},\hat{a}^{\dagger})\to{\cal O}(\alpha,\alpha^{*})$ and compute the
expectation value with the probability distribution $P(\alpha)$,

$\displaystyle\langle\hat{\cal O}(\hat{a},\hat{a}^{\dagger})\rangle\approx\int
d^{2}\alpha\,P(\alpha){\cal O}(\alpha,\alpha^{*}).$ (15)

The only subtlety is that the commutation of ladder operators does not vanish,
leading to additional terms in the above expression. However, those additional
contributions are suppressed by powers of the mean occupation number, and
therefore, they are subleading as long as the occupation number is large. The
error of the above approximation is ${\cal O}(1/\langle n\rangle)$.

Motivated by this, we take the same form of density matrix for the system of
the DM field $\hat{\phi}$. Since we decomposed the field as a linear
combination of different modes, the density operator is now given as product
of Eq. (13) for each mode $i$,

$\displaystyle\hat{\rho}=\bigg{[}\prod_{i}\int
d^{2}\alpha_{i}\,P(\alpha_{i})\bigg{]}|\\{\alpha_{i}\\}\rangle\langle\\{\alpha_{i}\\}|,$
(16)

where $P(\alpha_{i})$ is the probability distribution for each $\alpha_{i}$.
The form of $P(\alpha_{i})$ is the same as in Eq. (14). The coherent state
here is defined as
$\hat{a}_{j}|\\{\alpha_{i}\\}\rangle=\alpha_{j}|\\{\alpha_{i}\\}\rangle.$

The density operator completely determines the statistical properties of the
scalar field. We can compute the expectation value of any operators
constructed from $\hat{\phi}$. Similarly to the quantum mechanical harmonic
oscillator, one may replace $\hat{a}_{i}\to\alpha_{i}$, or $\hat{\cal
O}(\hat{a}_{i},\hat{a}_{i}^{\dagger})\to\hat{\cal
O}(\alpha_{i},\alpha^{*}_{i})$, and treat each $\alpha_{i}$ as a random
variable following the probability distribution $P(\alpha_{i})$. As in the
quantum mechanics example, there are additional terms arising from
nonvanishing commutation relation between creation and annihilation operators,
but such contributions are again suppressed by powers of the mean occupation
number. For instance, we find

$\displaystyle\langle a_{i}^{\dagger}a_{j}\rangle$ $\displaystyle=$
$\displaystyle\delta_{ij}\langle n_{i}\rangle,$ (17) $\displaystyle\langle
a_{j}a_{i}^{\dagger}\rangle$ $\displaystyle=$
$\displaystyle\delta_{ij}(\langle n_{i}\rangle+1)\approx\delta_{ij}\langle
n_{i}\rangle,$ (18)

where we see the quantum fluctuation is suppressed by $1/\langle
n_{i}\rangle$. This allows us to write the quantum field as

$\hat{\phi}(t,{\bf x})\to\phi(t,{\bf
x})=\sum_{i}\frac{1}{\sqrt{2mV}}\Big{[}\alpha_{i}\psi_{i}({\bf
x})e^{-i\omega_{i}t}+{\rm h.c.}\Big{]},$

where $\omega_{i}=m+k_{i}^{2}/2m$. We treat each $\alpha_{i}$ as a random
variable, following the probability distribution Eq. (14), while ignoring
possible quantum fluctuation suppressed by $1/\langle n\rangle$. In the plane
wave limit, $\psi_{{\bf k}}({\bf x})=e^{i{\bf k}\cdot{\bf x}}$, this approach
easily reproduces the statistical properties of the scalar field $\hat{\phi}$
obtained in different approaches [47, 48]. Including these stochastic
properties is important for data analysis of wave DM experiments as it may
affect interpretation of data [47, 48, 49, 50, 51].

####  III.2.2 Moments and spectrum

Equipped with the density matrix and the wave function, we proceed to compute
the ensemble averages, i.e. moments, of the field $\hat{\phi}$. We are
particularly interested in the second moment of the field for two reasons. The
second moment of the field $\langle\hat{\phi}^{2}\rangle$ is directly related
to the density contrast as
$1+\delta=\langle\hat{\phi}^{2}\rangle/\phi_{0}^{2}$, where $\phi_{0}$ is the
root mean square of the field value measured at asymptotically far away from
the source. From this, we can easily compare the density contrast of wave and
particle dark matter. Secondly, this quantity is, in many cases, directly
related to the observables in the terrestrial detectors. For instance, in
axion and axion-like particle searches for ${\cal L}\supset\phi
F\widetilde{F}$,111For instance, ADMX [52, 53, 54], HAYSTAC [55, 56],
ABRACADABRA [57, 58], DM-Radio [59], CAPP [60], MADMAX [61], ORGAN [62] and so
on. a total signal power is often proportional to
$\langle\hat{\phi}^{2}\rangle$ and signal power spectral density is related to
the spectral density of $\langle\hat{\phi}^{2}\rangle$. The same holds for
scalar dark matter searches on ${\cal L}\supset\phi\bar{\psi}\psi,\,\phi FF$
via the oscillations of fundamental constants.222For instance, atomic clocks
[63, 64, 65, 66, 67, 68, 69], atomic and molecular spectroscopy searches [70,
71], and interferometry [72, 73, 74, 75, 76].

Using a generalized expression of Eq. (15), we obtain the variance of the
field as

$\displaystyle\left\langle\hat{\phi}^{2}(t,{\bf x})\right\rangle$
$\displaystyle=$ $\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{\bf
k}}{m}|\psi_{\bf k}({\bf x})|^{2}$ (19) $\displaystyle=$
$\displaystyle\phi_{0}^{2}\int d^{3}v\;f({\bf v})|\psi_{\bf k}({\bf x})|^{2}\
,$

where we have taken the continuum limit in the first line. Here, $\phi_{0}$ is
the mean field value at asymptotically far away from the source and $n_{{\bf
k}}$ is the occupation number, which is the continuum version of $n_{i}$
defined in the density matrix. We have also defined the velocity probability
distribution function $f({\bf v})=(m/2\pi)^{3}n({\bf k})/(m\phi_{0}^{2})$. The
asymptotic field value $\phi_{0}$ is related to the asymptotic mass density as
$\phi_{0}=\sqrt{\rho_{0}}/m$. The variance of the field is directly related to
the density contrast as

$\displaystyle 1+\delta=\frac{\langle\phi^{2}\rangle}{\phi_{0}^{2}}=\int
d^{3}v\,f({\bf v})|\psi_{\bf v}({\bf x})|^{2}.$ (20)

This expression can be compared with the corresponding expression for particle
dark matter, Eq. (2). The Jacobian factor in Eq. (2) that describes the
particle flow is now replaced with the squared amplitude of the wave function.
For the Maxwell-Boltzmann distribution Eq. (3), the density contrast for the
wave dark matter is controlled by four parameters

$\displaystyle v_{\rm dm}/\sigma,\qquad
r/\bar{r},\qquad\hat{x}\cdot\hat{v}_{\rm dm},\qquad 2\pi\bar{r}/\lambda_{\rm
dB},$ (21)

where $\lambda_{\rm dB}=2\pi/mv$ is the de Broglie wavelength of dark matter,

$\lambda_{\rm dB}=\frac{2\pi}{mv}=10\,{\rm AU}\left(\frac{10^{-15}\,{\rm
eV}}{m}\right)\left(\frac{240\,{\rm km/sec}}{v}\right).$

As we will see in the next section, the de Broglie wavelength is a key
quantity that controls the gravitational response of the wave dark matter.

We also compute the spectrum of the field fluctuations in the frequency space.
The power and the power spectrum of the field $\hat{\phi}$ are defined as

$P_{\phi}=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}dt\;\hat{\phi}^{2}(t,{\bf
x})=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}S_{\phi}(\omega),$ (22)

where the power spectral density is

$\displaystyle
S_{\phi}(\omega)=\lim_{T\to\infty}\frac{1}{T}|\hat{\phi}(\omega)|^{2}.$ (23)

Here the Fourier component of the field is obtained via
$\hat{\phi}(\omega)=\int_{-\infty}^{+\infty}dt\,e^{i\omega t}\hat{\phi}(t)$.
Since the field is a random field, the power and power spectrum are also
random fields. The ensemble average of the spectrum is

$\displaystyle\langle
S_{\phi}(\omega)\rangle=\frac{2\pi\phi_{0}^{2}}{mv}\bar{f}(v)\bigg{|}_{v=\sqrt{2\left(\omega/m-1\right)}}$
(24)

where $\bar{f}(v)$ is given as

$\displaystyle\bar{f}(v)=\bigg{[}v^{2}\int d\Omega_{\hat{v}}\,f({\bf
v})|\psi_{\bf v}({\bf x})|^{2}\bigg{]}$ (25)

This $\bar{f}(v)$ might be interpreted as the speed distribution of wave dark
matter including the gravitational focusing effect.

In above, we focus on the variance of the field itself. Some of axion and
axion-like DM searches are based on ${\cal
L}\supset\partial_{\mu}\phi\,\bar{\psi}\gamma^{\mu}\gamma_{5}\psi$.333For
instance, CASPEr [77, 78], atomic magnetometers [79, 80], and an oscillating
neutron electric dipole moment experiment [81].  In such cases, the signal
power as well as the power spectrum would depend on the variance of the
gradient of the field. The variance of the gradient can be obtained in the
same way,

$\left\langle(\hat{n}\cdot\nabla\phi)^{2}\right\rangle=\phi_{0}^{2}\int
d^{3}v\;f({\bf v})|\hat{n}\cdot\bf{\nabla}\psi_{\bf v}({\bf x})|^{2}\ .$ (26)

where a sensitivity axis $\hat{n}$ is a unit vector determined by the specific
experimental setup. The spectrum can be obtained as

$\displaystyle\left\langle
S_{\nabla_{\hat{n}}\phi}(\omega)\right\rangle=\frac{2\pi\phi_{0}^{2}}{mv}\bigg{[}v^{2}\int
d\Omega\,f({\bf v})|\hat{n}\cdot\nabla\psi_{\bf v}({\bf x})|^{2}\bigg{]}.$
(27)

Same as the power spectrum of the field, the above expression should be
evaluated at $v=\sqrt{2(\omega/m-1)}$.

While we have mainly focused on computing variances at the same location, it
is straightforward to compute correlations of fields at different locations
with the same formalism. Such correlations contain additional information on
wave DM phase at different locations, which can be used to extract further
information, such as directionality of DM velocity distribution, and angular
distribution of potential DM substructures, from a network of wave DM
detectors [82, 83].

Figure 2: We show the density contrast for the wave dark matter in the
monochromatic limit, i.e. with vanishing velocity dispersion. We choose ${\bf
v}=(-240,\,0,\,0)$ km/sec and $m=10^{-15},\,10^{-14}\,{\rm eV}$. The source
mass is chosen to be the solar mass.

###  III.3 Monochromatic limit

To examine the gravitational response of the wave dark matter, we consider the
monochromatic limit. In this limit, the density contrast Eq. (20) becomes

$\displaystyle 1+\delta({\bf x})$ $\displaystyle=$ $\displaystyle|\psi_{{\bf
v}}({\bf x})^{2}|$ (28) $\displaystyle=$ $\displaystyle|\psi_{{\bf
v}}(0)|^{2}|M[i\beta,1,imvr(1-\hat{v}\cdot\hat{x})]|^{2}.$

Here $\psi_{{\bf v}}(0)$ is the wave function at the origin. The squared
amplitude at the origin is

$\displaystyle|\psi_{{\bf v}}(0)|^{2}=\frac{2\pi\beta}{1-e^{-2\pi\beta}}.$
(29)

The above is identical to the Sommerfeld enhancement for dark matter
annihilation when the interaction is mediated by the Coulomb interaction [84].

Figure 3: The density contrast at $\hat{x}\cdot\hat{v}=0$. The constant
density contrast is observed for small radii satisfying $m\tilde{v}r<1$. The
black line corresponds to $m\tilde{v}r=1$.  Figure 4: The ratio between the
density contrast averaged over the solid angle for wave and particle dark
matter. For large radii, $m\tilde{v}r\gtrsim 1$, the averaged density contrast
for wave dark matter approaches that of the particle dark matter. The solid
black line represents $m\tilde{v}r=1$.

In Figure 2, we show the density contrast of wave dark matter in the
monochromatic limit. We choose ${\bf v}=(-240,\,0,\,0)$ km/s, $m=10^{-15}$,
$10^{-14}$ eV, and $M=M_{\odot}$. We observe a similar density wake field
along the direction of dark matter velocity.

Some of wave properties in the density wake are revealed upon a closer
examination of Figures 2–4. Contrary to the particle dark matter, the density
contrast saturates to a constant value, $\delta\simeq|\psi_{{\bf
v}}(0)|^{2}-1$, for sufficiently small $r$. In addition, there is an
oscillating pattern of density contrast over a certain specific spatial scale.
Both of wave characteristics are controlled by de Broglie wavelength defined
as

$\displaystyle\tilde{\lambda}_{\rm dB}(r)=\frac{2\pi}{m\tilde{v}(r)}$ (30)

where $\tilde{v}$ is the speed at $r$, given as

$\displaystyle\tilde{v}^{2}(r)=v^{2}+v_{e}^{2}(r).$ (31)

Note that we have used the speed at a radial distance $r$, rather than the
velocity asymptotically far away from the source $M$.

To illustrate how this de Broglie wavelength controls the monochromatic
result, we show in Figure 3 the density contrast normalized with the central
value, $\delta({\bf x})/\delta(0)$, as a function of two parameters,
$r/\bar{r}$ and $2\pi\bar{r}/\lambda_{\rm dB}$. Here $\lambda_{\rm
dB}=2\pi/mv$ where $v$ is the dark matter velocity far away from the source.
The black line in the figure is $\tilde{\lambda}_{\rm dB}=2\pi r$, i.e.
$m\tilde{v}(r)r=1$. We observe that the density contrast saturated to
$\delta(0)$ for sufficiently small r, $m\tilde{v}(r)r\lesssim 1$. For
$m\tilde{v}(r)r\gtrsim 1$, the wave density contrast oscillates rapidly but on
average, it behaves similar to that of particle dark matter.

In Figure 4, we plot the ratio of density contrast averaged over the solid
angle for the wave and particle dark matter. For particle dark matter, we
choose $1+\delta_{\rm avg,\,p}=\sqrt{1+v_{e}^{2}/v^{2}}$. We observe that for
$m\tilde{v}(r)r\gtrsim 1$, the averaged density contrast for the wave dark
matter becomes increasingly similar to the particle result, while for
$m\tilde{v}(r)r\lesssim 1$, the averaged density contrast is suppressed as
$1+\delta_{\rm avg}\approx\pi(2\pi r/\tilde{\lambda}_{\rm
dB})(1+v_{e}^{2}/v^{2})^{1/2}$.

In both figures, we observe that a key parameter controlling the wave dark
matter density contrast is the relative size of the de Broglie wavelength
$\tilde{\lambda}_{\rm dB}$ and $r$. For $m\tilde{v}r\gtrsim 1$, the response
of wave dark matter to the gravitational potential is similar to that of
particle dark matter, while for $m\tilde{v}r\lesssim 1$, wave characteristics
appear. The limit $m\tilde{v}r\gg 1$ is the limit at which the wave dark
matter approaches to particle dark matter. This will be discussed in more
detail in Section V.2 and in Appendix B.2.

###  III.4 Comparison

The monochromatic results clearly show the role of de Broglie wavelength. We
now consider the density contrast of wave dark matter with velocities
distributed according to the Maxwell-Boltzmann distribution to compare the
wave dark matter with particle dark matter on an equal footing.

In Figure 5, we present two quantities, the density contrast along the
downstream (top) and the density contrast averaged over the solid angle
(bottom). For both plots, we choose $v_{\rm dm}/\sigma=1$. As in the
monochromatic case, particle-to-wave transition takes place when
$m\tilde{v}r\simeq 1$. Since we consider a medium with nonvanishing velocity
dispersion, we should be careful interpreting $\tilde{v}^{2}=v^{2}+v_{e}^{2}$
since the final density contrast is the result of the superposition of
different velocity components. For the following discussion, we redefine
$\tilde{v}^{2}=v_{\rm dm}^{2}+\sigma^{2}+v_{e}^{2}$ for the medium with
nonzero velocity dispersion. In the figure, the black dashed line represents
the result of particle dark matter, while colored lines represent the result
of wave dark matter with $m\sigma\bar{r}=0.5$ (cyan), $1$ (orange), $2$
(purple). For $m\tilde{v}r\lesssim 1$, the density contrast saturates to

$\displaystyle 1+\delta(0)=\int d^{3}v\,f({\bf v})|\psi_{{\bf
v}}(0)|^{2}\approx\frac{2\pi\langle\beta\rangle}{1-e^{-2\pi\langle\beta\rangle}}$
(32)

with $\langle\beta\rangle=(GMm/v_{\rm dm})\,{\rm erf}(v_{\rm dm}/2\sigma)$.
For $m\tilde{v}r\gtrsim 1$, the wave result converges to the particle dark
matter result.

Figure 5: (Top) We show the density contrast along the downstream. Each color
represents different values of $m\sigma\bar{r}$. The circle represents the
point $m\tilde{v}(r)r=1$. The black dashed line is the result of particle dark
matter obtained in Section II. We observe that the density contrast approaches
the particle dark matter result for $m\tilde{v}r\gtrsim 1$, while it saturates
to a constant value for $m\tilde{v}r\lesssim 1$. (Bottom) We show the averaged
density contrast. For both figures, we choose $v_{\rm dm}/\sigma=1$.

##  IV Dark matter in solar system

We apply our previous discussion to dark matter structures in the solar
system. We work in the Galactic coordinate system where the origin of
coordinates is at the location of the Sun. In the rectangular Cartesian
coordinate system, each component of $(X,Y,Z)$ points towards the Galactic
center, the direction of Galactic rotation, and the Galactic north pole. In
this coordinate system, we compute the density contrast at the position of
Earth in different times in a year for dark matter substructures with
different kinematic properties.

For the velocity distribution, we use the Maxwell-Boltzmann distribution, Eq
(3), for most of dark matter substructures. The distribution is therefore
characterized by the mean velocity $\langle{\bf v}\rangle={\bf v}_{\rm dm}$
and the one-dimensional velocity dispersion $\sigma$.

###  IV.1 Halo dark matter

We begin with the relaxed halo dark matter component. In the Galactic
coordinate, the halo dark matter can be represented by the standard halo model
with a mean velocity in the Galactic coordinate [85, 86],

$\displaystyle{\bf v}_{\rm dm}=-{\bf v}_{\odot}=-(11,\,241,\,7)\,{\rm km/sec}$
(33)

The velocity dispersion is

$\displaystyle\sigma=v_{c}(R_{\odot})/\sqrt{2}=162\,{\rm km/sec}$ (34)

where $v_{c}(R_{\odot})\simeq 229\,{\rm km/sec}$ is the circular velocity of
Milky Way at the position of the Sun [86].

Figure 6: Density contrast for halo dark matter at the position of Earth
during a year 2021. The black dashed line is for the particle dark matter,
while each colored line represents the wave dark matter result for
$\tilde{\lambda}_{\rm dB}=20,\,3,\,0.5$ AU. Corresponding wave DM mass with
respect to $m_{15}=10^{-15}\,{\rm eV}$ is also shown. The density contrast
flatten for $m\tilde{v}r\lesssim 1$, while it approaches the particle dark
matte result for $m\tilde{v}r\gtrsim 1$ similarly to the halo dark matter
example. Around $m\tilde{v}r\sim 1$, the wave density contrast becomes $\sim
25~{}\%$ larger than that of particle dark matter.

In Figure 6, we show the density contrast of wave dark matter at the position
of the Earth throughout the year of 2021 for different values of
$m\tilde{v}r$. Since $v_{e}\simeq 42\,{\rm km/sec}\ll\sigma,v_{\rm dm}$, we
expect that the overall amplitude of density contrast is $\delta\simeq
v_{e}^{2}/2\sigma^{2}\simeq 0.03$ at least for the particle dark matter, which
is confirmed by the black dashed line in the figure. The maximum of the
density contrast takes place around March 1st since
$\mu=\hat{x}\cdot\hat{v}_{\rm dm}$ takes the largest value over a year. This
pattern can be seen also from the wave dark matter density contrast with
varying degree of annual modulation strength. While the maximum and minimum
arise at the same time, we observe that the density contrast flattens for
$m\tilde{v}r\lesssim 1$ as the de Broglie wavelength becomes larger than an
astronomical unit. For $m\tilde{v}r\gg 1$, the result can be approximated as
that of particle dark matter density contrast. We also observe that the wave
density contrast is enhanced by $\sim 25\%$ compared to the particle density
contrast at $m\tilde{v}r\sim 1$.

Figure 7: (Top) The difference of speed distribution with and without the
gravitational focusing, computed at March 1st, 2021. (Bottom) Same as top
panel but computed at September 1st, 2021. Each colored line represents
$\tilde{\lambda}_{\rm dB}=20,\,3,\,0.5$ AU. The black dashed line is the same
result for the particle dark matter, and the black dotted line is the rescaled
speed distribution without gravitational focusing.

In Figure 7, we also show the modification of the speed distribution for the
wave dark matter. More specifically, we plot

$\displaystyle\Delta\bar{f}(v)=v^{2}\int d\Omega_{\hat{v}}f({\bf
v})(|\psi_{{\bf v}}({\bf x})|^{2}-1),$ (35)

which is the difference between the speed distribution with and without the
gravitational focusing. We choose March 1st and September 1st for the plots.
As it is clear from the figures, the effect of gravitational focusing is
prominent in the low velocity tail. Each colored line represent the result of
wave dark matter for different values of $m\tilde{v}r$. The black dashed line
is the modification of speed distribution in the case of particle dark matter.
The dotted line is rescaled speed distribution without gravitational focusing.

The variance of the projection of the field gradient along the sensitivity
axis of a detector for the halo DM component is shown in Appendix C.

###  IV.2 Gaia Sausage

We consider a dark matter substructure with anisotropic velocity distribution.
Such a dark matter substructure is motivated by the observation of a stellar
population with an anisotropic velocity distribution, referred to as Gaia-
Sausage or Gaia-Enceladus [25, 26]. It is inferred that this stellar component
originates from a relatively recent merger with a luminous satellite galaxy of
mass $M>10^{10}M_{\odot}$. It is naturally expected that not only stars but
also dark matter is accreted on to the Milky Way from the same event. Studies
have shown that the kinematic properties of accreted stars and dark matter are
similar [27, 28] and that it could consitute ${\cal O}(10\%)$ of local dark
matter [29].

We inspect the response of sausage-component dark matter to the solar
gravitational potential. Since the velocity structure is anisotropic, we use a
general three-dimensional normal distribution for the DM velocity
distribution,

$\displaystyle f({\bf v})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{3/2}\sqrt{\det\Sigma}}$ (36)
$\displaystyle\times\exp\left[-\frac{1}{2}({\bf v}-{\bf v}_{\rm
dm})\cdot\Sigma^{-1}\cdot({\bf v}-{\bf v}_{\rm dm})\right]$

where $\Sigma$ is the covariance matrix. We parameterize sausage dark matter
with ${\bf v}_{\rm dm}=-{\bf v}_{\odot}$ and

$\displaystyle\sigma_{r}=256\,{\rm
km/sec},\quad\sigma_{\theta}=\sigma_{\phi}=81\,{\rm km/sec},$ (37)

where $\Sigma={\rm
diag}(\sigma_{r}^{2},\sigma_{\theta}^{2},\sigma_{\phi}^{2})$ [31]. This
covariance matrix reflects the high velocity anisotropy of this substructure.

Figure 8:  Same as Figure 6 but for sausage dark matter component. Around
$\tilde{\lambda}_{\rm dB}\simeq 2$ AU, corresponding to $m\tilde{v}r\simeq 3$,
the wave dark matter contrast can be almost a factor two larger than the
particle density contrast. We also observe a slight underdensity around fall.

In Figure 8, we present the density contrast for the sausage dark matter
component. We observe a similar pattern as in the example of halo dark matter.
Overall amplitude is similar to the halo dark matter component, despite that
the velocity distribution is highly anisotropic. The maximum and the minimum
of the density contrast also take place at the same period of the time. A
slight underdensity is also observed around fall.

###  IV.3 Dark disk

As a next example, we consider dark matter in the form of disk. This dark
matter substructure is motivated by the thick stellar disk. To explain the
vertical structure of stellar disk, several scenarios have been suggested. For
instance, thick stellar disk is composed of stars that are accreted from
satellite galaxies [87] or stars in pre-existing thin disk are heated through
merger events [88, 89]. It has been suggested that such merger events could
naturally lead to the formation of thick dark disk as the accreted dark matter
is dragged to the disk plane through dynamical friction [32]. The resulting
dark matter disk is co-rotating with Galactic disk with slightly smaller
circular velocity by $\Delta v=50$ km/sec. While it is still not clear whether
the aforementioned merger scenario can solely explain the evolution of stellar
disk (see a review [90] and references therein), it is still worthwhile to
investigate the response of such dark matter component to the solar
gravitational potential.

Figure 9: Same as Figures 6 and 8, but for the dark disk. Due to the small
mean velocity and dispersion, the overall amplitude of density contrast is
larger than the previous DM structures.

For the investigation of dark disk, we parameterize the mea velocity and
velocity dispersion as

${\bf v}_{\rm dm}=(0,-50,0)~{}{\rm
km/sec}\quad\textrm{and}\quad\sigma=50\,{\rm km/sec}.$

While we refer this structure to dark disk, it could also represent any cold
dark matter substructure with a small mean velocity.444The velocity dispersion
of thick stellar disk is
$(\sigma_{r},\sigma_{\theta},\sigma_{\phi})\simeq(63,39,39)$ km/sec [91], but
we assume an isotropic distribution for dark matter for simplicity. In
addition, even in the merger scenario, the dark matter kinematic structure
depends on the merger details. Numerical studies have shown that the DM
kinematic structure follows that of accreted stars, and the accreted stars are
hotter than heated thin disk stars, indicating that the velocity dispersion of
dark matter is likely to be higher than that of stellar thick disk [32, 92].
Nonetheless, we choose $\sigma=50$ km/sec as this benchmark more clearly shows
the characteristic behavior of gravitational focusing.

In Figure 9, we show the density contrast for the dark disk. The pattern of
the density contrast is similar to previous cases, but with a larger amplitude
due to the small velocity and dispersion. In this particular example, the disk
dark matter can be focused and its density near the orbital trajectory of
Earth could be $\sim 30\%$ larger than its density far away from the Sun.

The variance of the field gradient in the case of a dark disk component is
shown in Appendix C.

###  IV.4 Stream

As a last example, we consider dark matter streams, whose kinematic structure
is similar to that of stellar stream. Stellar streams are stellar
substructures and they are coherent both spatially and kinematically. If
stellar streams originate from dwarf galaxies, it is likely that there are
associated dark matter streams with similar kinematic properties [28]. In the
rest frame of Milky Way dark matter halo, streams typically have large
streaming velocity with a small velocity dispersion. Various streams and
substructures have been observed in the inner halo in the past decades. See a
review [93] for details on stellar streams and substructures in Milky Way.

To investigate the gravitational response of fast-moving cold objects like
streams, we consider

$v_{\rm dm}=400~{}{\rm km/sec},\quad\sigma=30\,{\rm km/sec}.$

Note that the mean speed in the rest frame of the Sun might be smaller than
the above value for prograde streams. In such cases, the density contrast as
well as the spectrum would resemble that of dark disk. For the direction of
the stream, we consider two different inclination angles with respect to the
ecliptic plane, $\theta_{\rm dm}=0,\,\pi/6$.

Figure 10: (Top) The density contrast for the stream for the vanishing
inclination angle. Due to the hierarchy between the mean velocity and the
dispersion, the density contrast is highly localized near March 1st. The wave
dark matter contrast for $\tilde{\lambda}_{\rm dB}=0.5$ AU, corresponding to
$m\tilde{v}r=10$, is still quite different from the particle dark matter
result. See the main text for the explanation. (Bottom) The density contrast
for the stream for the inclination angle $\theta_{\rm dm}=\pi/6$. A small
velocity dispersion for stream does not completely erase the oscillations of
the wave function, and therefore, much faster oscillation of density contrast
is obtained in this example.

In Figure 10, we show that the density contrast for dark matter streams for
two different inclination angles. For the case of vanishing inclination angle
$\theta_{\rm dm}=0$, we observe interesting differences between the particle
dark matter and wave dark matter. Even for $m\tilde{v}r=10$, the wave contrast
is much broader than the particle contrast in the time-axis. We have argued
that the wave dark matter contrast approaches to the particle contrast for
$m\tilde{v}r\gg 1$, based on quantities like $\delta_{\rm avg}$, but note that
$\delta_{\rm avg}$ is the averaged over the solid angle. At the same time, we
have seen both for the particle and wave dark matter that there exists a
certain angular scale $\Delta\mu$ over which the density contrast takes
nonvanishing value. For the particle dark matter, this angular scale is
typically $\Delta\mu\sim(\sigma/v_{\rm dm})^{2}$ because of the non-vanishing
velocity dispersion, while for the wave dark matter, due to the de Broglie
wavelength, it is $\Delta\mu\sim\max[1/m\tilde{v}r,(\sigma/v_{\rm dm})^{2}]$.
The stream that we consider has a large hierarchy between $v_{\rm dm}$ and
$\sigma$, and hence, the particle density contrast is nonvanishing only for a
relatively small period of time in a year, i.e. $\Delta t\sim(\sigma/v_{\rm
dm})^{2}{\rm yr}$, while for the wave dark matter, for $m\tilde{v}r=10$, the
angular spread $\Delta\mu$ is still controlled by $1/m\tilde{v}r$ rather than
$(\sigma/v_{\rm dm})^{2}$. Therefore, the wave contrast in Figure 10 is less
localized in time (or in space) compared to the particle contrast.

This same feature can also be observed in the example of $\theta_{\rm
dm}=\pi/6$. In this case, the density contrast for the particle dark matter
takes nearly vanishing value, while for the wave dark matter, it takes still
non-vanishing value at the level of a few percent. Another interesting feature
is that the pattern of density contrast over a year yields much faster
oscillations, which is clearly distinct from other substructures. This can be
attributed to the small velocity dispersion; the wave behaves more a like
monochromatic wave and the characteristic oscillation patterns survive even
after the waves with different wavenumber are mixed together.

##  V Discussion

###  V.1 Coordinate transformation

We have worked in the Galactic coordinate system. This coordinate is useful
for the investigation of the effects of gravitational focusing since the
gravitational focusing effects can be easily isolated. Since terrestrial
experiments are bound to the Earth, it is worthwhile to mention how the
observables such as the density contrast and spectrum change in the detector
rest frame.

The density contrast does not change at all. Let us denote $x$ for the
Galactic coordinate and $x^{\prime}$ for the proper coordinate of detector.
The density contrast is defined from the two point function of the scalar
field $\langle\hat{\phi}^{2}(x)\rangle$, and since this quantity is a scalar
under the general coordinate transformation, it should be the same in both
frames, i.e. $\langle\hat{\phi}^{2}(x)\rangle=\langle\hat{\phi}^{\prime
2}(x^{\prime})\rangle$. Note that the density contrast is related to the total
power of the signal. This indicates that the density contrast obtained in the
previous section is directly related to the total power of the signal measured
in detectors on Earth.

The spectrum needs more careful discussion. The frequency $\omega$ is
associated with the coordinate time $t$ of the Galactic frame. The Fourier
transform in the Galactic frame is defined as $\hat{\phi}(\omega)=\int
dt\,e^{i\omega t}\hat{\phi}(x)$, whereas it is defined as
$\hat{\phi}^{\prime}(\omega^{\prime})=\int
dt^{\prime}\,e^{i\omega^{\prime}t^{\prime}}\hat{\phi}(x)$ in the detector
frame with a proper time $t^{\prime}$. If the detector is an inertial observer
moving in a constant velocity with respect to the Galactic frame, the proper
coordinate can be related to the Galactic coordinate as $x^{\prime}=\Lambda x$
with an appropriate Lorentz transformation $\Lambda$. In this case, it is
straightforward to show that the same expression for the spectrum, Eq. (24),
can be used in the detector frame upon a replacement of the velocity in the
distribution with the velocity observed in the detector frame,
$k\to\Lambda^{-1}k$.

The orbital motion of the Earth introduces complications to the above
discussion since a detector on Earth is not an inertial observer with a
constant velocity. The above procedure of replacing $k\to\Lambda^{-1}k$ to
obtain the spectrum in the detector frame is only valid if the time series is
extended over a time scale sufficiently shorter than the orbital time scale of
the Earth. If one attempts to Fourier-transform the signal over a time scale
longer than orbital time scale or the coherence time scale of wave dark
matter, this procedure would introduce side-bands in the frequency space due
to the frequency modulation induced by the orbital motion of Earth and the
oscillations of the wave function.555For experiments sensitive to the
gradient, the rotation of Earth introduces an amplitude modulation of the
signal, leading to side-bands in the spectrum. This also introduces additional
complications in the frequency space analysis [50].

###  V.2 Semiclassical limit

In the previous section, we have observed that the wave result approaches the
classical particle dark matter result when

$\displaystyle m\tilde{v}(r)r\gg 1.$ (38)

This is the limit at which the size of the system becomes much greater than
the de Broglie wavelength of the wave dark matter. We provide a brief
explanation why both results converge in this limit.

The discussion of particle dark matter is entirely based on the phase space
distribution modified by the gravitational potential of the source. To discuss
the limit at which particle and wave dark matter results converge, it is
convenient to find a quantity analogous to the classical phase space
distribution. The quantum mechanical analog of the classical phase space
distribution is given by the following quasi-probability distribution,

$\displaystyle f_{W}({\bf x},{\bf p})$ $\displaystyle=$ $\displaystyle\int
d^{3}y\,e^{i{\bf p}\cdot{\bf y}/\hbar}$ (39) $\displaystyle\times$
$\displaystyle\int d^{3}v\,f({\bf v})\psi^{*}_{\bf v}({\bf x}+{\bf
y}/2)\psi_{\bf v}({\bf x}-{\bf y}/2).$

which was first introduced by Wigner [94]. We reintroduce the Planck constant
$\hbar$. When it is integrated over the phase space $d^{3}p$, it reproduces
the variance of the field, which is directly related to the density contrast
as

$\displaystyle\int\frac{d^{3}p}{(2\pi)^{3}}f_{W}({\bf x},{\bf p})=\int
d^{3}v\,f({\bf v})|\psi_{{\bf k}}({\bf x})|^{2}=\langle\phi^{2}\rangle.$ (40)

Moreover, the quasi-probability distribution $f_{W}$ evolves according to

$\displaystyle\frac{\partial f_{W}}{\partial t}$ $\displaystyle+$
$\displaystyle{\bf v}\cdot\frac{\partial f_{W}}{\partial{\bf x}}$
$\displaystyle+$ $\displaystyle\frac{i}{\hbar}\left[V\big{(}{\bf
x}+\frac{i\hbar}{2}\nabla_{{\bf p}}\big{)}-V\big{(}{\bf
x}-\frac{i\hbar}{2}\nabla_{{\bf p}}\big{)}\right]f_{W}=0.$

This can be derived from the definition of $f_{W}$, Eq. (39), and the
Schrödinger equation with potential $V$. We see that the leading term in the
semiclassical limit ($\hbar\to 0$) is identical to the Boltzmann equation.
From above, we can conclude that the wave and particle dark matter result
should converge in the limit $\hbar\to 0$. This semiclassical limit is
identical to Eq. (38), which can be easily shown from the Schrödinger
equation. See Appendix B.2 for more details.

##  VI Conclusion

We have discussed the gravitational response of wave dark matter. After a
brief review on the gravitational response of particle dark matter, we have
discussed and provided a simple formalism to investigate wave dark matter in
the presence of the gravitational potential due to astrophysical object of
mass $M$. We specifically discuss a spinless bosonic field, minimally coupled
to the gravity. As in the canonical formulation of quantum field theories, we
have expanded the dark matter field in terms of the creation and annihilation
operators with a mode function. We solve the Schrödinger equation to find the
response of each wave mode and then specify the density operator which
completely sets the statistical properties of the field. This procedure allows
us to consider both the gravitational response and the stochastic nature of
the wave dark matter.

A comparison of particle and wave dark matter reveals interesting similarity
and at the same time differences between them. It is already clear from the
monochromatic result that the wave nature becomes manifest only for
$m\tilde{v}r\lesssim 1$, while for $m\tilde{v}r\gg 1$, it practically
converges to the result of particle dark matter. The characteristics of wave
dark matter is the appearance of an additional scale, the de Broglie
wavelength, $\tilde{\lambda}_{\rm dB}=2\pi/m\tilde{v}$. The spatial structure
of the wave density contrast is smoothened on the scale of its de Broglie
wavelength, which is the feature commonly observed in the studies of wave dark
matter.

We have applied the formalism to the dark matter substructures as well as halo
dark matter in the solar system. Different substructures are considered to
examine the gravitational response of wave dark matter as a function of the
mean velocity and the velocity dispersion. Substructures considered here are
motivated either by recent astrometric data or by a theoretical possibility to
explain the stellar substructures in the Milky Way. The wave characteristics
are clearly seen in all of the provided examples: for $m\tilde{v}r\gg 1$, the
wave response is essentially the same as the particle response, while for
$m\tilde{v}r\ll 1$, the wave result becomes independent on the location of the
Earth. Interestingly, for $m\tilde{v}r\sim 1$, we observe that the wave
density contrast could be at most a factor two larger than that of particle
dark matter. Moreover, for streams with small velocity dispersion, the angular
distribution of the density contrast could still be different from that of
particle density contrast. Although the gravitational focusing effects remain
subleading in many examples we studied, they could provide important clues for
the local dark matter structures upon the detection of dark matter.

###### Acknowledgements.

This work is supported by the Deutsche Forschungsgemeinschaft under Germany’s
Excellence Strategy - EXC 2121 Quantum Universe - 390833306. We thank the
hospitality of INFN Galileo Galilei Institute for Theoretical Physics during
the workshop ‘New Physics from The Sky’.

##  Appendix A Particle gravitational focusing

We discuss details of the gravitational response of particle dark matter. In
the main text, we have provided two approximations for (i) the density
contrast along the downstream Eq. (5)

$1+\delta_{\rm
ds}=1+\delta|_{\mu=1}\approx\Big{(}1+\frac{v_{e}^{2}}{\sigma^{2}}\Big{)}^{\frac{1}{2}}$

and (ii) the density contrast averaged over the solid angle Eq. (6)

$1+\delta_{\rm
avg}=\int\frac{d\Omega_{\hat{x}}}{4\pi}(1+\delta)\approx\Big{(}1+\frac{v_{e}^{2}}{\sigma^{2}+v_{\rm
dm}^{2}}\Big{)}^{\frac{1}{2}}.$

In Figure 11, we compare these approximations and numerical results. In the
top panel, we plot $\delta_{\rm ds}$ for $v_{\rm dm}/\sigma\in[0.1,10]$. The
result shows that $\delta_{\rm ds}$ does not crucially depend on the value of
$v_{\rm dm}/\sigma$. The dashed line is the approximation given by Eq. (5),
and the numerical result matches well to the approximation within a factor of
two for a wide range of $v_{\rm dm}/\sigma$. In the bottom panel, we plot
$\delta_{\rm avg}$ for the same range of $v_{\rm dm}/\sigma$. The dashed line
is the approximation given by Eq. (6). The result $\delta_{\rm avg}$ is only
sensitive to $r(\sigma^{2}+v_{\rm dm}^{2})/GM$.

Figure 11: (Top) We show the density contrast along the downstream
$\mu=\hat{x}\cdot\hat{v}_{\rm dm}=1$ for different values of $v_{\rm
dm}/\sigma\in[0.1,10]$ (red band). The lower and upper boundary correspond to
$v_{\rm dm}/\sigma=0.1$, $10$, respectively. The dashed line corresponds to
Eq. (5), $1+\delta_{\rm ds}=\sqrt{1+v_{e}^{2}/\sigma^{2}}$. The numerical
result is consistent to the approximation Eq. (5) within a factor two for a
wide range of $v_{\rm dm}/\sigma$. (Bottom) We show the density contrast
averaged over the solid angle for $v_{\rm dm}/\sigma\in[0.1,10]$ (red band).
The lower and upper boundary correspond to $v_{\rm dm}/\sigma=0.1,\,10$,
respectively. The dashed line corresponds to the approximation Eq. (6). The
numerical result is consistent to the approximation Eq. (6) within a factor
two for a wide range of $v_{\rm dm}/\sigma$.

In the limit $v_{e}\ll\sigma$, an analytic expression for the density contrast
can be obtained [95],

$\delta({\bf x})\approx\frac{v_{e}^{2}}{2\sigma^{2}}\exp\left[-\frac{v_{\rm
dm}^{2}}{2\sigma^{2}}\sin^{2}\theta\right]{\rm erfc}\left(-\frac{v_{\rm
dm}}{\sqrt{2}\sigma}\cos\theta\right)$

where $\cos\theta=\hat{x}\cdot\hat{v}_{\rm dm}$. From this, one may obtain the
approximations Eqs. (5)–(6) for $v_{e}/\sigma<1$. This expression can be used
for an investigation of dark matter substructures with $\sigma>v_{e}$.

##  Appendix B Wave gravitational focusing

###  B.1 Average density contrast in monochromatic limit

In Figure 12, we present the density contrast averaged over the solid angle
$\hat{x}$ in the monochromatic limit for different values of $r/\bar{r}$. The
averaged density contrast is normalized with $\delta_{\rm
avg,p}(r)=\sqrt{1+v_{e}^{2}(r)/v^{2}}+1$. We observe that for $m\tilde{v}r$
the wave dark matter result is suppressed relative to the particle dark matter
result as $\delta_{\rm avg,w}/\delta_{\rm avg,p}\propto\pi m\tilde{v}r$. The
behavior of $\delta_{\rm avg,w}/\delta_{\rm avg,p}$ does not sensitively
depend on the values of $r/\bar{r}$.

Figure 12: We present the density contrast averaged over the solid angle in
the monochromatic case. Different solid lines correspond to different values
of $r/\bar{r}$. All of them approach to the classical particle result for
$m\tilde{v}r\gg 1$. For $m\tilde{v}r\lesssim 1$, the averaged density contrast
is suppressed.

###  B.2 Semiclassical limit

In this appendix, we discuss the semiclassical limit in the Schrödinger
equation. We reintroduce the Planck constant $\hbar$. We begin with the
Schrödinger equation, $i\hbar\dot{\psi}=(-\hbar^{2}\nabla^{2}/2m+V)\psi$. With
$\psi(t,{\bf x})\to e^{-iEt}\psi({\bf x})$, the Schrödinger equation becomes

$\displaystyle 0=\left[\nabla^{2}+\frac{2m(E-V)}{\hbar^{2}}\right]\psi.$ (41)

The energy eigenvalue is $E=k^{2}/2m$ where ${\bf k}$ is the wave number
asymptotically far way from the source $M$. The wave function is characterized
by ${\bf k}$. We perform Legendre expansion of the field,

$\psi_{{\bf k}}({\bf
x})=\frac{1}{r}\sum_{\ell=0}^{\infty}\chi_{k\ell}(r)P_{\ell}(\cos\theta),$

where $\cos\theta=\hat{x}\cdot\hat{v}$. This procedure leads to the
Schrödinger equation for the radial wave function $\chi_{k\ell}$,

$\displaystyle\chi_{k\ell}^{\prime\prime}+\frac{2m(E-V_{\ell})}{\hbar^{2}}\chi_{k\ell}=0,$
(42)

where the prime denotes the derivative with respect to $r$ and the potential
$V_{\ell}$ is defined as

$\displaystyle
V_{\ell}(r)=V(r)+\frac{\hbar^{2}}{2m}\frac{\ell(\ell+1)}{r^{2}}=V(r)+\frac{L^{2}}{2mr^{2}}.$
(43)

We have defined $L^{2}=\hbar^{2}\ell(\ell+1)$.

We seek for the solution of the following form,

$\displaystyle\chi_{k\ell}(r)=e^{i\sigma/\hbar}=\exp\left[\frac{i}{\hbar}\sum\hbar^{n}\sigma^{(n)}_{\ell}\right].$
(44)

The Schrödinger becomes $\sigma^{\prime
2}-i\hbar\sigma^{\prime\prime}=2m(E-V_{\ell})$. By solving the Schrödinger
equation order by order in $\hbar$, we find

$\displaystyle\sigma_{\ell}^{(0)}(r)$ $\displaystyle=$
$\displaystyle\int^{r}dr^{\prime}\,\sqrt{2m[E-V_{\ell}]}\equiv\int^{r}dr^{\prime}\,p_{\ell}(r^{\prime}),$
(45) $\displaystyle\sigma_{\ell}^{(1)}(r)$ $\displaystyle=$
$\displaystyle\frac{i}{2}\ln p_{\ell}(r).$ (46)

We have defined $p_{\ell}(r)\equiv\sqrt{2m[E-V_{\ell}(r)]}$, which is
identical to the radial momentum in $1/r$-potential in classical mechanics.

This semi-classical or WKB approximation is valid when $|\hbar
p_{\ell}^{\prime}/p_{\ell}^{2}|\ll 1$. This can be written as

$\displaystyle\frac{d}{dr}\frac{\hbar}{p_{\ell}}\ll 1.$ (47)

Naively speaking, this condition is approximately
$\hbar/p_{\ell}r\sim\lambda_{\rm dB}/r\ll 1$, indicating that the semi-
classical limit is valid only when the size of the system is much greater than
the de Broglie wavelength, $\lambda_{\rm dB}\sim 1/p$. More specifically, the
above condition can be written as

$\displaystyle\frac{\hbar m|F_{\ell}|}{p_{\ell}^{3}}\ll 1$ (48)

with

$\displaystyle
F_{\ell}=p_{\ell}^{\prime}=\frac{1}{2mr^{3}}\left[(mv_{e}r)^{2}-2L^{2}\right]$
(49)

For a given set of parameters, the dominant contribution in the series
expansion of wave function arises at $\ell\sim m\tilde{v}r$. With
$m\tilde{v}r>mv_{e}r$, it is straightforward to show that this condition Eq.
(48) can be approximated as $m\tilde{v}r\gg\hbar$.

Figure 13: Daily modulation of the variance of the field gradient projected on
the sensitivity axis of a detector located in Hamburg, Germany and pointing in
the east direction for halo DM (top) and dark disk (bottom) on March 1st,
2021. Colored lines are the results with gravitational focusing focusing for
different values of $\tilde{\lambda}_{\rm dB}$ as in Figure 6 (top) and 9
(bottom), respectively. The black dotted line denotes the asymptotic result.
The bottom panels shows the contrast of focused result with respect to
unfocused ones.

##  Appendix C Field gradient daily modulation

In Figure 13, we show the daily modulation of the field gradient calculated
with Eq. (26). As a benchmark we consider a detector located in Hamburg,
Germany $(54^{\circ}{\rm N},10^{\circ}{\rm E})$ with sensitivity axis pointing
towards the east direction on March 1st, 2021. We show the result for halo DM
and dark disk component, with same benchmark choices for $m\tilde{v}r$
exploited in Figure 6 (top) and 9 (bottom). The dotted black line shows the
result with no focusing, asymptotically far away from the solar system. In the
bottom panels we show the contrast $\delta_{\nabla}$ between the focused lines
and the asymptotic value. Note that in both cases the magnitude of the
contrast is comparable to the density contrast, signaling the effect is of
order $\mathcal{O}(v_{e}^{2}/\sigma^{2})$. Another interesting effect due to
gravitational focusing is the modification of the position of the peaks. This
effect is due to the velocity component acquired by DM from the gravitation
potential of the Sun, which changes the direction of the total DM velocity
vector, hence the time at which the peak value is reached.

## References

  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie [1983] L. F. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120, 137 (1983).
  * Peccei and Quinn [1977] R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38, 1440 (1977).
  * Weinberg [1978] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40, 279 (1978).
  * Kim [1979] J. E. Kim, Weak Interaction Singlet and Strong CP Invariance, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ [1980] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Can Confinement Ensure Natural CP Invariance of Strong Interactions?, Nucl. Phys. B 166, 493 (1980).
  * Dine  _et al._ [1981] M. Dine, W. Fischler, and M. Srednicki, A Simple Solution to the Strong CP Problem with a Harmless Axion, Phys. Lett. B 104, 199 (1981).
  * Zhitnitsky [1980] A. R. Zhitnitsky, On Possible Suppression of the Axion Hadron Interactions. (In Russian), Sov. J. Nucl. Phys. 31, 260 (1980).
  * Graham  _et al._ [2015a] P. W. Graham, D. E. Kaplan, and S. Rajendran, Cosmological Relaxation of the Electroweak Scale, Phys. Rev. Lett. 115, 221801 (2015a), arXiv:1504.07551 [hep-ph] .
  * Arvanitaki  _et al._ [2017] A. Arvanitaki, S. Dimopoulos, V. Gorbenko, J. Huang, and K. Van Tilburg, A small weak scale from a small cosmological constant, JHEP 05, 071, arXiv:1609.06320 [hep-ph] .
  * Banerjee  _et al._ [2019] A. Banerjee, H. Kim, and G. Perez, Coherent relaxion dark matter, Phys. Rev. D 100, 115026 (2019), arXiv:1810.01889 [hep-ph] .
  * Arkani-Hamed  _et al._ [2021] N. Arkani-Hamed, R. Tito D’agnolo, and H. D. Kim, Weak scale as a trigger, Phys. Rev. D 104, 095014 (2021), arXiv:2012.04652 [hep-ph] .
  * Tito D’Agnolo and Teresi [2021a] R. Tito D’Agnolo and D. Teresi, Sliding Naturalness, (2021a), arXiv:2106.04591 [hep-ph] .
  * Tito D’Agnolo and Teresi [2021b] R. Tito D’Agnolo and D. Teresi, Sliding Naturalness: Cosmological Selection of the Weak Scale, (2021b), arXiv:2109.13249 [hep-ph] .
  * Svrcek and Witten [2006] P. Svrcek and E. Witten, Axions In String Theory, JHEP 06, 051, arXiv:hep-th/0605206 .
  * Arvanitaki  _et al._ [2010] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String Axiverse, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Hui [2021] L. Hui, Wave Dark Matter, (2021), arXiv:2101.11735 [astro-ph.CO] .
  * Graham  _et al._ [2015b] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Experimental Searches for the Axion and Axion-Like Particles, Ann. Rev. Nucl. Part. Sci. 65, 485 (2015b), arXiv:1602.00039 [hep-ex] .
  * Drukier  _et al._ [1986] A. K. Drukier, K. Freese, and D. N. Spergel, Detecting Cold Dark Matter Candidates, Phys. Rev. D 33, 3495 (1986).
  * Freese  _et al._ [1988] K. Freese, J. A. Frieman, and A. Gould, Signal Modulation in Cold Dark Matter Detection, Phys. Rev. D 37, 3388 (1988).
  * Vogelsberger  _et al._ [2009] M. Vogelsberger, A. Helmi, V. Springel, S. D. M. White, J. Wang, C. S. Frenk, A. Jenkins, A. D. Ludlow, and J. F. Navarro, Phase-space structure in the local dark matter distribution and its signature in direct detection experiments, Mon. Not. Roy. Astron. Soc. 395, 797 (2009), arXiv:0812.0362 [astro-ph] .
  * Kuhlen  _et al._ [2010] M. Kuhlen, N. Weiner, J. Diemand, P. Madau, B. Moore, D. Potter, J. Stadel, and M. Zemp, Dark Matter Direct Detection with Non-Maxwellian Velocity Structure, JCAP 02, 030, arXiv:0912.2358 [astro-ph.GA] .
  * Belokurov  _et al._ [2018] V. Belokurov, D. Erkal, N. W. Evans, S. E. Koposov, and A. J. Deason, Co-formation of the disc and the stellar halo, MNRAS 478, 611 (2018), arXiv:1802.03414 [astro-ph.GA] .
  * Helmi  _et al._ [2018] A. Helmi, C. Babusiaux, H. H. Koppelman, D. Massari, J. Veljanoski, and A. G. A. Brown, The merger that led to the formation of the Milky Way’s inner stellar halo and thick disk, Nature 563, 85 (2018), arXiv:1806.06038 [astro-ph.GA] .
  * Necib  _et al._ [2018a] L. Necib, M. Lisanti, and V. Belokurov, Inferred Evidence For Dark Matter Kinematic Substructure with SDSS-Gaia 10.3847/1538-4357/ab095b (2018a), arXiv:1807.02519 [astro-ph.GA] .
  * Necib  _et al._ [2018b] L. Necib, M. Lisanti, S. Garrison-Kimmel, A. Wetzel, R. Sanderson, P. F. Hopkins, C.-A. Faucher-Giguère, and D. Kereš, Under the Firelight: Stellar Tracers of the Local Dark Matter Velocity Distribution in the Milky Way 10.3847/1538-4357/ab3afc (2018b), arXiv:1810.12301 [astro-ph.GA] .
  * Evans  _et al._ [2019] N. W. Evans, C. A. J. O’Hare, and C. McCabe, Refinement of the standard halo model for dark matter searches in light of the Gaia Sausage, Phys. Rev. D 99, 023012 (2019), arXiv:1810.11468 [astro-ph.GA] .
  * O’Hare  _et al._ [2018] C. A. J. O’Hare, C. McCabe, N. W. Evans, G. Myeong, and V. Belokurov, Dark matter hurricane: Measuring the S1 stream with dark matter detectors, Phys. Rev. D 98, 103006 (2018), arXiv:1807.09004 [astro-ph.CO] .
  * O’Hare  _et al._ [2020] C. A. J. O’Hare, N. W. Evans, C. McCabe, G. Myeong, and V. Belokurov, Velocity substructure from Gaia and direct searches for dark matter, Phys. Rev. D 101, 023006 (2020), arXiv:1909.04684 [astro-ph.GA] .
  * Read  _et al._ [2008] J. I. Read, G. Lake, O. Agertz, and V. P. Debattista, Thin, thick and dark discs in $\Lambda$CDM, MNRAS 389, 1041 (2008), arXiv:0803.2714 [astro-ph] .
  * Danby and Camm [1957] J. M. A. Danby and G. L. Camm, Statistical dynamics and accretion, MNRAS 117, 50 (1957).
  * Danby and Bray [1967] J. M. A. Danby and T. A. Bray, Density of interstellar matter near a star, AJ 72, 219 (1967).
  * Griest [1988] K. Griest, Effect of the Sun’s Gravity on the Distribution and Detection of Dark Matter Near the Earth, Phys. Rev. D 37, 2703 (1988).
  * Sikivie and Wick [2002] P. Sikivie and S. Wick, Solar wakes of dark matter flows, Phys. Rev. D 66, 023504 (2002), arXiv:astro-ph/0203448 .
  * Alenazi and Gondolo [2006] M. S. Alenazi and P. Gondolo, Phase-space distribution of unbound dark matter near the sun, Phys. Rev. D 74, 083518 (2006).
  * Lee  _et al._ [2014] S. K. Lee, M. Lisanti, A. H. G. Peter, and B. R. Safdi, Effect of Gravitational Focusing on Annual Modulation in Dark-Matter Direct-Detection Experiments, Phys. Rev. Lett. 112, 011301 (2014), arXiv:1308.1953 [astro-ph.CO] .
  * Hu  _et al._ [2000] W. Hu, R. Barkana, and A. Gruzinov, Cold and fuzzy dark matter, Phys. Rev. Lett. 85, 1158 (2000), arXiv:astro-ph/0003365 .
  * Hui  _et al._ [2017] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Ultralight scalars as cosmological dark matter, Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] .
  * Bar-Or  _et al._ [2019] B. Bar-Or, J.-B. Fouvry, and S. Tremaine, Relaxation in a Fuzzy Dark Matter Halo, Astrophys. J. 871, 28 (2019), arXiv:1809.07673 [astro-ph.GA] .
  * Lancaster  _et al._ [2020] L. Lancaster, C. Giovanetti, P. Mocz, Y. Kahn, M. Lisanti, and D. N. Spergel, Dynamical Friction in a Fuzzy Dark Matter Universe, JCAP 01, 001, arXiv:1909.06381 [astro-ph.CO] .
  * Landau and Lifshits [1991] L. D. Landau and E. M. Lifshits,  _Quantum Mechanics: Non-Relativistic Theory_ , Course of Theoretical Physics, Vol. v.3 (Butterworth-Heinemann, Oxford, 1991).
  * Loudon [1973] R. Loudon,  _The quantum theory of light._ (1973).
  * Glauber [1963] R. J. Glauber, Coherent and incoherent states of the radiation field, Phys. Rev. 131, 2766 (1963).
  * Sudarshan [1963] E. C. G. Sudarshan, Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams, Phys. Rev. Lett. 10, 277 (1963).
  * Derevianko [2018] A. Derevianko, Detecting dark-matter waves with a network of precision-measurement tools, Phys. Rev. A 97, 042506 (2018), arXiv:1605.09717 [physics.atom-ph] .
  * Foster  _et al._ [2018] J. W. Foster, N. L. Rodd, and B. R. Safdi, Revealing the Dark Matter Halo with Axion Direct Detection, Phys. Rev. D 97, 123006 (2018), arXiv:1711.10489 [astro-ph.CO] .
  * Centers  _et al._ [2019] G. P. Centers _et al._ , Stochastic fluctuations of bosonic dark matter, (2019), arXiv:1905.13650 [astro-ph.CO] .
  * Lisanti  _et al._ [2021] M. Lisanti, M. Moschella, and W. Terrano, Stochastic properties of ultralight scalar field gradients, Phys. Rev. D 104, 055037 (2021), arXiv:2107.10260 [astro-ph.CO] .
  * Gramolin  _et al._ [2022] A. V. Gramolin, A. Wickenbrock, D. Aybas, H. Bekker, D. Budker, G. P. Centers, N. L. Figueroa, D. F. J. Kimball, and A. O. Sushkov, Spectral signatures of axionlike dark matter, Phys. Rev. D 105, 035029 (2022), arXiv:2107.11948 [hep-ph] .
  * Du  _et al._ [2018] N. Du _et al._ (ADMX), A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment, Phys. Rev. Lett. 120, 151301 (2018), arXiv:1804.05750 [hep-ex] .
  * Braine  _et al._ [2020] T. Braine _et al._ (ADMX), Extended Search for the Invisible Axion with the Axion Dark Matter Experiment, Phys. Rev. Lett. 124, 101303 (2020), arXiv:1910.08638 [hep-ex] .
  * Bartram  _et al._ [2021] C. Bartram _et al._ (ADMX), Search for ”Invisible” Axion Dark Matter in the $3.3\text{-}4.2~{}{\mu}$eV Mass Range, (2021), arXiv:2110.06096 [hep-ex] .
  * Brubaker  _et al._ [2017] B. M. Brubaker _et al._ , First results from a microwave cavity axion search at 24 $\mu$eV, Phys. Rev. Lett. 118, 061302 (2017), arXiv:1610.02580 [astro-ph.CO] .
  * Zhong  _et al._ [2018] L. Zhong _et al._ (HAYSTAC), Results from phase 1 of the HAYSTAC microwave cavity axion experiment, Phys. Rev. D 97, 092001 (2018), arXiv:1803.03690 [hep-ex] .
  * Kahn  _et al._ [2016] Y. Kahn, B. R. Safdi, and J. Thaler, Broadband and Resonant Approaches to Axion Dark Matter Detection, Phys. Rev. Lett. 117, 141801 (2016), arXiv:1602.01086 [hep-ph] .
  * Ouellet  _et al._ [2019] J. L. Ouellet _et al._ , First Results from ABRACADABRA-10 cm: A Search for Sub-$\mu$eV Axion Dark Matter, Phys. Rev. Lett. 122, 121802 (2019), arXiv:1810.12257 [hep-ex] .
  * Chaudhuri  _et al._ [2018] S. Chaudhuri, K. Irwin, P. W. Graham, and J. Mardon, Optimal Impedance Matching and Quantum Limits of Electromagnetic Axion and Hidden-Photon Dark Matter Searches, (2018), arXiv:1803.01627 [hep-ph] .
  * Kwon  _et al._ [2021] O. Kwon _et al._ (CAPP), First Results from an Axion Haloscope at CAPP around 10.7 $\mu$eV, Phys. Rev. Lett. 126, 191802 (2021), arXiv:2012.10764 [hep-ex] .
  * Caldwell  _et al._ [2017] A. Caldwell, G. Dvali, B. Majorovits, A. Millar, G. Raffelt, J. Redondo, O. Reimann, F. Simon, and F. Steffen (MADMAX Working Group), Dielectric Haloscopes: A New Way to Detect Axion Dark Matter, Phys. Rev. Lett. 118, 091801 (2017), arXiv:1611.05865 [physics.ins-det] .
  * McAllister  _et al._ [2017] B. T. McAllister, G. Flower, E. N. Ivanov, M. Goryachev, J. Bourhill, and M. E. Tobar, The ORGAN Experiment: An axion haloscope above 15 GHz, Phys. Dark Univ. 18, 67 (2017), arXiv:1706.00209 [physics.ins-det] .
  * Arvanitaki  _et al._ [2015] A. Arvanitaki, J. Huang, and K. Van Tilburg, Searching for dilaton dark matter with atomic clocks, Phys. Rev. D 91, 015015 (2015), arXiv:1405.2925 [hep-ph] .
  * Van Tilburg  _et al._ [2015] K. Van Tilburg, N. Leefer, L. Bougas, and D. Budker, Search for ultralight scalar dark matter with atomic spectroscopy, Phys. Rev. Lett. 115, 011802 (2015), arXiv:1503.06886 [physics.atom-ph] .
  * Hees  _et al._ [2016] A. Hees, J. Guéna, M. Abgrall, S. Bize, and P. Wolf, Searching for an oscillating massive scalar field as a dark matter candidate using atomic hyperfine frequency comparisons, Phys. Rev. Lett. 117, 061301 (2016), arXiv:1604.08514 [gr-qc] .
  * Aharony  _et al._ [2021] S. Aharony, N. Akerman, R. Ozeri, G. Perez, I. Savoray, and R. Shaniv, Constraining Rapidly Oscillating Scalar Dark Matter Using Dynamic Decoupling, Phys. Rev. D 103, 075017 (2021), arXiv:1902.02788 [hep-ph] .
  * Kennedy  _et al._ [2020] C. J. Kennedy, E. Oelker, J. M. Robinson, T. Bothwell, D. Kedar, W. R. Milner, G. E. Marti, A. Derevianko, and J. Ye, Precision Metrology Meets Cosmology: Improved Constraints on Ultralight Dark Matter from Atom-Cavity Frequency Comparisons, Phys. Rev. Lett. 125, 201302 (2020), arXiv:2008.08773 [physics.atom-ph] .
  * Savalle  _et al._ [2021] E. Savalle, A. Hees, F. Frank, E. Cantin, P.-E. Pottie, B. M. Roberts, L. Cros, B. T. Mcallister, and P. Wolf, Searching for Dark Matter with an Optical Cavity and an Unequal-Delay Interferometer, Phys. Rev. Lett. 126, 051301 (2021), arXiv:2006.07055 [gr-qc] .
  * Campbell  _et al._ [2021] W. M. Campbell, B. T. McAllister, M. Goryachev, E. N. Ivanov, and M. E. Tobar, Searching for Scalar Dark Matter via Coupling to Fundamental Constants with Photonic, Atomic and Mechanical Oscillators, Phys. Rev. Lett. 126, 071301 (2021), arXiv:2010.08107 [hep-ex] .
  * Antypas  _et al._ [2019] D. Antypas, O. Tretiak, A. Garcon, R. Ozeri, G. Perez, and D. Budker, Scalar dark matter in the radio-frequency band: atomic-spectroscopy search results, Phys. Rev. Lett. 123, 141102 (2019), arXiv:1905.02968 [physics.atom-ph] .
  * Oswald  _et al._ [2021] R. Oswald _et al._ , Search for oscillations of fundamental constants using molecular spectroscopy, (2021), arXiv:2111.06883 [hep-ph] .
  * Arvanitaki  _et al._ [2018] A. Arvanitaki, P. W. Graham, J. M. Hogan, S. Rajendran, and K. Van Tilburg, Search for light scalar dark matter with atomic gravitational wave detectors, Phys. Rev. D 97, 075020 (2018), arXiv:1606.04541 [hep-ph] .
  * Grote and Stadnik [2019] H. Grote and Y. V. Stadnik, Novel signatures of dark matter in laser-interferometric gravitational-wave detectors, Phys. Rev. Res. 1, 033187 (2019), arXiv:1906.06193 [astro-ph.IM] .
  * Vermeulen  _et al._ [2021] S. M. Vermeulen _et al._ , Direct limits for scalar field dark matter from a gravitational-wave detector, (2021), arXiv:2103.03783 [gr-qc] .
  * Aiello  _et al._ [2021] L. Aiello, J. W. Richardson, S. M. Vermeulen, H. Grote, C. Hogan, O. Kwon, and C. Stoughton, Constraints on scalar field dark matter from co-located Michelson interferometers, (2021), arXiv:2108.04746 [gr-qc] .
  * Badurina  _et al._ [2022] L. Badurina, D. Blas, and C. McCabe, Refined ultralight scalar dark matter searches with compact atom gradiometers, Phys. Rev. D 105, 023006 (2022), arXiv:2109.10965 [astro-ph.CO] .
  * Graham and Rajendran [2013] P. W. Graham and S. Rajendran, New Observables for Direct Detection of Axion Dark Matter, Phys. Rev. D 88, 035023 (2013), arXiv:1306.6088 [hep-ph] .
  * Budker  _et al._ [2014] D. Budker, P. W. Graham, M. Ledbetter, S. Rajendran, and A. Sushkov, Proposal for a Cosmic Axion Spin Precession Experiment (CASPEr), Phys. Rev. X 4, 021030 (2014), arXiv:1306.6089 [hep-ph] .
  * Bloch  _et al._ [2020] I. M. Bloch, Y. Hochberg, E. Kuflik, and T. Volansky, Axion-like Relics: New Constraints from Old Comagnetometer Data, JHEP 01, 167, arXiv:1907.03767 [hep-ph] .
  * Bloch  _et al._ [2021] I. M. Bloch, G. Ronen, R. Shaham, O. Katz, T. Volansky, and O. Katz (NASDUCK), NASDUCK: New Constraints on Axion-like Dark Matter from Floquet Quantum Detector (2021) arXiv:2105.04603 [hep-ph] .
  * Abel  _et al._ [2017] C. Abel _et al._ , Search for Axionlike Dark Matter through Nuclear Spin Precession in Electric and Magnetic Fields, Phys. Rev. X 7, 041034 (2017), arXiv:1708.06367 [hep-ph] .
  * Foster  _et al._ [2021] J. W. Foster, Y. Kahn, R. Nguyen, N. L. Rodd, and B. R. Safdi, Dark Matter Interferometry, Phys. Rev. D 103, 076018 (2021), arXiv:2009.14201 [hep-ph] .
  * Chen  _et al._ [2021] Y. Chen, M. Jiang, J. Shu, X. Xue, and Y. Zeng, Dissecting Axion and Dark Photon with A Network of Vector Sensors, (2021), arXiv:2111.06732 [hep-ph] .
  * Arkani-Hamed  _et al._ [2009] N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, and N. Weiner, A Theory of Dark Matter, Phys. Rev. D 79, 015014 (2009), arXiv:0810.0713 [hep-ph] .
  * Schönrich  _et al._ [2010] R. Schönrich, J. Binney, and W. Dehnen, Local kinematics and the local standard of rest, MNRAS 403, 1829 (2010), arXiv:0912.3693 [astro-ph.GA] .
  * Eilers  _et al._ [2019] A.-C. Eilers, D. W. Hogg, H.-W. Rix, and M. K. Ness, The Circular Velocity Curve of the Milky Way from 5 to 25 kpc, ApJ 871, 120 (2019), arXiv:1810.09466 [astro-ph.GA] .
  * Abadi  _et al._ [2003] M. G. Abadi, J. F. Navarro, M. Steinmetz, and V. R. Eke, Simulations of Galaxy Formation in a $\Lambda$ Cold Dark Matter Universe. II. The Fine Structure of Simulated Galactic Disks, ApJ 597, 21 (2003), arXiv:astro-ph/0212282 [astro-ph] .
  * Kazantzidis  _et al._ [2008] S. Kazantzidis, J. S. Bullock, A. R. Zentner, A. V. Kravtsov, and L. A. Moustakas, Cold Dark Matter Substructure and Galactic Disks. I. Morphological Signatures of Hierarchical Satellite Accretion, ApJ 688, 254 (2008), arXiv:0708.1949 [astro-ph] .
  * Villalobos and Helmi [2008] Á. Villalobos and A. Helmi, Simulations of minor mergers - I. General properties of thick discs, MNRAS 391, 1806 (2008).
  * Rix and Bovy [2013] H.-W. Rix and J. Bovy, The Milky Way’s stellar disk. Mapping and modeling the Galactic disk, A&A Rev. 21, 61 (2013), arXiv:1301.3168 [astro-ph.GA] .
  * Soubiran  _et al._ [2003] C. Soubiran, O. Bienaymé, and A. Siebert, Vertical distribution of Galactic disk stars. I. Kinematics and metallicity, A&A 398, 141 (2003), arXiv:astro-ph/0210628 [astro-ph] .
  * Purcell  _et al._ [2009] C. W. Purcell, J. S. Bullock, and M. Kaplinghat, The Dark Disk of the Milky Way, ApJ 703, 2275 (2009), arXiv:0906.5348 [astro-ph.GA] .
  * Helmi [2020] A. Helmi, Streams, Substructures, and the Early History of the Milky Way, ARA&A 58, 205 (2020), arXiv:2002.04340 [astro-ph.GA] .
  * Wigner [1932] E. P. Wigner, On the quantum correction for thermodynamic equilibrium, Phys. Rev. 40, 749 (1932).
  * Binney and Tremaine [2008] J. Binney and S. Tremaine,  _Galactic Dynamics: Second Edition_ (2008).