# Riding the dark matter wave: Novel limits on general dark photons from LISA
Pathfinder

Jonas Frerick<EMAIL_ADDRESS>Joerg Jaeckel<EMAIL_ADDRESS>heidelberg.de Felix Kahlhoefer<EMAIL_ADDRESS>Kai Schmidt-Hoberg
<EMAIL_ADDRESS>
###### Abstract

We note the possibility to perform a parametrically improved search for gauged
baryon ($B$) and baryon minus lepton ($B-L$) Dark Photon Dark Matter (DPDM)
using auxiliary channel data from LISA Pathfinder. In particular we use the
measurement of the differential movement between the test masses (TMs) and the
space craft (SC) which is nearly as sensitive as the tracking between the two
TMs. TMs and SC are made from different materials and therefore have different
charge-to-mass ratios for both $B-L$ and $B$. Thus, the surrounding DPDM field
induces a relative acceleration of nearly constant frequency. For the case of
$B-L$, we find that LISA Pathfinder can constrain previously unexplored
parameter space, providing the world leading limits in the mass range $4\cdot
10^{-19}\,\text{eV}<m<3\cdot 10^{-17}\,\text{eV}$. This limit can easily be
recast also for dark photons that arise from gauging other global symmetries
of the SM.

###### keywords:

dark photon dark matter , direct detection , gravitational wave interferometry

††journal: Physics Letters B

[desy]organization= Deutsches Elektronen-Synchrotron DESY, addressline=
Notkestr. 85, postcode=22607, city=Hamburg, country=Germany

[heidelberg]organization=Institut für Theoretische Physik, Universität
Heidelberg, addressline= Philosophenweg 16, postcode=69120 ,city=Heidelberg,
country=Germany

[kit]organization=Institute for Theoretical Particle Physics (TTP), Karlsruhe
Institute of Technology (KIT), postcode= 76128,city =Karlsruhe,
country=Germany

##  1 Introduction

The existence of Dark Matter (DM) is a well-established observational fact
[1]. For a long time the WIMP paradigm has dominated the quest for DM [2] but
with null observations in the increasingly sensitive direct detection (DD)
experiments [3, 4, 5] there is increased interest in alternative DM
candidates. One particularly well-motivated class of such particles are ultra-
light and weakly coupled bosons (see, e.g., [6, 7, 8] for reviews). These
include, among others, the axion [9, 10, 11, 12] or general axion-like
particles (ALPs) [13] as well as new vector bosons (cf., e.g., [14, 15, 16,
17]), often referred to as dark photons (DPs).111In this work we will refer to
any light new vector boson as a dark photon, allowing for couplings that gauge
global symmetry groups of the SM. Our method is unfortunately insensitive to
the “canonical”, kinetically mixed DP.

In this work, we will focus on ultra-light DPs as a DM candidate. Small DP
masses $m$ can be generated either by the Stückelberg [18, 19] or by the Higgs
mechanism (where the latter often causes additional constraints). For
$m\lesssim 30\,$eV the local DM halo behaves like a classical wave as the
spacing between particles becomes smaller than the de Broglie wave length
[20]. We will stay agnostic about the details of the production mechanism.
While heavier DPs can be produced from the thermal SM bath [21], very light
DPs require a non-thermal mechanism to ensure that the DM is cold. To this
end, there are many gravitational or extended dark sector solutions that
provide the correct relic density [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32,
33, 34, 35]. Furthermore, we will assume that the DP is gauged under a
combination of baryon number ($B$) and lepton number ($L$), with a particular
focus on the difference ($B-L$).

Astrophysical objects like the Sun can efficiently produce light DPs which, in
turn, results in impressive limits on the existence of DPs without requiring
them to be DM [36, 37, 38, 39, 40]. Additionally, lab experiments testing the
equivalence principle are perfect candidates to look for this kind of new
physics that induces long-range forces beyond electromagnetism and gravity
[41, 42]. Planetary [43] and asteroidal [44, 45] orbits are sensitive to new
long-range forces as well. Finally, the gauge anomaly associated with baryon
number leads to strong constraints from meson decays [15] for this specific
gauge group.

If the DPs are also DM, new tests become available, e.g. with accelerometers
as proposed in Ref. [46] and realized in [47]. A particularly interesting
possibility is to search for these DM candidates directly at gravitational
wave observatories using the method outlined in Ref. [48]. This original
analysis also focused on new vector bosons gauged under $B$ and $B-L$, and is
applicable to both ground-based and space-based laser interferometers.
However, it makes use of the differential acceleration between the two equal
test masses and is therefore limited by the arm length. Here, we point out
that the setup of LISA Pathfinder (LPF) also offers the possibility to use the
differential acceleration between a test mass and the satellite carrying the
interferometer itself. As we will argue this type of search is not limited by
the arm length, significantly increasing the sensitivity in the relevant mass
range compared to a previous analysis [49]. That said, we want to point out
that the use of auxiliary channels was already proposed for KAGRA in Ref.
[50].

We briefly discuss the signal prediction and analysis method in sec. 2 and
point out the similarities and differences to previous DPDM interferometer
limits. This is followed by an introduction to LPF including a discussion of
the sensitivity in sec. 3. Finally, we estimate the sensitivity of this
instrument to DPDM and discuss the necessary steps for a refined analysis in
sec. 4 before we conclude in sec. 5. Throughout this letter we work in natural
units $\hbar=c=1$.

##  2 Calculation of the signal

Let us begin by introducing the DP Lagrangian [51, 52, 16]

$\mathcal{L}\supset-\frac{1}{4}F_{\mu\nu}^{\prime}F^{\prime\mu\nu}-\frac{\epsilon_{\text{KM}}}{2}F_{\mu\nu}^{\prime}F^{\mu\nu}+\frac{m^{2}}{2}A_{\mu}^{\prime}A^{\prime\mu}-\epsilon_{g}eA_{\mu}^{\prime}J^{\mu}_{g}\;,$
(1)

which contains the renormalisable interactions of the DP field
$A^{\prime}_{\mu}$ in full generality. It takes into account both the kinetic
mixing $\epsilon_{\text{KM}}$ between the field strength tensors
$F_{\mu\nu}^{(\prime)}$ of the SM photon and the dark photon, and an explicit
coupling $\epsilon_{g}$ to a current $J^{\mu}_{g}$ associated with a gauge
group $g$. Note that we have rescaled the gauge coupling $g_{g}$ to the
electromagnetic coupling $e$, i.e. $g_{g}=\epsilon_{g}e$. From now on, we will
assume the kinetic mixing to be negligible and focus on the explicit couplings
with $g=B$ or $g=B-L$ unless mentioned otherwise. In sec. 4 we will discuss
how to generalize our analysis to arbitrary gauge groups.

Under the aforementioned assumptions, any piece of baryonic matter is directly
charged under both gauge groups. Thus, in a background field of DPDM these
charges behave in full analogy to electric charges in an electric field.
Assuming the DM to be cold and have mass $m$, the field will be nearly
monochromatic with a linewidth suppressed by the non-relativistic velocity
$v\sim 10^{-3}$ of the halo [53]

$\mathbf{A}(t,x)=\mathbf{A}_{\text{DM}}e^{-i\omega t+\phi(x)}\;,$ (2)

where $\omega=m+\mathcal{O}(v^{2})$ denotes the non-relativistic particle
energy, $\mathbf{A}_{\text{DM}}$ is the 3-vector of the DPDM field, and
$\phi(x)=i\;\mathbf{k}\cdot\mathbf{x}+\phi_{0}$ is a weakly position dependent
phase where $\mathbf{k}\approx m\mathbf{v}$ denotes the momentum of the wave
and $\phi_{0}$ is a constant phase.222For our purposes, weakly dependent means
that $|\mathbf{k}|L\ll 1$ where $L$ denotes the size of the experiment. We can
obtain the temporal component of the 4-potential from the “Lorenz condition”

$\partial_{\mu}A^{\mu}=0\Rightarrow
A^{0}(t)=-\frac{\mathbf{k}\cdot\mathbf{A}(t)}{\omega}\approx\mathbf{v}\cdot\mathbf{A}(t)\ll|\mathbf{A}(t)|\;,$
(3)

which has to be fulfilled for a massive vector boson as dictated by the
equations of motion. Due to the weak spatial dependence of eq. (2) we have
omitted the $x$ in the argument of all vector components. We observe that the
temporal component is generically velocity suppressed for non-relativistic
DPDM.333Indeed, for transversely polarized DPs, i.e.
$\mathbf{k}\cdot\mathbf{A}_{\text{DM}}=0$, the component vanishes exactly.

Within a coherence patch the signal can be treated as being monochromatic. The
coherence length is given by the wavelength $\lambda_{c}\simeq 2\pi/(mv)$ and
the coherence time is $t_{c}\simeq 2\pi/(mv^{2})$. It is important to realize
that neither the amplitude nor the direction of the field changes within a
coherence patch. LPF covers a frequency range from a few Hz down to around
$10^{-5}\,$Hz (see sec. 3). Especially for the lowest frequencies this implies
an extremely long coherence time due to the non-relativistic velocities. In
fact, even for the highest frequencies in LPF’s sensitivity range the
coherence time is more than a week so that even long-term searches for
monochromatic signals suffer at most weakly from the decoherence of the
signal, thus enhancing the limits significantly without employing new
technology.

For the sake of clarity, we will give a minimal description of LPF here and
use it as an example for the main idea behind our analysis but the following
arguments are equally valid for two generic objects made from different
materials. In LPF, two (almost) identical test masses (TMs) are enclosed
separately in a space craft (SC) and the relative motion between the TMs
themselves and between TMs and the SC is tracked. A more detailed description
will follow in sec. 3. The dominant effect of the gauged DPDM on a charged
object is analogous to the electric component of the Lorentz force (see [46]
for a similar calculation as we do in the following). Therefore, we need to
determine the “electric field”

$\displaystyle\mathbf{E}_{g}=-\partial_{t}\mathbf{A}(t)=i\omega\mathbf{A}_{\text{DM}}e^{-i\omega
t+\phi(x)}\;.$ (4)

For our purposes, we can completely ignore the phase of the field
$\phi(x)\approx 0$ and consider it spatially constant over the size of the
experiment. Therefore, the field is oscillating at a single frequency as long
as we consider only coherent time scales.

This field then exerts a force on all objects charged under the given gauge
group. Therefore, we find the following acceleration

$\displaystyle\mathbf{a}(t)\simeq
i\omega\epsilon_{g}e\frac{q}{M}\mathbf{A}_{\text{DM}}e^{-i\omega
t}=i\frac{q}{M}\sqrt{2\rho_{\text{DM}}}\;\hat{\mathbf{e}}_{A}\;e^{-i\omega
t}\;,$ (5)

for an object of mass $M$ and charge $q$ (under $g$). Furthermore, we used the
well-known relation between the average energy density and amplitude of wave-
like DM

$\rho_{\text{DM}}=\frac{1}{2}\omega^{2}|\mathbf{A}_{\text{DM}}|^{2}\;.$ (6)

Finally, $\hat{\mathbf{e}}_{A}$ is the unit vector in direction of the DPDM
field within a coherence patch.

To estimate the signal-to-noise ratio in LPF we need the amplitude of the
relative acceleration between the SC center of mass and the TMs for the three
SC axes $i=x,y,z$. For a monochromatic signal, we obtain those values by
taking the real part and dropping the harmonic behavior of eq. (5)

$\displaystyle\Delta a_{i}=$
$\displaystyle\epsilon_{g}e\left(\Delta\frac{q}{M}\right)\sqrt{2\rho_{\text{DM}}}\cos\theta_{A,i}\;.$
(7)

We note that a one-dimensional setup or even a typical planar interferometer
can in principle be totally insensitive to this effect if the polarization is
orthogonal to the plane of the experiment. LPF offers the advantage that it
features a full 3D sensitivity of the SC motion w.r.t. the TMs.

In general, there a two different polarization models in the literature. One
of them assumes that the DPs have the same polarization everywhere in space
while the other one assumes that each coherence patch has a different
polarization which is distributed uniformly on the unit sphere. From our
discussion of the coherence time we conclude that in the ultra-low mass regime
we will not be able to tell the difference as all experiments with realistic
lifetimes will only observe a single coherence patch.444At the higher end of
the frequency range the situation might be more promising if we allow for an
observation time of several years. This can get even better if the sensitivity
can be extended to higher frequencies. In contrast, for larger masses
measuring at different times corresponds to measuring different polarization
of the DPs. Combining this with the fact that LPF has a non-trivial orbit and
orientation will result in a very complex scheme required to perform a
rigorous analysis. Nevertheless, applying this information which in principle
is known might provide additional constraining power as demonstrated in Ref.
[54]. We will treat this issue in more detail in sec. 4. To conclude this
discussion we emphasize that the position and orientation of the SC will not
change significantly on the time scales of the used observations [55].

Let us quickly compare our result, using the auxiliary channels between SC and
TMs, to the case where $\Delta\frac{q}{M}=0$ which corresponds to two bodies
made from the same material. To the best of our knowledge, the TMs for the
main interferometers in all GW searches including LPF fulfill this criterion.
Additionally, any elemental impurities that could break this degeneracy are
kept extremely small in order to improve the performance of the
interferometer. Therefore, in this case, we have to look for subleading
effects, e.g. from the phase in eq. (4) which introduces both an arm-length
and a velocity suppression. This was the standard approach in the literature
following Ref. [48] until it was refined later using the finite light-
traveling time of the laser [56] which leads to an improvement if the length
scale associated with the DP mass $1/m$ coincides with the arm length of the
interferometer $L$. This analysis method and the originally proposed
decoherence effect from the small inhomogeneity of the field will give an
observable relative acceleration even for strictly equal charge-to-mass
ratios. Nevertheless, this acceleration is suppressed by $\max\left\\{(\omega
L)^{2},v\omega L\right\\}$. For full scale interferometers where the arm
length is on the scale of $1/m\sim 1/\omega$ by construction this may not be
big a problem. But, for LPF with its very limited arm length of $\sim 40\,$cm
these effects will substantially suppress all limits derived following the
standard methods in the literature as shown in Ref. [49]. Therefore, looking
for auxiliary channels between TMs and the SC that feature
$\Delta\frac{q}{M}\neq 0$ is promising. Indeed, for a different gravitational
wave interferometer, KAGRA, this observation was already utilized in Ref.
[50], which enhanced the limits in the low frequency region significantly.

With the prediction for the acceleration amplitude, it is straightforward to
estimate the signal-to-noise ration (SNR) of an interferometer with a given
relative acceleration amplitude spectral density (ASD) $S_{a}^{1/2}(f)$ via

$\text{SNR}=\frac{\Delta a_{i}}{S_{a}^{1/2}(f)}\sqrt{T_{\text{eff}}}\;,$ (8)

where $T_{\text{eff}}$ depends on observation time $T_{\text{obs}}$ and
coherence time $t_{c}$ via

$\displaystyle T_{\text{eff}}=\begin{cases}T_{\text{obs}}\ \ ,&\quad
T_{\text{obs}}\leq t_{c}\\\ \sqrt{T_{\text{obs}}t_{c}}\ \ ,&\quad
T_{\text{obs}}>t_{c}\end{cases}\;,$ (9)

as outlined e.g. in the appendix of Ref. [57].555Conventionally, results are
quoted as ASDs for the relative displacement instead of acceleration as used
by us. The translation from our ASD to this so-called strain sensitivity is
fairly simple as it just requires a rescaling factor of $\sim\omega^{2}/L$.

##  3 LISA Pathfinder sensitivity

LPF [58] was a precursor mission to the planned space-borne gravitational wave
interferometer LISA [59]. The mission’s objective was to demonstrate that the
technology developed for LISA will be able to perform as predicted under
realistic space conditions. For this purpose, a SC containing two TMs was sent
to the first Lagrange point of the Sun-Earth system. These TMs are $2\,$kg
Gold-Platinum alloy cubes with a side length of $\sim 5\,$cm and they were
placed in two separate electrode housings with an optical bench placed in
between. The main aim was to keep the noise in the relative acceleration
between the two free falling TMs at a level that would verify the
applicability of this technology for LISA. Indeed, the test was successful and
performed even better than expected [60, 61]. Such a high-precision instrument
requires more scrutiny than just a single interferometer measuring the
relative TM displacement. Therefore, LPF contained a radiation monitor [62],
additional interferometers [63] and capacitive sensing [64]. Several of these
auxiliary channels are used to avoid a collision between “TM1”, the reference
test mass, and the SC. The choice of a preferred TM is required as two TMs on
their respective geodesics within a single SC cannot coexist without a
collision. Therefore, there are measures in place to correct the trajectory of
the second TM w.r.t. to the reference TM.

The TMs are aligned on what the collaboration labeled the x-axis. On this
axis, there is the so called $x_{12}$ interferometer which is the central
instrument on-board as it is used to measure the relative acceleration between
the TMs. For our purposes, we want to focus on another instrument, the $x_{1}$
interferometer, which controls the SC position w.r.t. to TM1. This auxiliary
interferometer will be the best channel to search for DPDM over a large
mass/frequency range.

As we have pointed out in sec. 2 we require knowledge about the charge-to-mass
ratio and therefore the elemental composition of both the TMs and the SC.
Unfortunately, the SC itself is made up of a collection of different materials
and components but they are on average expected to be at much lower atomic
number than Au or Pt and thus they will have different charge-to-mass ratios.
This result is intuitive for $B-L$ as the total charge of an atom is given by
the neutron number and the neutron-to-proton ratio tends to increase with
atomic number yielding different charge-to-mass ratios for light and heavy
elements. This effect is more subtle for $B$. The difference mainly comes from
the variation in binding energies and the small mass difference between proton
and neutron. This immediately provides us with an estimate for the suppression
of the charge-to-mass ratio w.r.t. to the $B-L$ result: both the binding
energies and the nucleon mass difference are of order MeV compared to the
total nucleon masses which are at the GeV scale. Naively, this suggests a
suppression factor $\sim 10^{-3}$ which turns out to be quite accurate, as we
will see in Sec. 4.

In the composition of the SC, the second most important contribution after the
technology package enclosing the TMs arises from the structure of the SC which
is made mostly from carbon and aluminium [65]. Indeed, the SC contains many
different sub-components of similar mass and some of them will also contain
elements with atomic number much larger than C or Al. A detailed analysis of
the SC composition is beyond the scope of this letter and thus we will simply
use a lower bound on the charge-to-mass ratio of the SC. To arrive at this
conservative estimate, we will assume that all components are made from the
same material as the TMs except for the SC structure which we assume to be
entirely made from carbon. Using this approximation and table 1 from Ref.
[65], we conclude that the $450\,$kg SC has an $83\,$kg C component and the
remaining material will have a charge-to-mass ratio equal to that of gold.

For a better understanding of the geometry we show an exploded view of LPF in
fig. 1. An important factor in the sensitivity analysis of the SC motion
against the TMs is that not just the x-axis but also the y- and z-axis are
tracked where the z-axis points from the TMs to the solar array. These axes
are measured via capacitive sensing which in general is less precise than the
interferometers for most frequencies. Nevertheless, we have the advantage of
being able to analyze the relative acceleration ASDs for all SC axes [66] and
we show these results in fig. 2. They represent the simulated, data-backed
sensitivities to the relative acceleration of the SC w.r.t. to the TM(s) which
is exactly what we are interested in for eq. (8).666In fact, the y- and
z-direction is tracked w.r.t. to the average of both TM coordinates while for
the x-axis only the relative motion w.r.t. TM1 is measured via the additional
interferometer. These results were obtained from a 6.5 day noise-only run in
April 2016 [60]. The curves explicitly account for all known noise on SC and
TMs and therefore they present the best estimate for the stability of the SC
w.r.t. to the TM(s). To derive limits, we will set cuts at $1\,$Hz and
$10^{-4}\,$Hz as a careful evaluation of the highest and lowest frequencies is
beyond the scope of this work. Nevertheless, a detailed analysis of the data
will most likely lead to interesting constraints in these extremal regimes.

Figure 1: Exploded view of LPF showing the science module containing the test
masses in its center as well as the propulsion module. Image by ESA/ATG
medialab (with permission).

Figure 2: Sensitivity of the LPF SC acceleration w.r.t. the TM(s). Red shows
the interferometer sensitivity while green and orange are found from the
capacitive sensors in the housing averaged over both TMs; data from [66]. The
dashed black line shows the maximum sensitivity of the LPF TM-TM measurement
based on the same data set; taken from [60]

.

Before we calculate the limits from these ASDs, let us briefly discuss their
behavior. For higher frequencies down to $\sim 10^{-3}\,$Hz, the sensitivity
is limited mostly by the so-called out-of-loop noise which describes several
external influences on the SC. The bump at the low frequency part of the
spectrum is due to the star-tracker noise which comes from imperfections in
the determination of the position of the SC. Ref. [66] argues that this low-
frequency noise will most likely be mitigated in the LISA mission pointing out
an interesting avenue for future investigations of gauged DPDM. We will
discuss this in more detail in sec. 4. In the extreme low frequency region we
observe additional loss in sensitivity from the capacitive actuation noise
experienced by the TMs. In this regime, the simulation also predicts a
significantly better sensitivity than the data as shown in Ref. [66] further
justifying the cuts introduced above. At peak sensitivity, the x-axis almost
reaches the TM1-TM2 result, cf. the dashed black line in fig. 2 taken from
Ref. [60] which is based on the same data sample.777Unfortunately, the full
frequency range is not shown in that work. At this point of best sensitivity
the other axes perform comparably worse as the capacitive sensing cannot
compete with the interferometer on the x-axis.

##  4 Results

Material |  | Au | C | SC-TM
---|---|---|---|---
$\left(\frac{q}{M}\right)_{B-L}\ \text{in GeV}^{-1}$ |  | $0.64$ | $0.54$ | $0.018$
$\left(\frac{q}{M}\right)_{B}\ \text{in GeV}^{-1}$ |  | $1.0736$ | $1.0737$ | $1.8\cdot 10^{-5}$
Table 1: Charge-to-mass ratios for Au, C, and the difference between these two
elements rescaled to our estimate for the SC composition.

Now let us piece together our detailed knowledge of the LPF sensitivity with
our signal prediction. Table 1 shows the elemental charge-to-mass ratios for
carbon and gold [67]. We ignore the Pt contribution to the TMs as its charge-
to-mass ratio is close to the one of Au. We can derive the simple relation for
the SC-TM difference of charge-to-mass ratio under the conservative
assumptions about the SC composition of sec. 3:

$\displaystyle\left|\Delta\left(\frac{q}{M}\right)\right|=$
$\displaystyle\left|\left(\frac{q}{M}\right)_{\text{TM}}-\left(\frac{q}{M}\right)_{\text{SC}}\right|\approx
0.18\left|\left(\frac{q}{M}\right)_{\text{Au}}-\left(\frac{q}{M}\right)_{\text{C}}\right|\;.$
(10)

The last column of Table 1 shows the corresponding absolute value of this
difference.888It is an unfortunate coincidence that the baryon charge-to-mass
ratio is so similar for Au and C. Taking into account the true composition of
the SC will alleviate this suppression. Finally, we observe that our initial
estimate of the suppression in the $B$ charge-to-mass ratio is in good
agreement with the actual calculation.

Demanding that the SNR in eq. (8) is at most unity we get a good estimate for
the LPF sensitivity on the coupling strength of the DP to the chosen gauge
group. For the DM density, we assume $\rho_{\text{DM}}\simeq 0.4\,$GeV/cm3
[68].

As noted earlier, there is a rigorous way to combine the different axes but it
requires taking into account a proper convolution of SC position and
orientation with all possible DP polarizations. Whereas this procedure
necessitates knowledge of the exact orbit of LPF it will provide even stronger
limits if one follows the detailed guide provided in Ref. [54]. LPF offers the
advantageous feature that it is sensitive in all three spatial dimensions
which means that our results cannot suffer from a “blindness” due to an
unfortunate orientation of the polarization. We can always set a conservative
estimate from taking the least sensitive axis for every frequency according to
fig. 2.

In fig. 3, we show our main result for $B-L$ as solid lines for the individual
axes following the color-coding of fig. 2. Here we assume for each axis
separately that the polarization is exactly aligned with the given axis, i.e.
setting $\cos\theta_{A,i}=1$ in eq. (7). We see that we get rather similar
constraints from all axes except for the better peak sensitivity of the
x-axis. Following the above argument by taking the upper envelope, i.e. to
just consider the weakest limit for every mass, we can obtain a conservative
combination of the limits.

Keeping this in mind, we will nevertheless opt for a more optimistic way to
simplify the visualization of additional forecasts and later the results for
$B$. Following Ref. [48] we perform an average over all possible velocities
and polarizations. While this is technically not the most conservative
assumption for these long coherence times, we adopt this approach to
facilitate comparison with previous LPF limits [49] and LISA projections [48,
56] (see also appendix A of Ref. [69]). As our limits are independent of the
velocity of the DPs, the resulting “geometry factor” is $1/\sqrt{3}$ as
compared to the usual result of $1/3$. Then, instead of taking the upper
envelope we will use the lower envelope, i.e. the strongest limit for every
mass, multiplied by this suppression factor. We will refer to this as the
envelope simplification.999Most of the limits we show with this method are
more optimistic projections anyways. Only for the solid blue line of the $B$
limits in fig. 4 and the blue region in fig. 5 we should keep the shape of all
three axes in mind.

Using this approach we also include an estimate of the improved reach of LPF
as a dashed blue line taking into account the whole data set and using the
improved understanding of the detector noise and reduction of Brownian noise
in the later stages of the mission [61].101010For this we set the observation
time to 1 month and assume an improvement in noise-reduction by a factor 3. We
demonstrate the impact of the observation time by also including the blue
dotted line which assumes the same sensitivity as the solid lines but we set
the observation time to the coherence time for each frequency. This explains
why the high frequency sensitivity is similar to the "fixed observation time"
scenario: for the highest frequencies available, i.e. around $1\,$Hz, the
coherence time is about $10^{6}\,$s which is roughly on the time scale of a
week, coinciding with the real observation time used to model the sensitivity
curve. On the low frequency end of the spectrum, the coherence times approach
millennia scales making the dashed line much stronger than the other limits.

To demonstrate the power of our approach we compare it to three major results
from the literature. We see that our analysis is able to cover new parameter
space beyond the otherwise dominant limits set by the fifth-force search
interpretation of the MICROSCOPE experiment [42]. Furthermore, it is
immediately clear that our analysis can easily outperform previous LPF limits
[49] just because there is no need to rely on the decoherence of the field
which is the dominant effect if one only considers the two test masses. In
fact, the improvement of our limits over the naive results that can be
obtained from the decoherence method evaluated around our peak sensitivity at
$\sim 5\cdot 10^{-18}\,$eV is given by

$\displaystyle\frac{\epsilon_{B-L,\text{sat}}}{\epsilon_{B-L,\text{dec}}}\sim\frac{\Delta\frac{q}{m}}{\frac{q}{M}}\left(mvL\right)^{-1}\sim
3\cdot 10^{12}\;,$ (11)

ignoring the small difference in sensitivity between the $x_{12}$ and the
$x_{1}$ sensitivity at this mass. The first factor takes into account that our
method suffers from a mild charge-to-mass ratio suppression w.r.t. the
decoherence method whereas the second factor comes from smallness of the
decoherence on a length scale of $40\,$cm. We note that the limits found in
Ref. [49] are better than naively expected from our analysis method,
presumably because of their more sophisticated statistical analysis. This
observation makes us confident in the potential reach of our approach for
future analyses using all the available data. The third literature result is a
LISA forecast using the conventional analysis method [56].

Figure 3: Limits on the rescaled coupling to $B-L$, $\epsilon_{B-L}$, of DPDM.
In grey we show the DM-independent limits from searches for violation of the
Equivalence Principle [42] and the dark red filled region shows the LPF limits
derived from decoherence in Ref. [49]. The dark red dotted line shows the
forecast from LISA [56]. In red, green and orange we show the main result of
this paper. Forecasts for similar analyses are shown in blue using the
envelope simplification explained in the text.

Before turning to the LISA projections let us discuss the results for $B$
shown in fig. 4. We lose around 3 orders of magnitude in sensitivity which can
be explained by the stronger charge-to-mass ratio suppression in eq. (11) for
$B$. This becomes immediately clear in the comparison of our results to the
decoherence limits from LPF which do not suffer from this issue as they scale
with the total charge-to-mass ratio. Nevertheless, the decrease in sensitivity
of the Equivalence Principle limits due to the same effect still allows to
probe a small region of new parameter space and makes an extended study of LPF
(and LISA) auxiliary channels very attractive as it will cover a significant
amount of new parameter space. Additionally, we added limits from the baryon
number anomaly [15] which are non-existent for $B-L$.

Figure 4: Limits on the rescaled coupling to $B$, $\epsilon_{B}$, of DPDM.
While most of the limits in fig. 3 are quite similar for baryon number, there
are additional limits from the anomalous nature of this gauge group [15] shown
as a solid black line.

The limits shown as solid lines are rather robust and include conservative
estimates on several different levels. Now we will take a more optimistic
point of view and focus especially on the future LISA mission. As noted
earlier in eqs. (8) & (9), longer observation times up to one coherence time
are extremely efficient to enhance the limits. With the peak sensitivity of
LPF lying at around $10^{-3}\,$Hz it would be ideal to have data for around 30
years which of course is far beyond the actual lifetime of the mission.
Nevertheless, the LISA mission may take data for up to 10 years [59] which
means that it naively maximizes the efficiency for frequencies around $\sim
3\cdot 10^{-3}\,$Hz. Together with a general decrease in the noise this will
allow for probing the ultra-low frequency parameter space complementary to the
previous LISA forecasts.111111Optimistically, we will assume a factor 10
improvement from the LPF sensitivity in 2016 and mitigation of the star
tracker noise for our projections.

Previous projections using the planned arm length of around $2.5\cdot
10^{6}\,$km significantly cut into unexplored parameter space as shown in
figs. 3 & 4. These limits are based on looking for TM-TM displacements using
the light-traveling time method and they are strongest around masses of
$10^{-16}\,$eV. However, decreasing the mass by just one order of magnitude
already introduces a decline of the limits by a factor of at least 100. In
contrast to that, our method is well-suited for the lowest frequencies
available because eq. (8) does not depend on the arm length at all. Thus,
there is no suppression of the constraints for low frequencies, i.e. large
coherence lengths, except for the intrinsic sensitivity loss of the
instrument. Indeed, the enhanced reach of our limits at small masses agrees
very well with the findings of Ref. [50] using the KAGRA auxiliary channels.
Even though these auxiliary channels are at best as sensitive as the main
interferometer they clearly outperform the conventional limits in the low mass
region. In conclusion, the LISA mission will provide us with a very powerful
tool to constrain the interactions of DPDM when combining the main channel
analysis with the auxiliary channel analysis. These limits, spanning several
orders of magnitude in mass, will reach deeply into unprobed parameter space.

As mentioned in sec. 2 our approach is not limited to $B-L$ and $B$. In fact,
there is a plethora of additional gauge groups $g$ that will have very similar
limits. These limits just require a proper rescaling procedure depending on
the type of coupling. As noted earlier, there are essentially two types of
couplings in our problem when it comes to analyzing observations involving
different elements. The first one (“$B-L$-like”) is essentially sensitive to
different neutron-to-proton ratios of the different elements while the second
one (“$B$-like”) relies on the smaller differences in binding energies for
different nuclei. Limits that instead depend on the total charge-to-mass ratio
do not suffer from this “binding energy suppression” as can be seen from the
small changes between the LISA projections and the previous LPF limits in fig.
3 to fig. 4. For arbitrary gauge groups with a given combination of baryon and
lepton number $\alpha B-\beta L$ we find that the charge-to-mass ratio can
change from element to element.121212For simplicity, we take $L=L_{e}$ here as
$L_{\mu}$ and $L_{\tau}$ will give no contribution. Therefore, ignoring
different isotopes and changing to nuclear physics notation we find for an
element with atomic number $Z=L$ and mass number $A=B$

$\displaystyle\left(\frac{q}{M}\right)_{\alpha B-\beta L}\simeq$
$\displaystyle\frac{\alpha A-\beta
Z}{Am_{p}}=\alpha\frac{1}{m_{p}}-\beta\frac{Z/A}{m_{p}}\;,$ (12)

where $m_{p}$ denotes the proton mass.

If instead we are interested in the difference between two elements we find

$\displaystyle\Delta\left(\frac{q}{M}\right)_{\alpha B-\beta
L}=\begin{cases}\alpha\Delta\left(\frac{q}{M}\right)_{B}\ \ ,&\beta=0\\\
\beta\Delta\left(\frac{q}{M}\right)_{B-L}\ \ ,&\text{else}\end{cases}\;,$ (13)

using our results for the total charge-to-mass ratios from before. This makes
the distinction between “$B-L$-like” and “$B$-like” immediately clear. Only a
gauge group without coupling to electron number will suffer from the binding
energy suppression. The interesting observation is that the calculation for
$\beta\neq 0$ is already enough to rescale our limits to all possible gauge
groups fulfilling this criterion. We present a selection of groups in table 2.
The second column shows the rescaling for the relative charge-to-mass ratio
and the third column shows the rescaling for the total charge-to-mass ratio.
We note that the decoherence/light-traveling limit rescalings are technically
only valid for Au and the $\beta=0$ are only valid for Au-C systems.
Nevertheless, the rescalings for these cases will still give solid
approximations for the true rescaling factor to arbitrary elements.

Limits | Type of coupling | Relative Acceleration | Decoherence/Light-traveling time (Au)
---|---|---|---
$L_{e}-L_{\mu}$ | $B-L$ | 1 | $79/(197-79)\approx 0.67$
$L_{e}-L_{\tau}$ | $B-L$ | 1 | $79/(197-79)\approx 0.67$
$B-3L_{e}$ | $B-L$ | 3 | $(3\cdot 79-197)/(197-79)\approx 0.34$
$B-3L_{\mu}$ | $B$ | 1 | 1
$B-3L_{\tau}$ | $B$ | 1 | 1
Table 2: Recipe to rescale the limits for more gauge groups. We give the type
of coupling and the corresponding rescaling factor.

In sec. 2 we neglected any contribution from kinetic mixing. Let us briefly
discuss the main reasons why this is well justified. First of all, in-medium
effects lead to an effective suppression of the kinetic mixing if the plasma
mass $\omega_{p}=\sqrt{4\pi\alpha n/m_{e}}$ is larger than the DP mass [36],
i.e. $\epsilon_{\text{KM,eff}}\propto m^{2}/\omega_{p}^{2}$. $m_{e}$ denotes
the electron mass and $n\sim 5e^{-}/$cm3 denotes the electron density in the
interplanetary medium close to Earth [70] implying a plasma mass of $\sim
10^{-10}\,$eV which is much larger than our mass range of interest. Secondly,
both the TMs and the SC are essentially electrically neutral [71]. Finally,
the SC acts like a Faraday cage for the TMs [72]. Of course, the plasma will
also interact with gauged DPs but if we consider the plasma mass of the DPs
due to their direct coupling to SM particles
$\omega_{p,g}\sim\epsilon_{g}\omega_{p}$ we see that the effects are very
small.

Finally, let us put the constraints derived in this letter into larger context
for gauged $B-L$ using the excellent collection of limits from Ref. [73] shown
in fig. 5. Note that the y-axis shows the gauge coupling
$g_{B-L}=\epsilon_{B-L}e$. In addition to the limits shown above, one can also
consider equivalence principle violation searches as direct detection
experiments in a similar mass range [47]. These limits are quite similar to
our work as they also search for a monochromatic DPDM signal on a
“$B-L$-dipole” test mass. Several additional projections are shown in this
plot coming from asteroids [69], atomic interferometry [74], space-based
quantum sensors [75], and future torsion balance experiments [46]. We see that
neither LPF nor LISA is expected to have the best sensitivity in the long run
but as LPF already has available data, this makes it the leading limit over
almost two orders of magnitude in mass and at peak sensitivity it outperforms
the current limits by more than two orders of magnitude in the gauge coupling.
Furthermore, we have outlined why and how a detailed analysis of the LPF data
can push the sensitivity providing excellent motivation for further work.

Figure 5: Limits on the gauge coupling of $B-L$ DPDM. These now include
additional projections and limits as presented in [73].

##  5 Conclusion

In this work we demonstrated how to improve existing DP limits based on the
LPF data. The novel idea is that there is the option to use auxiliary
measurements for the acceleration between the SC and the TMs to constrain the
coupling strength of gauged $B-L$ and $B$ DPDM in analogy to the use of
auxiliary arms in KAGRA [56]. The main advantage is the different atomic
compositions of the test masses and the space craft leading to a relative
acceleration. Relying only on the measurement between the two TMs will lead to
extremely suppressed limits as the TMs react identically to the DPDM field.
The existing literature focused on decoherence and light-traveling time
effects which weakly break this degeneracy at the cost of a massive
suppression for low frequencies where the arm length is much smaller than the
wavelength. While our new limits are also moderately suppressed by the similar
charge-to-mass ratios they are free from any arm length suppression and can
therefore rely on the auxiliary channels working at almost full sensitivity.
For LPF the auxiliary channels at their peak frequency are not significantly
more noisy than the main ($x_{12}$) channel which is an important advantage
for our work. Furthermore, we can cover all three spatial dimension with the
auxiliary channels which will prevent a potential blindness towards specific
DP polarizations.

We showed that even conservative estimates of the LPF results are already able
to probe much new parameter space in the $B-L$ case and at least a small
region for $B$ considering masses around $5\cdot 10^{-18}\,$eV. Our approach
offers an enhancement of up to about 12 orders of magnitude over the most
naive analysis of $B-L$. It is therefore likely that a detailed analysis of
the whole data set of LPF will set even better and thus world-leading limits
over a considerable mass range. Additionally, this work motivates a rigorous
analysis of the reach of LISA using auxiliary channels as our approach might
be highly complementary to the previous forecasts.

## Acknowledgements

We would like to thank Andreas Ringwald for useful discussions. JF is grateful
to Fermilab for its hospitality. This project has received funding from the
European Union’s Horizon Europe research and innovation programme under the
Marie Skłodowska-Curie Staff Exchange grant agreement No 101086085 – ASYMMETRY
and the ITN HIDDeN grant agreement No 860881 – HIDDeN, and from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306 and through the
Emmy Noether Grant No. KA 4662/1-2.

## References

  * [1] G. Bertone, D. Hooper, History of dark matter, Rev. Mod. Phys. 90 (4) (2018) 045002.  arXiv:1605.04909, doi:10.1103/RevModPhys.90.045002.
  * [2] T. Lin, Dark matter models and direct detection, PoS 333 (2019) 009.  arXiv:1904.07915, doi:10.22323/1.333.0009.
  * [3] E. Aprile, et al., First Dark Matter Search with Nuclear Recoils from the XENONnT Experiment, Phys. Rev. Lett. 131 (4) (2023) 041003.  arXiv:2303.14729, doi:10.1103/PhysRevLett.131.041003.
  * [4] J. Aalbers, et al., First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment, Phys. Rev. Lett. 131 (4) (2023) 041002.  arXiv:2207.03764, doi:10.1103/PhysRevLett.131.041002.
  * [5] J. Liu, The first results of PandaX-4T, Int. J. Mod. Phys. D 31 (05) (2022) 2230007.  doi:10.1142/S0218271822300075.
  * [6] J. Jaeckel, A. Ringwald, The Low-Energy Frontier of Particle Physics, Ann. Rev. Nucl. Part. Sci. 60 (2010) 405–437.  arXiv:1002.0329, doi:10.1146/annurev.nucl.012809.104433.
  * [7] A. Ringwald, Ultralight Particle Dark Matter, in: 25th Rencontres de Blois on Particle Physics and Cosmology, 2013.  arXiv:1310.1256.
  * [8] B. Döbrich, Looking for dark matter on the light side, J. Phys. Conf. Ser. 632 (1) (2015) 012004.  arXiv:1501.03274, doi:10.1088/1742-6596/632/1/012004.
  * [9] R. D. Peccei, H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38 (1977) 1440–1443.  doi:10.1103/PhysRevLett.38.1440.
  * [10] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40 (1978) 223–226.  doi:10.1103/PhysRevLett.40.223.
  * [11] F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40 (1978) 279–282.  doi:10.1103/PhysRevLett.40.279.
  * [12] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643 (2016) 1–79.  arXiv:1510.07633, doi:10.1016/j.physrep.2016.06.005.
  * [13] E. Masso, Axions and axion like particles, Nucl. Phys. B Proc. Suppl. 114 (2003) 67–73.  arXiv:hep-ph/0209132, doi:10.1016/S0920-5632(02)01893-5.
  * [14] J. Jaeckel, A force beyond the Standard Model - Status of the quest for hidden photons, Frascati Phys. Ser. 56 (2012) 172–192.  arXiv:1303.1821.
  * [15] J. A. Dror, R. Lasenby, M. Pospelov, New constraints on light vectors coupled to anomalous currents, Phys. Rev. Lett. 119 (14) (2017) 141803.  arXiv:1705.06726, doi:10.1103/PhysRevLett.119.141803.
  * [16] M. Bauer, P. Foldenauer, J. Jaeckel, Hunting All the Hidden Photons, JHEP 07 (2018) 094.  arXiv:1803.05466, doi:10.1007/JHEP07(2018)094.
  * [17] M. Fabbrichesi, E. Gabrielli, G. Lanfranchi, The Dark Photon (5 2020).  arXiv:2005.01515, doi:10.1007/978-3-030-62519-1.
  * [18] E. C. G. Stueckelberg, Die wechselwirkungskräfte in der elektrodynamik und in der feldtheorie der kernkräfte. teil ii und iii, Helvetica Physica Acta 11 (IV) (1938) 299.  doi:10.5169/seals-110855.
URL https://doi.org/10.5169/seals-110855

  * [19] H. Ruegg, M. Ruiz-Altaba, The Stueckelberg field, Int. J. Mod. Phys. A 19 (2004) 3265–3348.  arXiv:hep-th/0304245, doi:10.1142/S0217751X04019755.
  * [20] L. Hui, Wave Dark Matter, Ann. Rev. Astron. Astrophys. 59 (2021) 247–289.  arXiv:2101.11735, doi:10.1146/annurev-astro-120920-010024.
  * [21] J. Redondo, M. Postma, Massive hidden photons as lukewarm dark matter, JCAP 02 (2009) 005.  arXiv:0811.0326, doi:10.1088/1475-7516/2009/02/005.
  * [22] A. E. Nelson, J. Scholtz, Dark Light, Dark Matter and the Misalignment Mechanism, Phys. Rev. D 84 (2011) 103501.  arXiv:1105.2812, doi:10.1103/PhysRevD.84.103501.
  * [23] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, A. Ringwald, WISPy Cold Dark Matter, JCAP 06 (2012) 013.  arXiv:1201.5902, doi:10.1088/1475-7516/2012/06/013.
  * [24] P. W. Graham, J. Mardon, S. Rajendran, Vector Dark Matter from Inflationary Fluctuations, Phys. Rev. D 93 (10) (2016) 103520.  arXiv:1504.02102, doi:10.1103/PhysRevD.93.103520.
  * [25] M. Bastero-Gil, J. Santiago, L. Ubaldi, R. Vega-Morales, Vector dark matter production at the end of inflation, JCAP 04 (2019) 015.  arXiv:1810.07208, doi:10.1088/1475-7516/2019/04/015.
  * [26] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, F. Takahashi, Relic Abundance of Dark Photon Dark Matter, Phys. Lett. B 801 (2020) 135136.  arXiv:1810.07188, doi:10.1016/j.physletb.2019.135136.
  * [27] R. T. Co, A. Pierce, Z. Zhang, Y. Zhao, Dark Photon Dark Matter Produced by Axion Oscillations, Phys. Rev. D 99 (7) (2019) 075002.  arXiv:1810.07196, doi:10.1103/PhysRevD.99.075002.
  * [28] J. A. Dror, K. Harigaya, V. Narayan, Parametric Resonance Production of Ultralight Vector Dark Matter, Phys. Rev. D 99 (3) (2019) 035036.  arXiv:1810.07195, doi:10.1103/PhysRevD.99.035036.
  * [29] A. J. Long, L.-T. Wang, Dark Photon Dark Matter from a Network of Cosmic Strings, Phys. Rev. D 99 (6) (2019) 063529.  arXiv:1901.03312, doi:10.1103/PhysRevD.99.063529.
  * [30] K. Nakayama, Vector Coherent Oscillation Dark Matter, JCAP 10 (2019) 019.  arXiv:1907.06243, doi:10.1088/1475-7516/2019/10/019.
  * [31] Y. Nakai, R. Namba, Z. Wang, Light Dark Photon Dark Matter from Inflation, JHEP 12 (2020) 170.  arXiv:2004.10743, doi:10.1007/JHEP12(2020)170.
  * [32] A. Ahmed, B. Grzadkowski, A. Socha, Gravitational production of vector dark matter, JHEP 08 (2020) 059.  arXiv:2005.01766, doi:10.1007/JHEP08(2020)059.
  * [33] E. W. Kolb, A. J. Long, Completely dark photons from gravitational particle production during the inflationary era, JHEP 03 (2021) 283.  arXiv:2009.03828, doi:10.1007/JHEP03(2021)283.
  * [34] B. Salehian, M. A. Gorji, H. Firouzjahi, S. Mukohyama, Vector dark matter production from inflation with symmetry breaking, Phys. Rev. D 103 (6) (2021) 063526.  arXiv:2010.04491, doi:10.1103/PhysRevD.103.063526.
  * [35] H. Firouzjahi, M. A. Gorji, S. Mukohyama, B. Salehian, Dark photon dark matter from charged inflaton, JHEP 06 (2021) 050.  arXiv:2011.06324, doi:10.1007/JHEP06(2021)050.
  * [36] J. Redondo, Helioscope Bounds on Hidden Sector Photons, JCAP 07 (2008) 008.  arXiv:0801.1527, doi:10.1088/1475-7516/2008/07/008.
  * [37] H. An, M. Pospelov, J. Pradler, New stellar constraints on dark photons, Phys. Lett. B 725 (2013) 190–195.  arXiv:1302.3884, doi:10.1016/j.physletb.2013.07.008.
  * [38] E. Hardy, R. Lasenby, Stellar cooling bounds on new light particles: plasma mixing effects, JHEP 02 (2017) 033.  arXiv:1611.05852, doi:10.1007/JHEP02(2017)033.
  * [39] S.-P. Li, X.-J. Xu, Production rates of dark photons and $Z^{\prime}$ in the Sun and stellar cooling bounds (4 2023).  arXiv:2304.12907.
  * [40] J. Frerick, F. Kahlhoefer, K. Schmidt-Hoberg, A’ view of the sunrise: boosting helioscopes with angular information, JCAP 03 (2023) 001.  arXiv:2211.00022, doi:10.1088/1475-7516/2023/03/001.
  * [41] T. A. Wagner, S. Schlamminger, J. H. Gundlach, E. G. Adelberger, Torsion-balance tests of the weak equivalence principle, Class. Quant. Grav. 29 (2012) 184002.  arXiv:1207.2442, doi:10.1088/0264-9381/29/18/184002.
  * [42] P. Touboul, et al., MICROSCOPE Mission: Final Results of the Test of the Equivalence Principle, Phys. Rev. Lett. 129 (12) (2022) 121102.  arXiv:2209.15487, doi:10.1103/PhysRevLett.129.121102.
  * [43] T. Kumar Poddar, S. Mohanty, S. Jana, Constraints on long range force from perihelion precession of planets in a gauged $L_{e}-L_{\mu,\tau}$ scenario, Eur. Phys. J. C 81 (4) (2021) 286.  arXiv:2002.02935, doi:10.1140/epjc/s10052-021-09078-9.
  * [44] Y.-D. Tsai, Y. Wu, S. Vagnozzi, L. Visinelli, Novel constraints on fifth forces and ultralight dark sector with asteroidal data, JCAP 04 (2023) 031.  arXiv:2107.04038, doi:10.1088/1475-7516/2023/04/031.
  * [45] Y.-D. Tsai, D. Farnocchia, M. Micheli, S. Vagnozzi, L. Visinelli, Constraints on fifth forces and ultralight dark matter from OSIRIS-REx target asteroid Bennu (9 2023).  arXiv:2309.13106.
  * [46] P. W. Graham, D. E. Kaplan, J. Mardon, S. Rajendran, W. A. Terrano, Dark Matter Direct Detection with Accelerometers, Phys. Rev. D 93 (7) (2016) 075029.  arXiv:1512.06165, doi:10.1103/PhysRevD.93.075029.
  * [47] E. A. Shaw, M. P. Ross, C. A. Hagedorn, E. G. Adelberger, J. H. Gundlach, Torsion-balance search for ultralow-mass bosonic dark matter, Phys. Rev. D 105 (4) (2022) 042007.  arXiv:2109.08822, doi:10.1103/PhysRevD.105.042007.
  * [48] A. Pierce, K. Riles, Y. Zhao, Searching for Dark Photon Dark Matter with Gravitational Wave Detectors, Phys. Rev. Lett. 121 (6) (2018) 061102.  arXiv:1801.10161, doi:10.1103/PhysRevLett.121.061102.
  * [49] A. L. Miller, L. Mendes, First search for ultralight dark matter with a space-based gravitational-wave antenna: LISA Pathfinder, Phys. Rev. D 107 (6) (2023) 063015.  arXiv:2301.08736, doi:10.1103/PhysRevD.107.063015.
  * [50] Y. Michimura, T. Fujita, S. Morisaki, H. Nakatsuka, I. Obata, Ultralight vector dark matter search with auxiliary length channels of gravitational wave detectors, Phys. Rev. D 102 (10) (2020) 102001.  arXiv:2008.02482, doi:10.1103/PhysRevD.102.102001.
  * [51] B. Holdom, Two U(1)’s and Epsilon Charge Shifts, Phys. Lett. B 166 (1986) 196–198.  doi:10.1016/0370-2693(86)91377-8.
  * [52] K. S. Babu, C. F. Kolda, J. March-Russell, Implications of generalized Z - Z-prime mixing, Phys. Rev. D 57 (1998) 6788–6792.  arXiv:hep-ph/9710441, doi:10.1103/PhysRevD.57.6788.
  * [53] N. W. Evans, C. A. J. O’Hare, C. McCabe, Refinement of the standard halo model for dark matter searches in light of the Gaia Sausage, Phys. Rev. D 99 (2) (2019) 023012.  arXiv:1810.11468, doi:10.1103/PhysRevD.99.023012.
  * [54] A. Caputo, A. J. Millar, C. A. J. O’Hare, E. Vitagliano, Dark photon limits: A handbook, Phys. Rev. D 104 (9) (2021) 095029.  arXiv:2105.04565, doi:10.1103/PhysRevD.104.095029.
  * [55] L. Giulicchi, S.-F. Wu, T. Fenal, Attitude and orbit control systems for the lisa pathfinder mission, Aerospace Science and Technology 24 (1) (2013) 283–294, vFE-2.  doi:https://doi.org/10.1016/j.ast.2011.12.002.
URL https://www.sciencedirect.com/science/article/pii/S127096381100201X

  * [56] S. Morisaki, T. Fujita, Y. Michimura, H. Nakatsuka, I. Obata, Improved sensitivity of interferometric gravitational wave detectors to ultralight vector dark matter from the finite light-traveling time, Phys. Rev. D 103 (5) (2021) L051702.  arXiv:2011.03589, doi:10.1103/PhysRevD.103.L051702.
  * [57] D. Budker, P. W. Graham, M. Ledbetter, S. Rajendran, A. Sushkov, Proposal for a Cosmic Axion Spin Precession Experiment (CASPEr), Phys. Rev. X 4 (2) (2014) 021030.  arXiv:1306.6089, doi:10.1103/PhysRevX.4.021030.
  * [58] L. P. collaboration, LISA Pathfinder: First steps to observing gravitational waves from space, J. Phys. Conf. Ser. 840 (1) (2017) 012001.  doi:10.1088/1742-6596/840/1/012001.
  * [59] P. Amaro-Seoane, et al., Laser Interferometer Space Antenna (2 2017).  arXiv:1702.00786.
  * [60] M. Armano, et al., Sub-Femto- g Free Fall for Space-Based Gravitational Wave Observatories: LISA Pathfinder Results, Phys. Rev. Lett. 116 (23) (2016) 231101.  doi:10.1103/PhysRevLett.116.231101.
  * [61] M. Armano, et al., Beyond the Required LISA Free-Fall Performance: New LISA Pathfinder Results down to 20 $\mu$Hz, Phys. Rev. Lett. 120 (6) (2018) 061101.  doi:10.1103/PhysRevLett.120.061101.
  * [62] M. Armano, et al., Measuring the Galactic Cosmic Ray Flux with the LISA Pathfinder Radiation Monitor, Astropart. Phys. 98 (2018) 28–37.  arXiv:1711.07427, doi:10.1016/j.astropartphys.2018.01.006.
  * [63] H. E. Audley, LISA Pathfinder: Optical Metrology System monitoring during operations, J. Phys. Conf. Ser. 840 (1) (2017) 012034.  doi:10.1088/1742-6596/840/1/012034.
  * [64] M. Armano, et al., Capacitive sensing of test mass motion with nanometer precision over millimeter-wide sensing gaps for space-borne gravitational reference sensors, Phys. Rev. D 96 (6) (2017) 062004.  doi:10.1103/PhysRevD.96.062004.
  * [65] G. D. Racca, P. W. McNamara, The LISA Pathfinder mission: Tracing Einstein’s geodesics in space, Space Sci. Rev. 151 (2010) 159–181.  doi:10.1007/s11214-009-9602-x.
  * [66] M. Armano, et al., LISA Pathfinder Platform Stability and Drag-free Performance, Phys. Rev. D 99 (8) (2019) 082001.  arXiv:1812.05491, doi:10.1103/PhysRevD.99.082001.
  * [67] https://www.nist.gov/pml/atomic-weights-and-isotopic-compositions-relative-atomic-masses.
  * [68] P. F. de Salas, K. Malhan, K. Freese, K. Hattori, M. Valluri, On the estimation of the Local Dark Matter Density using the rotation curve of the Milky Way, JCAP 10 (2019) 037.  arXiv:1906.06133, doi:10.1088/1475-7516/2019/10/037.
  * [69] M. A. Fedderke, A. Mathur, Asteroids for ultralight dark-photon dark-matter detection, Phys. Rev. D 107 (4) (2023) 043004.  arXiv:2210.09324, doi:10.1103/PhysRevD.107.043004.
  * [70] K. Issautier, C. Perche, S. Hoang, C. Lacombe, M. Maksimovic, J.-L. Bougeret, C. Salem, Solar wind electron density and temperature over solar cycle 23: Thermal noise measurements on wind, Advances in Space Research 35 (12) (2005) 2141–2146, the Dynamic Heliosphere.  doi:https://doi.org/10.1016/j.asr.2005.04.085.
URL https://www.sciencedirect.com/science/article/pii/S0273117705006022

  * [71] M. Armano, et al., Charging of free-falling test masses in orbit due to cosmic rays: Results from LISA Pathfinder, Phys. Rev. D 107 (6) (2023) 062007.  arXiv:2211.09309, doi:10.1103/PhysRevD.107.062007.
  * [72] S. Chaudhuri, P. W. Graham, K. Irwin, J. Mardon, S. Rajendran, Y. Zhao, Radio for hidden-photon dark matter detection, Phys. Rev. D 92 (7) (2015) 075012.  arXiv:1411.7382, doi:10.1103/PhysRevD.92.075012.
  * [73] C. O’Hare, cajohare/axionlimits: Axionlimits, https://cajohare.github.io/AxionLimits/ (Jul. 2020).  doi:10.5281/zenodo.3932430.
  * [74] M. Abe, et al., Matter-wave Atomic Gradiometer Interferometric Sensor (MAGIS-100), Quantum Sci. Technol. 6 (4) (2021) 044003.  arXiv:2104.02835, doi:10.1088/2058-9565/abf719.
  * [75] P. Wolf, et al., STE-QUEST: Space Time Explorer and QUantum Equivalence principle Space Test (11 2022).  arXiv:2211.15412.


# Reinterpretation of the Fermi acceleration of cosmic rays in terms of the
ballistic surfing acceleration in supernova shocks

Krzysztof Stasiewicz

###### Abstract

The applicability of first-order Fermi acceleration in explaining the cosmic
ray spectrum has been reexamined using recent results on shock acceleration
mechanisms from the Multiscale Magnetospheric mission in Earth’s bow shock. It
is demonstrated that the Fermi mechanism is a crude approximation of the
ballistic surfing acceleration (BSA) mechanism. While both mechanisms yield
similar expressions for the energy gain of a particle after encountering a
shock once, leading to similar power-law distributions of the cosmic ray
energy spectrum, the Fermi mechanism is found to be inconsistent with
fundamental equations of electrodynamics. It is shown that the spectral index
of cosmic rays is determined by the average magnetic field compression rather
than the density compression, as in the Fermi model. It is shown that the knee
observed in the spectrum at an energy of $5\times 10^{15}$ eV could correspond
to ions with a gyroradius comparable to the size of shocks in supernova
remnants. The BSA mechanism can accurately reproduce the observed spectral
index $s=-2.5$ below the knee energy, as well as a steeper spectrum, $s=-3$,
above the knee. The acceleration time up to the knee, as implied by BSA, is on
the order of 300 years. First-order Fermi acceleration does not represent a
physically valid mechanism and should be replaced by ballistic surfing
acceleration in applications or models related to quasi-perpendicular shocks
in space. It is noted that BSA, which operates outside of shocks, was
previously misattributed to shock drift acceleration (SDA), which operates
within shocks.

###### keywords:

Cosmic rays , Shock waves , Acceleration of particles , Supernova remnants

††journal:

[label1]organization=Space Research Centre, Polish Academy of
Sciences,addressline=Bartycka 18A, city=Warszawa, postcode=00-716,
country=Poland label2label2footnotetext: Corresponding author:
<EMAIL_ADDRESS>
##  1 Introduction

Cosmic rays exhibit a power-law energy distribution with a spectral index of
$s\approx-2.5$ below the knee at an energy of $K\approx 5{\times}10^{15}$ eV,
and a steeper slope, $s\approx-3$, above the knee and below the ankle at
$5{\times}10^{18}$ eV (Hillas, 1984; Helder et al., 2012). Understanding the
spectrum and acceleration mechanisms operating at such high energies has been
an important, yet not fully explained problem in astrophysics. Fermi (1949)
proposed that the cosmic ray spectrum could correspond to ions accelerated by
bouncing off magnetic clouds in the interstellar medium. It is now believed
that these clouds represent magnetic turbulence responsible for second-order
Fermi acceleration. On the other hand, first-order Fermi acceleration is
generally regarded as the major mechanism responsible for the formation of
cosmic rays via the diffusive shock acceleration (DSA) process (Bell, 1978;
Blandford and Eichler, 1987; Longair, 2011).

The recent Magnetospheric Multiscale (MMS) mission (Burch et al., 2016),
consisting of four satellites flying through the Earth’s bow shock, has
provided the best measurements of collisionless shocks to date. The MMS
satellites fly occasionally with separation distances of 20 km, which is
smaller than the thermal proton gyroradius of 100 km. In this mission the
electric, $\bm{E}$ and magnetic, $\bm{B}$ fields are sampled at the rate of
8192 s-1, while particles distribution functions are measured with time
resolution of 30 ms for electrons and 150 ms for ions.

Based on MMS data, it has been established that thermalization, heating, and
acceleration of ions and electrons in collisionless shocks are related to four
plasma processes: SWE (stochastic wave energization), TTT (transit time
thermalization), BSA (ballistic surfing acceleration), and QAH (quasi-
adiabatic heating), explained further in the text (see also list of acronyms).
Understanding processes in the bow shock can help understand more powerful
astrophysical shocks, such as those created during supernova explosions, which
are most likely involved in the acceleration of cosmic rays.

The processes SWE and TTT are stochastic in nature, related to deterministic
chaos. They rely on strong gradients of the electric and magnetic fields that
lead to randomization of particle orbits and efficient stochastic heating on
the timescale of one gyroperiod (Stasiewicz, 2023b). TTT thermalizes streaming
ions on magnetic field gradients within a fraction of the gyroperiod and does
not require any waves, instabilities, or anomalous collisions. SWE works on
electric field gradients and can accelerate protons to a few hundred keV
(Stasiewicz and Eliasson, 2021; Stasiewicz et al., 2021; Stasiewicz and Kłos,
2022). SWE of electrons is responsible both for heating and formation of flat-
top electron distributions (Stasiewicz, 2023a). Quasi-adiabatic heating (QAH)
operates on particles with a conserved first adiabatic invariant
($v_{\perp}^{2}/B=\mathrm{constant}$), which requires the gyroradius $r_{c}$
to be much smaller than the width of the shock, $D\equiv
B|\bm{\nabla}B|^{-1}$.

Ballistic surfing acceleration (BSA) operates on particles with a gyroradius
$r_{c}\gg D$, which applies to superthermal ions and mildly relativistic
electrons (Stasiewicz, 2023b). It bears resemblance to the shock drift
acceleration (SDA) previously discussed by various authors (Jokipii, 1982;
Jones and Ellison, 1991; Zank et al., 1996). Another similar process is shock
surfing acceleration (SSA), where low energy particles drift along the shock
front due to surface waves (Shapiro et al., 2001; Lever et al., 2001; Hoshino
and Shimada, 2002). In all these processes, energization is due to particle
motion along the convection electric field. However, the term SDA implies
$\bm{\nabla}B$ drift acceleration, which occurs within the ramp. In contrast,
BSA occurs outside the ramp, where particles with large gyroradii engage in
ballistic surfing and do not experience $\bm{\nabla}B$ or wave effects.

As is well known, plasma can also be effectively heated by waves at resonance
frequencies: cyclotron, lower hybrid, upper hybrid, and plasma frequencies.
Resonant heating is commonly used in laboratory plasmas, but does not appear
to be important in shocks.

A question arises, what is the relation of the Fermi mechanism and the DSA
process with the above-described processes identified at the bow shock. In the
next section we shall study heating of energetic ions in a model shock,
searching for signatures of Fermi acceleration. We find that the first order
Fermi acceleration represents a crude approximation of BSA. We show that the
cosmic ray spectrum and the presence of the knee can be explained exclusively
by BSA.

##  2 Ballistic surfing acceleration

To understand acceleration of high energy particles we consider motion of test
particles with velocity $v$ much larger than the thermal velocity $v_{Ti}$ of
ions that maintain the shock. Trajectories of particles with rest mass $m_{0}$
and charge $q$ are described by the momentum equation

$\frac{\mathrm{d}\bm{p}}{\mathrm{d}t}=q(\bm{E}+\bm{v}\times\bm{B}),$ (1)

where $\bm{p}=\gamma m_{0}\bm{v},\;\gamma=(1-v^{2}/c^{2})^{-1/2}$, and $c$ is
the speed of light. Here, we utilize the shock reference frame in a geometry
where the $\hat{\bm{x}}$ axis is in the negative direction to the shock
normal, $\hat{\bm{y}}$ is in the direction of the convection electric field,
and the magnetic field is in the $[B_{x0},0,B_{z}(x)]$ plane. The kinetic
energy increase of a particle implied by this equation is

$\Delta K=q\int\bm{E}\cdot\bm{v}\,dt.$ (2)

There are three types of electric fields that determine particle dynamics and
energization processes in the shock frame: $\hat{\bm{y}}E_{y}$ – the
convection field, constant across one-dimensional structures,
$\hat{\bm{x}}E_{S}(x)$ – the cross-shock electric field, maintained by the
electron pressure gradient, and $\tilde{\bm{E}}(\bm{r},t)$ – the wave electric
field. In the reference frame moving with the convection velocity, $E_{y}$
vanishes, but acceleration is facilitated by the inductive electric field,
$\bm{\nabla}\times\bm{E}=-\partial\bm{B}/\partial t$ observed by a particle
moving in a time varying magnetic field. SWE works best on thermal particles
and is ineffective for high-speed particles, so we exclude electrostatic waves
from consideration in this paper.

The model shock is described in Stasiewicz (2023b) and represents a magnetic
ramp with compression $c_{B}=B_{d}/B_{u}$ between downstream and upstream
values. A normal component $B_{x0}$ implies the upstream field angle
$\cos\eta=B_{x0}/B_{u}$ from the normal direction. Magnetic turbulence is not
included in the shock model, but it can be easily implemented similarly to
electrostatic waves. It is expected to lead to the isotropisation of particle
distributions and could also scatter particles back to the shock.

Parameters that enter the model are as follows: the upstream sonic Mach
number, $M=V_{u}/v_{Ti}=8$, the ratio of the thermal ion gyroradius $r_{ci}$
to the width of the shock ramp, $r_{ci}/D=1$, and compression $c_{B}=4$. The
transit time thermalization parameter is set to $\chi_{B}=8$ and the
stochastic parameter for the cross-shock electric field to $\chi_{S}=1.5$; see
Stasiewicz (2023b) for definitions.

Two ions are injected into an oblique shock $\eta=85^{\circ}$ at $x=-40$ in
units of the shock width and followed with differential equations described in
Stasiewicz (2023b). Fig. 1(a) shows ion trajectories in plane $(x,y)$. The
initial velocity consists of the $E{\times}B$ drift normalized by the upstream
thermal speed, $\bm{u}_{E{\times}B}=\bm{V}_{E{\times}B}/v_{Ti}$, with
additional velocities: $u_{x}=+20$ (blue), and $u_{y}=-20$ (red). Panel (b)
shows the total kinetic energies $u^{2}$ of ions along the respective
trajectories, panel (c) shows gyration energies only
$(\bm{u}_{\perp}-\bm{u}_{E{\times}B})^{2}$, and panel (d) shows
$u_{\parallel}^{2}$. Adiabatic projection of the initial perpendicular energy
of the red ion is shown as the black curve $b(x)u_{0\perp}^{2}$, where
$b(x)=B(x)/B_{u}$. Both ions behave non-adiabatically and do not follow the
adiabatic projection.

We observe that the blue ion makes three crossings of the shock at $x=0$,
while the red ion makes seven crossings before being transmitted downstream.
Each crossing is associated with an increase in gyration energy, as seen in
panel (c). The kinetic energy in panel (b) correlates with particle movements
in the $\pm y$ direction, along $E_{y}$, which has been defined as the BSA by
Stasiewicz (2023b). This obvious acceleration process operating outside the
shock was previously misattributed to SDA, which is related to $\bm{\nabla}B$
drift and operates within the shock. However, particles undergoing BSA in the
gradient-free zone are unaware of the existence of the shock, and therefore
cannot experience shock drift acceleration.

Temperature variations in panel (c) and the transfer of energy into the
parallel motion shown in panel (d) occur only during shock transitions at
$x\approx 0$ and are caused by TTT, which is related to the inductive electric
field or to the changing magnetic field direction in the shock frame.

Figure 1: Ion trajectories and heating by TTT and BSA in an oblique shock with
a shock angle of $\eta=85^{\circ},\;\chi_{B}=8$, the cross-shock electric
field $\chi_{S}=1.5$, compression $c_{B}=4$, and thickness ratio $r_{ci}/D=1$,
without waves. (a): Trajectories of two ions injected at $x=-40$ with a sonic
Mach number $M=8$ and additional velocities: $u_{x}=+20$ (blue) and
$u_{y}=-20$ (red) in units of the upstream ion thermal velocity $v_{Ti}$. (b):
Total kinetic energies of ions along the respective trajectories. (c): Thermal
(gyration) energies of ions along trajectories. The black curve shows
theoretical adiabatic heating for the shock profile given by $b(x)\equiv
B(x)/B_{u}$. (d): Parallel energies of ions along trajectories. Typical values
upstream of the bow shock are: $T_{i}\approx 20$ eV, $v_{Ti}\approx 60$ km s-1
, $r_{ci}\approx 100$ km, $B_{u}=5$ nT.  Figure 2: BSA and TTT of ions in the
same format as in Fig. 1 but for a shock angle $\eta=70^{\circ}$. The red ion
is reflected upstream. Ion gyration energy is $u_{\perp}^{2}=400$ and the
drift energy $M^{2}=64$ at the starting position $x=-40$.

By decreasing the shock angle to $\eta=70^{\circ}$ in Fig. 2 we observe the
reflection of the red ion back into the upstream region. Particle reflections
are more common for smaller angles $\eta$, and are not related to the magnetic
mirror force, which does not apply for large gyroradius particles, $r_{c}\gg
D$.

For particles reflected upstream, the energy gain is directed into parallel
motion rather than perpendicular heating. The transfer of energy into parallel
motion, as shown in panels (d), is not related to the cross-shock electric
field, which has a parallel component $E_{\parallel}=E_{S}(x)\cos\eta$, as the
patterns remain the same when $E_{S}$ is set to zero. The cross-shock
potential is $e\Delta\Phi\approx 2(c_{B}-1)=6$ in normalized energy units of
$T_{eu}\sim T_{iu}$, so it has negligible effect even for zero temperature
ions with energy $u^{2}=M^{2}=64$; see Stasiewicz (2023b), not to mention the
high-energy ions, $u^{2}\sim 1000$ in the simulations here. Furthermore, the
cross-shock potential energy does not accumulate, but cancels out after each
gyration around the shock.

One should not confuse cyclotron turning points with shock reflections.
Cyclotron turning points recur periodically, once during every gyroperiod,
while shock reflections are singular events that occur only once. Turning
points can manifest anywhere, contingent upon the initial conditions and the
gyroradius. Sensitivity to initial conditions is a hallmark of deterministic
chaos. Figures 1–2b illustrate that neither the process of reflection nor the
transitions through the shock at $x\approx 0$ are associated with significant
changes in particle energy.

Figure 3: The explanation of ballistic surfing acceleration: Three ions with a
perpendicular velocity of $v_{\perp}=50v_{Ti}$ are injected at incident angles
of $0^{\circ}$ (black), $+10^{\circ}$ (blue), and $-10^{\circ}$ (red) into the
perpendicular shock located at position $x=0$. The parameters are set as
follows: $c_{B}=4$, Mach number $M=3,\;r_{ci}/D=1$, and the cross-shock
electric field $\chi_{S}=1.5$. The energy gain after one gyration across the
shock is $\Delta K\approx 2qE_{y}(r_{cu}-r_{cd})$, which leads directly to Eq.
(4). SDA occurs within the ramp denoted by vertical lines, where $|x|/D<1$,
while BSA operates in the gradient-free zone for $|x|/D\geq 1$.

##  3 BSA and formation of the cosmic ray spectrum

As observed in the previous section, energetic particles are accelerated
outside shocks through the ballistic surfing acceleration process by the
convection electric field $E_{y}$, rather than through reflections as
suggested by Fermi (1949). A necessary condition for BSA is that the
particle’s gyroradius is larger than the shock width. The BSA process is
elucidated in Fig. 3, where it can be observed that ions moving in the $+y$
direction in the upstream part of the orbit increase kinetic energy
($qE_{y}v_{y}>0$), while those moving in the $-y$ direction in the downstream
orbit decrease energy ($qE_{y}v_{y}<0$). Because $B_{d}>B_{u}$, the downstream
gyroradius $r_{cd}$ is smaller than the upstream one $r_{cu}$, and particles
always gain energy. Within the shock ramp denoted by vertical lines, the
particles experience gyrocenter drifts induced by $\bm{\nabla}B$ and are
subject to the combined convection and cross-shock electric fields,
$\hat{\bm{y}}E_{y}+\hat{\bm{x}}E_{S}(x)$. The energy gain after a full
rotation around the shock consists of three terms:

$\Delta K=\int_{|x|\geq
D}^{(\mathrm{BSA})}qE_{y}v_{y}dt+\int_{|x|<D}^{(\mathrm{SDA})}qE_{y}v_{y}dt+\int_{|x|<D}qE_{S}v_{x}dt$
(3)

where the first term corresponds to ballistic surfing acceleration outside the
shock ramp, the second term to shock drift acceleration within the ramp, and
the third term to acceleration by the cross-shock electric field. The first
BSA term is the only significant one for high-energy particles because the
integration domain is much larger than in the other two integrals ($2\pi
r_{c}\gg 2D$). It can be concluded that BSA is the primary mechanism for
accelerating particles with large gyroradii ($r_{c}\gg D$), while the effects
of SDA are negligible due to the very short interaction time compared to the
gyroperiod. Conversely, for small gyroradii particles ($r_{c}\ll D$), SDA
serves as the primary mechanism for the acceleration and BSA does not apply.
However, this process is equivalent to adiabatic heating, so SDA can be
considered a subset of the more general QAH (Stasiewicz and Eliasson, 2023).

After one gyration across the shock the energy gain implied by Eq. (3) is
$\Delta K\approx 2qE_{y}(r_{cu}-r_{cd})$, as illustrated in Fig. 3. For
relativistic particles with kinetic energy $K\approx pc$ and gyroradius
$r_{c}=p_{\perp}/qB$, the energy gain is

$\Delta K\approx
2qE_{y}\frac{1}{\pi}\int_{0}^{\pi}(r_{cu}-r_{cd})\sin\theta\,d\theta=K\frac{(c_{B}-1)}{c_{B}}\frac{V_{u}}{c}$
(4)

where $V_{u}=E_{y}/B_{u}$ is the upstream convection velocity, and the
integral is the average over pitch angles, assuming isotropic distribution.
Averaging over incident angles of particles entering the shock would introduce
a numerical factor $g\lesssim 1$ to this equation; see Fig. 3. Such a factor
could be incorporated into a slightly modified value of $c_{B}$, which is a
free parameter of the model. The effective compression will be then
$c_{B}^{\prime}=c_{B}/[c_{B}-g(c_{B}-1)]$.

Equation (4) is analogous to the expression derived through a wholly distinct
approach by Bell (1978). This formulation is commonly referred to as the
first-order Fermi acceleration:

$\Delta K_{F}=K\frac{(c_{N}-1)}{3}\frac{V_{u}}{c}.$ (5)

This equation was derived from the difference in particle energy between two
inertial frames, characterized by the velocity ratio
$V_{u}/V_{d}=n_{d}/n_{u}\equiv c_{N}$; see Longair (2011) for derivation.
Here, $c_{N}$ denotes the density compression, typically set to the standard
value $c_{N}=4$. It’s important to note that this equation was mistakenly
linked to energization. However, the energy difference between two inertial
frames is simply a scalar value derived from the Lorentz transformation and
doesn’t relate to acceleration, which is defined by Eq. (2), and must be
computed within the same reference frame. This yields the accurate expression
(4), which depends on $c_{B}$ but remains independent of $c_{N}$. Figs. 1-2b
show that particles undergo energy gain over an extended period by ballistic
surfing in the upstream region and lose energy in the downstream region, as
described by $\mathrm{d}K/\mathrm{d}t=qE_{y}v_{y}$. Contrary to suppositions
of the Fermi model, shock transitions, $|x|<D$, are not associated with
significant energy changes, as can be seen in the aforementioned figures and
in Eq. (3).

Cumulative energy increases caused by multiple encounters of particles with
shocks, combined with scattering by turbulence outside the ramp represent
diffusion in velocity space, described as DSA. Because the BSA energy increase
in Eq. (4) is similar to the first-order Fermi acceleration (5), we can follow
the standard approach to DSA, as described, for example, by Bell (1978) and
Longair (2011). Equation (4) implies that the particle energy after one
interaction with a shock is

$K_{1}=hK_{0};\quad h=1+(1-1/c_{B})V_{u}/c.$ (6)

This acceleration is not related to the reflection or to the shock crossing
process but to the full gyration across the shock ramp with $c_{B}>1$. Let $P$
be the probability that particles remain in the shock region after one
interaction or gyration. Then, after $j$ interactions there are $N=N_{0}P^{j}$
particles with energies $K=K_{0}h^{j}$. Eliminating $j$, one obtains
$N/N_{0}=(K/K_{0})^{\ln P/\ln h}$, where $N$ is the number of particles that
reached energy $K$ and can be accelerated further. Using Longair (2011) value
for probability, $P\approx 1-V_{u}/c$, we find the spectrum of cosmic rays
predicted by the BSA model

$\frac{\mathrm{d}N}{\mathrm{d}K}\propto K^{s}$ (7)

where

$s=\frac{\ln P}{\ln h}-1=-\frac{2c_{B}-1}{c_{B}-1}$ (8)

which can be compared with the spectral index predicted by the Fermi/DSA model
(Longair, 2011)

$s_{F}=\frac{\ln P}{\ln h_{F}}-1=-\frac{c_{N}+2}{c_{N}-1}.$ (9)

In deriving this equation, $h_{F}=1+\frac{1}{3}(c_{N}-1)V_{u}/c$ was used from
incorrect Eq. (5). The Fermi/DSA formulation neglects the impact of the
electric field, which is crucial for particle energization. Moreover, it
introduces an erroneous dependence on density compression, inconsistent with
the fundamental equations (1) and (2) which determine heating, but remain
independent of plasma density. Eq. (5) is clearly non-physical, because it
diverges for large values of $c_{N}$, while the numerical factor in the
correct Eq. (4) approaches unity for $c_{B}\gg 1$.

While not being physically valid, Eq. (9) provides one correct value
$s_{F}=-2.5$ for $c_{N}=3$, coinciding with Eq. (8) for $c_{B}=3$, and both
agree with the measured index $s$. This single, incidentally correct value of
the Fermi/DSA model can possibly explain its success and popularity despite
erroneous physical assumptions. The error has not been disclosed earlier,
probably due to the observed relation $n\propto B$ in shocks, which yields
$c_{N}\approx c_{B}$. The physical dependence of particle energy gain on
$c_{B}$ has been confused with a circumstantial dependence on $c_{N}$.

The BSA model would reproduce the observed index $s\approx-2.5$ for the
effective shock compression $c_{B}\approx 3$. For the standard value
$c_{B}=4$, implied from Rankine-Hugoniot relations (Kennel, 1988), we obtain
$s=-2.3$. Smaller compressions $c_{B}=2$ would lead to $s=-3$ observed above
the knee, while magnetic walls, $c_{B}\gg 1$ would lead to the asymptotic
value $s=-2$.

Although BSA would function effectively at arbitrarily high energy, the
acceleration in the upstream region will be counteracted by deceleration in
the downstream region as the diameter of the orbit on the compressed side
approaches the shock length $L$, see Fig. 3. The condition $2r_{cd}\sim L$
would result in a knee in the spectrum located at energy

$K_{L}\sim\frac{qc}{2}\langle LB_{d}\rangle.$ (10)

Interestingly, this does not depend on shock velocity or particle mass. The
observed distribution of supernova remnant sizes, $L_{\mathrm{SNR}}$, ranges
from 1 pc ($=3\times 10^{16}$ m) to 200 pc (Badenes et al., 2010). The
observed knee energy, $K_{L}\approx 5\times 10^{15}$ eV, is derived from Eq.
(10) for $\langle LB_{d}\rangle\sim 1$ nT pc, which could correspond to, for
example, $B_{d}\approx 1$ nT on the inner (compressed) side of a spherically
expanding supernova shock and $L\sim 1$ pc. Shocks with length $L<1$ pc inside
supernova remnants with an effective compression $c_{B}=3$ would lead to
$s=-2.5$, which could explain the energy spectrum below the knee.

However, even particles with a gyroradius much larger than the shock length
can undergo further acceleration. This scenario may occur when downstream
flows become stagnant at some distance from the shock, with $V_{d}\sim 0$,
leading to the vanishing of the electric field in this region, $E_{yd}\sim 0$.
Particles gyrating in the downstream region would not experience deceleration,
as depicted in Figs. 1-2b, but they would still be subject to acceleration in
the upstream region, $\mathrm{d}K/\mathrm{d}t=qE_{y}v_{y}>0$. This process
could yield cosmic ray energies up to the ankle at $5\times 10^{18}$ eV in
shocks shorter than the gyroradius, $L<r_{c}$. The spectral index $s\approx-3$
in this energy range can be achieved with $c_{B}\approx 2$; however, other
factors may also play a significant role in determining the observed index.

##  4 Discussion

BSA requires that the gyroradius of particles is larger than the shock width,
$r_{c}/D>1$, which is fulfilled at the bow shock by protons with energy higher
than 100 eV and by electrons with energy higher than 180 keV. Initial heating
of cold ions in shocks is accomplished by TTT and SWE mechanisms, which can be
seen in Fig. 4 of Stasiewicz (2023b), where streaming protons with temperature
$T_{i}\approx 20$ eV, are TTT thermalized and SWE energized to 400 $T_{i}=8$
keV within a gyroperiod. In quasi-parallel shocks ions are accelerated by the
SWE mechanism to $\sim$ 200 keV which corresponds to the $E{\times}B$ velocity
in the wave electric field, $V_{E{\times}B}=\tilde{E}_{\perp}/B$ with
$\tilde{E}_{\perp}\sim 100$ mV m-1 being much larger than the convection field
$E_{y}\sim 5$ mV m-1 (Stasiewicz and Eliasson, 2021; Stasiewicz et al., 2021).
Further acceleration of these ions could continue by means of BSA in
subsequent shock encounters.

Let us assume that the injection energy for BSA is $K_{0}=10$ keV, and the
shock parameters are $c_{B}=4$, and $V_{u}=10,000$ km s-1 . Using the exact
form of Eq. (4) with the kinetic energy defined as
$K\equiv(m_{0}^{2}c^{4}+p^{2}c^{2})^{1/2}-m_{0}c^{2}$, we find that protons
reach the knee at $K{=5\times}10^{15}$ eV after 657 BSA interactions without
energy losses. In case of faster shock speeds, $V_{u}=20,000$ km s-1 the final
energy would be reached after 334 BSA, while $V_{u}=40,000$ km s-1 requires
only 172 BSA interactions. We have assumed here that transmitted particles
could be scattered back to the shock by downstream turbulence, while particles
reflected upstream can encounter another shock front, or can be turned back by
upstream waves.

The physical picture of acceleration in the shock reference frame adopted in
this paper is easier to understand than the traditional (Fermi) analysis of
ion acceleration in the plasma reference frame, where the convection electric
field is removed. This approach ignores the fundamental laws of physics that
particle trajectories are determined by the Lorentz force equation (1), as
illustrated in Figs. 1-2, and the energization is determined by Eq. (2), and
not by the energy transformation between two inertial systems as in the
Fermi/DSA model. Furthermore, Fermi (1949) relies on the concept of
reflections by the magnetic mirror force, which is valid for particles with
$r_{c}<D$, and clearly not applicable for high-energy, large gyroradius
particles.

##  5 Conclusions

It is shown that ions and electrons with gyroradius $r_{c}\gg D$ are
accelerated by the convection electric field $E_{y}$ in a process described as
ballistic surfing acceleration. BSA operates outside of shocks on particles
from approximately 100 eV for protons up to very high energies observed in the
cosmic ray spectrum. BSA predicts a knee in the spectrum when the gyroradius
becomes comparable to the size of the shock, which determines the knee energy
given by Eq. (10). The spectral index in BSA is determined by the effective
shock compression $c_{B}$, and can accurately reproduce the observed index
$s\approx-2.5$ below the knee energy as well as a steeper spectrum,
$s\approx-3$, above the knee. It is demonstrated that the Fermi/DSA model,
which yields the spectral index (9), is inconsistent with the fundamental
equations of electrodynamics (1) and (2), and is therefore not valid. The
popularity of the Fermi model was unjustified, primarily due to the
coincidental agreement of one spectral index value with the correct BSA model
and with observations.

It is found that to reach the knee energy of 5$\times$1015 eV, a proton
starting from 10 keV in a collisionless environment needs only 334 BSA
interactions in shocks moving with velocity $V_{u}=20,000$ km s-1 . For
protons moving in the average magnetic field of $\langle B\rangle=$ 1 nT, the
net acceleration time of 334 gyroperiods would correspond to 227 years of
sidereal time, accounting for time dilation during each gyroperiod.

List of Acronyms
BSA – ballistic surfing acceleration
DSA – diffusive shock acceleration
QAH – quasi adiabatic heating
SDA – shock drift acceleration (a subset of QAH)
SSA – shock surfing acceleration
SWE – stochastic wave energization
TTT – transit time thermalization

## Data/software availability

The mathematical shock model used to make Figs. 1-3 was published in:
https://doi.org/10.1093/mnrasl/slad071.

## Acknowledgements

This work has been supported by Narodowe Centrum Nauki (NCN), Poland, through
grant No. 2021/41/B/ST10/00823.

## References

  * Badenes et al. (2010) Badenes, C., Maoz, D., Draine, B.T., 2010.  On the size distribution of supernova remnants in the Magellanic Clouds.  MNRAS 407, 1301–1313.  doi:10.1111/j.1365-2966.2010.17023.x.
  * Bell (1978) Bell, A.R., 1978.  The acceleration of cosmic rays in shock fronts - I.  MNRAS 182, 147–156.  doi:10.1093/mnras/182.2.147.
  * Blandford and Eichler (1987) Blandford, R., Eichler, D., 1987\.  Particle acceleration at astrophysical shocks: A theory of cosmic ray origin.  Physics Reports 154, 1 – 75.  doi:10.1016/0370-1573(87)90134-7.
  * Burch et al. (2016) Burch, J.L., Moore, R.E., Torbert, R.B., Giles, B., 2016\.  Magnetospheric multiscale overview and science objectives.  Space Sci. Rev. 199, 1–17.  doi:10.1007/s11214-015-0164-9.
  * Fermi (1949) Fermi, E., 1949.  On the origin of the cosmic radiation.  Phys. Rev. 75, 1169–1174.  doi:10.1103/PhysRev.75.1169.
  * Helder et al. (2012) Helder, E., Vink, J., Bykov, A., Ohira, Y., Raymond, J., Terrier, R., 2012\.  Observational signatures of particle acceleration in supernova remnants.  Space Sci. Rev. 173, 369–431.  doi:10.1007/s11214-012-9919-8.
  * Hillas (1984) Hillas, A.M., 1984.  The origin of ultra-high-energy cosmic rays.  Annual Review of Astronomy and Astrophysics 22, 425–444.  doi:10.1146/annurev.aa.22.090184.002233.
  * Hoshino and Shimada (2002) Hoshino, M., Shimada, N., 2002\.  Nonthermal electrons at high mach number shocks: Electron shock surfing acceleration.  The Astrophysical Journal 572, 880\.  doi:10.1086/340454.
  * Jokipii (1982) Jokipii, J.R., 1982.  Particle drift, diffusion and acceleration at shocks.  The Astrophysical Journal 255, 716–720.  doi:10.1086/159870.
  * Jones and Ellison (1991) Jones, F.C., Ellison, D.C., 1991\.  The plasma physics of shock acceleration.  Space Sci. Rev. 58, 259–346.  doi:10.1007/BF01206003.
  * Kennel (1988) Kennel, C.F., 1988.  Shock structure in classical magnetohydrodynamics.  Journal of Geophysical Research: Space Physics 93, 8545–8557.  doi:10.1029/JA093iA08p08545.
  * Lever et al. (2001) Lever, E.L., Quest, K.B., Shapiro, V.D., 2001.  Shock surfing vs. shock drift acceleration.  Geophysical Research Letters 28, 1367–1370.  doi:10.1029/2000GL012516.
  * Longair (2011) Longair, M.S., 2011.  High Energy Astrophysics.  University of Cambridge, Cambridge.
  * Shapiro et al. (2001) Shapiro, V.D., Lee, M.A., Quest, K.B., 2001.  Role of lower hybrid turbulence in surfing acceleration at perpendicular shocks.  Journal of Geophysical Research: Space Physics 106, 25023–25030.  doi:10.1029/1999JA000384.
  * Stasiewicz (2023a) Stasiewicz, K., 2023a.  Origin of flat-top electron distributions at the earth’s bow shock.  MNRAS 527, L71–L75.  doi:10.1093/mnrasl/slad146.
  * Stasiewicz (2023b) Stasiewicz, K., 2023b.  Transit time thermalization and the stochastic wave energization of ions in quasi-perpendicular shocks.  MNRAS 524, L50–L54.  doi:10.1093/mnrasl/slad071.
  * Stasiewicz and Eliasson (2021) Stasiewicz, K., Eliasson, B., 2021\.  Ion acceleration to 100 keV by the ExB wave mechanism in collisionless shocks.  MNRAS 508, 1888–1896.  doi:10.1093/mnras/stab2739.
  * Stasiewicz and Eliasson (2023) Stasiewicz, K., Eliasson, B., 2023\.  Electron heating mechanisms at the bow shock - revisited with Magnetospheric Multiscale measurements.  MNRAS 520, 3238–3244.  doi:10.1093/mnras/stad361.
  * Stasiewicz et al. (2021) Stasiewicz, K., Eliasson, B., Cohen, I.J., Turner, D.L., Ergun, R.E., 2021.  Local acceleration of protons to 100 keV by the ExB wave mechanism in a quasi-parallel bow shock.  Journal of Geophysical Research: Space Physics 126, e2021JA029477.  doi:10.1029/2021JA029477.
  * Stasiewicz and Kłos (2022) Stasiewicz, K., Kłos, Z., 2022\.  On the formation of quasi-parallel shocks, magnetic and electric field turbulence, and ion energisation mechanism.  MNRAS 513, 5892–5899.  doi:10.1093/mnras/stac1193.
  * Zank et al. (1996) Zank, G.P., Pauls, H.L., Cairns, I.H., Webb, G.M., 1996\.  Interstellar pickup ions and quasi-perpendicular shocks: Implications for the termination shock and interplanetary shocks.  Journal of Geophysical Research: Space Physics 101, 457–477.  doi:10.1029/95JA02860.