# Dark Solar Wind

Jae Hyeok Chang<EMAIL_ADDRESS>Department of Physics and Astronomy, Johns
Hopkins University, Baltimore, MD 21218, USA  Maryland Center for Fundamental
Physics, Department of Physics, University of Maryland, College Park, MD
20742, USA  David E. Kaplan<EMAIL_ADDRESS>Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  Surjeet
Rajendran<EMAIL_ADDRESS>Department of Physics and Astronomy, Johns Hopkins
University, Baltimore, MD 21218, USA  Harikrishnan Ramani
<EMAIL_ADDRESS>Stanford Institute for Theoretical Physics, Stanford
University, Stanford, CA 94305, USA  Erwin H. Tanin<EMAIL_ADDRESS>Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD
21218, USA

###### Abstract

We study the solar emission of light dark sector particles that self-interact
strongly enough to self-thermalize. The resulting outflow behaves like a fluid
which accelerates under its own thermal pressure to highly relativistic bulk
velocities in the solar system. Compared to the ordinary non-interacting
scenario, the local outflow has at least $\sim 10^{3}$ higher number density
and correspondingly at least $\sim 10^{3}$ lower average energy per particle.
We show how this generic phenomenon arises in a dark sector comprised of
millicharged particles strongly self-interacting via a dark photon. The
millicharged plasma wind emerging in this model has novel yet predictive
signatures that encourages new experimental directions. This phenomenon
demonstrates how a small step away from the simplest models can lead to
radically different outcomes and thus motivates a broader search for dark
sector particles.

††preprint: UMD-PP-022-05

Light particles with some coupling to the Standard Model (SM) can be produced
in the Sun with $\sim\text{keV}$ energies. The luminosity of such particles is
strongly limited by stellar cooling arguments Raffelt (1990). If these
particles simply free stream away as soon as they are produced in the Sun, the
outcome is an outflow of $\sim\text{keV}$ energy particles whose particle-
number flux is currently too low to be detected near the Earth. However,
interactions within the dark sector are in general poorly constrained, and, as
we will show, the story can change dramatically if one takes them into
account.

We focus on the predictive scenario where the interactions within the dark
sector allow these particles to locally thermalize via number-changing
processes. There are two natural outcomes of this scenario. (1) While the
dark-particle luminosity is still limited by cooling bounds, the self-
thermalization of these particles has the effect of enhancing the resulting
particle-number flux at the expense of lowering the average energy per
particle. Furthermore, (2) once local thermal equilibrium can be established,
the mean free path of these particles can be microscopically small and on
macroscopic scales they collectively display hydrodynamic behaviour which
further modifies the properties of the outflow. The actual dynamics is a
mixture of both effects and this gives rise to novel experimental and
astrophysical signatures.

As long as these particles are relativistic, their thermal pressure will
continually convert thermal energy into bulk fluid motion. In a way
mathematically analogous to the Parker solar wind model Parker (1965), this
eventually leads to a steady outflow of the fast dark-particle fluid, which we
refer to as the dark solar wind.

We elaborate these points in the remainder of this Letter. For concreteness,
we adopt a model of dark fermions interacting via dark photons as a
representative of a self-interacting dark sector. In order to have this sector
produced in the Sun, the dark photon is assumed to have a small kinetic mixing
with the SM photon Holdom (1986), thus making the dark fermions effectively
millicharged. This property also enables this dark solar wind to be detectable
on Earth.

### Model.—

Dark sector particles need to satisfy the following conditions to form a dark
solar wind:

  1. 1.

There exist a coupling between the dark sector particles and SM particles so
that they are produced in the Sun, and this coupling must be weak enough to
avoid the stellar cooling bound.

  2. 2.

Their self-interaction is sufficiently strong to allow them to self-thermalize
after being produced.

  3. 3.

They are light enough to remain relativistic even after self-thermalization.

There are various dark sector models to meet these conditions, but in this
paper we focus on a specific model, a dark photon and a dark fermion which
becomes a millicharged particle after photon-dark photon mixing.

Our Lagrangian is described by

$\displaystyle\mathcal{L}_{D}$ $\displaystyle=$
$\displaystyle-\frac{1}{4}F^{\prime}_{\mu\nu}F^{\prime\mu\nu}-\frac{\epsilon}{2}F^{\prime}_{\mu\nu}F^{\mu\nu}$
(1)
$\displaystyle+\,\,\bar{\chi}\left(i\gamma^{\mu}\partial_{\mu}+g_{D}\gamma^{\mu}A^{\prime}_{\mu}-m_{\chi}\right)\chi\,,$

where $\chi$ is a dark fermion with mass $m_{\chi}$, $A^{\prime}$ is the
massless dark photon, $\epsilon$ is the mixing angle between the SM photon and
the dark photon, and $g_{D}=\sqrt{4\pi\alpha_{D}}$ is the dark gauge coupling.
Once the kinetic terms are diagonalized, $\chi$ is effectively millicharged
and couples to the SM photon with electric charge $\epsilon g_{D}$. This model
has been studied in various contexts (see e.g. Alexander  _et al._ (2016);
Battaglieri  _et al._ (2017) for reviews), including the recent works
Arvanitaki  _et al._ (2021); Dvorkin  _et al._ (2019); Mathur  _et al._
(2020).

In the main part of the analysis, we limit ourselves to the parameter space
where $m_{\chi}$ is light enough that the resulting phenomenology is
equivalent to that of massless fermions. We will clarify the boundary of this
regime and comment on how the phenomenology would change for heavier
$m_{\chi}$.

### Production in the Sun.—

Due to their non-zero electric charge, $\chi$ particles can be pair created
from the hot and dense SM plasma inside the Sun. For a dark fermion $\chi$
with mass much smaller than the solar temperature ($m_{\chi}\ll T_{\odot}$),
the dominant production mechanism is via transverse plasmon decays. The
number-density production rate $\dot{n}$ and the power per unit volume
$\dot{Q}$ in the form of $\chi\bar{\chi}$ pairs through this channel were
found analytically in Chang  _et al._ (2021). The calculation details are
shown in the Supplemental Material. To obtain their numerical values, we pair
this analysis with the solar temperature and density profiles of Bahcall  _et
al._ (2005). For the crude estimates in this Letter, we will use the following
values of $\dot{n}$ at the center of the Sun and the $\chi$ production
luminosity

$\displaystyle\dot{n}_{c}$ $\displaystyle\equiv$ $\displaystyle\,\dot{n}(r=0)$
(2) $\displaystyle\sim$ $\displaystyle\,2\times
10^{9}\left(\frac{\epsilon}{10^{-15}}\right)^{2}\left(\frac{\alpha_{D}}{1}\right)\,\textrm{cm}^{-3}\textrm{s}^{-1}\,,$
$\displaystyle L_{\chi}$ $\displaystyle\equiv$
$\displaystyle\,{\displaystyle\int_{V_{\odot}}}\dot{Q}dV$ (3)
$\displaystyle\sim$ $\displaystyle\,8\times
10^{36}\left(\frac{\epsilon}{10^{-15}}\right)^{2}\left(\frac{\alpha_{D}}{1}\right)\,\textrm{MeV}\,\textrm{s}^{-1}\,.$

If the $\chi$ electric charge $\epsilon g_{D}$ is too large, the $\chi$
particles produced in stars carry away anomalously large amounts of energy
thereby changing the evolutionary history of the stars. Bounds from the non-
observation of such an anomalous evolution in red giants sets the most
stringent stellar bound on the $\chi$ electric charge: Vogel and Redondo
(2014)

$\epsilon\alpha_{D}^{1/2}\lesssim 2\times 10^{-15}$ (4)

### Self-thermalization.—

The initial population of $\chi$ particles produced in the Sun is in a state
far from (local) thermal equilibrium. Here, we work out a sufficient condition
for these particles to achieve thermalization.

The newly pair-created $\chi$ particles from plasmon decays in the Sun have a
“hard” energy spectrum, with a typical energy $E_{\rm hard}$ roughly given by
the temperature at the core of the Sun ($T_{\odot}$)

$E_{\rm hard}\sim T_{\odot}\sim 1\text{ keV}\,.$ (5)

Conservatively we start with the lowest possible abundance of these hard
particles arising from the free-streaming regime. Since these particles are
produced relativistically, they typically stay inside the solar core radius
$r_{\rm core}$, which we take as $r_{\rm core}\approx 0.2r_{\odot}$111Note
that the SM plasma, and hence $\dot{n}$, is roughly uniform within the Sun’s
core radius $r_{\rm core}\approx 0.2r_{\odot}$ and decays rapidly outside of
it. For this reason it is more appropriate to use $r_{\rm core}$ instead of
the full solar radius $r_{\odot}$ in our estimates., for a period of $\sim
r_{\rm core}$. Hence, their starting number density is

$\displaystyle n_{\rm hard}$ $\displaystyle\sim$
$\displaystyle\dot{n}_{c}r_{\rm core}$ (6) $\displaystyle\sim$ $\displaystyle
8\times
10^{8}\left(\frac{\epsilon}{10^{-15}}\right)^{2}\left(\frac{\alpha_{D}}{1}\right)\text{
cm}^{-3}$

These $E_{\rm hard}$ and $n_{\rm hard}\ll E_{\rm hard}^{3}$ are a much higher
average energy and a much lower number density compared to their would-be
thermal equilibrium values for the same energy density. In order to
thermalize, these particles must decrease their average energy and increase
their number density.

The high-energy, under-occupied initial state of the $\chi$ particles produced
in the Sun in our scenario resembles that of the products of perturbative
inflaton decay in the early universe, a well-studied scenario Davidson and
Sarkar (2000); Harigaya and Mukaida (2014); Mukaida and Yamada (2016). The
subsequent thermalization of such particles proceeds dominantly through
inelastic processes. While kinematics forbids a single $\chi$ particle from
spontaneously emitting a dark photon $\gamma_{D}$, particle production instead
proceeds through soft bremsstrahlung. The rate of such a $2\rightarrow 3$
process is roughly the rate of the enabling soft scattering process
$\Gamma_{2\rightarrow 2}^{\rm soft}$ multiplied by a factor of $\alpha_{D}$
for the $\gamma_{D}$ emission. Despite the extra $\alpha_{D}$ suppression, the
rate $\Gamma_{2\rightarrow 3}$ is enhanced due to the fact that it is
dominated by soft momentum exchanges, whose cross-section is large. This rate
$\Gamma_{2\rightarrow 3}$ is also limited by the formation time of the emitted
$\gamma_{D}$, i.e. the so-called Landau-Pomeranchuk-Migdal (LPM) effect Arnold
_et al._ (2001); Kurkela and Moore (2011), but even so it is still much faster
than other effects, such as the large-angle elastic scattering rate
$\Gamma_{2\rightarrow 2}^{\Delta\theta\sim 1}$.

The soft $2\rightarrow 2$ scattering rate is infrared divergent and given by

$\Gamma_{2\rightarrow 2}^{\rm soft}\sim\frac{\alpha_{D}^{2}n_{\rm
hard}}{k_{\rm min}^{2}}$ (7)

where $k_{\rm min}$ is the IR cutoff of the dark photon momentum, $k>k_{\rm
min}$. In our setup, $k_{\rm min}$ is set by the pre-thermalization dark-
sector Debye scale $\omega_{D}^{\rm pre}$, which is initially

$\displaystyle k_{\rm min}\sim\omega_{D}^{\rm pre}$ $\displaystyle\sim$
$\displaystyle\left(\frac{\alpha_{D}n_{\rm hard}}{E_{\rm hard}}\right)^{1/2}$
(8) $\displaystyle\sim$ $\displaystyle 7\times
10^{-5}\left(\frac{\epsilon}{10^{-15}}\right)\left(\frac{\alpha_{D}}{1}\right)\,\textrm{eV}$

A $\chi$ particle with an incoming momentum $p_{\rm in}\sim E_{\rm hard}$ can
scatter with another $\chi$ particle, become off-shell, and emit an extra
$\gamma_{D}$ with a momentum $k\lesssim p_{\rm in}$. The rate for such a
$2\rightarrow 3$ process is given by Garny  _et al._ (2019)

$\displaystyle\Gamma_{2\rightarrow 3}$ $\displaystyle\sim$
$\displaystyle\alpha_{D}\text{ min}\left(\Gamma_{2\rightarrow 2}^{\rm
soft},t_{\rm form}^{-1}\right)$ (9) $\displaystyle\sim$
$\displaystyle\alpha_{D}\Gamma_{2\rightarrow 2}^{\rm
soft}\text{min}\left(1,\sqrt{\frac{k}{k_{\rm LPM}}}\right)$

where $t_{\rm form}^{-1}$ sets an upper bound on the splitting rate due to the
fact that only one dark photon can be emitted in the timescale $t_{\rm form}$
it takes to resolve the dark photon. This leads to a suppression of the
splitting rate known as the LPM effect. In our case, $t_{\rm
form}\sim\sqrt{E_{\rm hard}^{2}/\alpha_{D}^{2}n_{\rm hard}k}$ Garny  _et al._
(2019) and on the second line we re-expressed it in terms of the threshold
momentum of the emitted $\gamma_{D}$ below which the LPM effect becomes
relevant

$k_{\rm LPM}\sim\frac{E_{\rm hard}^{2}\Gamma_{2\rightarrow
2}^{2}}{\alpha_{D}^{2}n_{\rm hard}}\sim\alpha_{D}\left(\frac{E_{\rm
hard}}{k_{\rm min}}\right)^{2}E_{\rm hard}\gg E_{\rm hard}$ (10)

The finding that $k_{\rm LPM}\gg E_{\rm hard}$ tells us that the LPM effect is
always important in our setup and the $2\rightarrow 3$ process is the fastest
for $k\sim E_{\rm hard}$, which boils down to the rate being
$\Gamma_{2\rightarrow 3}\sim\alpha_{D}^{3/2}\omega_{D}^{\rm pre}$. In the
subsequent $2\rightarrow 3$ splittings, $E_{\rm hard}$ goes down, $n_{\rm
hard}$ goes up, $k_{\rm min}$ goes up, which means $\Gamma_{2\rightarrow 3}$
will keep increasing. Through processes such as
$\gamma_{D}\chi\rightarrow\bar{\chi}\chi\chi$, the abundance of $\chi$ pairs
increases with the abundance of the dark photons $\gamma_{D}$ leading to an
acceleration to thermal equilibrium in the dark sector. Hence, the bottleneck
lies in the beginning and the requirement for achieving thermalization with
this process is $\Gamma_{2\rightarrow 3}^{\rm initial}r_{\rm core}\gtrsim 1$,
or

$\epsilon\alpha_{D}^{5/2}\gtrsim 2\times 10^{-26}$ (11)

We have identified a sufficient condition for the dark sector particles to
thermalize via bremsstrahlung processes. Other inelastic processes not
considered here are expected to bring the dark sector phase-space distribution
closer to thermalization. We have checked that elastic processes such as
scatterings within the dark sector, scatterings with SM particles, and the
bending of trajectories due to the solar magnetic field are negligible in the
pre-thermalization phase. Even if these effects were not weak they would only
help trap the $\chi$ particles, making it even easier for them to thermalize.

### Dark Solar Wind.—

If the dark sector particles manage to completely self-thermalize, their mean
free path would generically be far smaller than the length scales that dictate
their collective macroscopic dynamics222This is trivially satisfied for the
model (1) we are considering if the condition (11) is met. More specifically,
we have checked a posteriori that the mean free time of the dark particles as
seen in the Sun’s frame is easily much smaller than the timescale for, e.g.,
the density of particles to change by an $O(1)$ factor due to the fluid
expansion.. In that limit, such particles behave like a perfect fluid whose
properties can be described in a largely model-independent way. Given the
billion-year age of the Sun, it is likely that this fluid has relaxed by now
to a steady state, described by the following time-independent energy and
momentum equations333Numerical simulations Stellingwerf and Buff (1978); Keto
(2020) of analogous systems Parker (1965); Garlick (1979); Petterson  _et al._
(1980); Theuns and David (1992); Ray and Bhattacharjee (2002) suggest that the
steady state solution is an attractor and the relevant relaxation timescale is
given by the hydrodynamical timescale, which for our solar system is roughly
the flow time to 1 AU, namely $\sim 10$ minutes.

$\displaystyle\frac{1}{r^{2}}\partial_{r}\left[r^{2}\gamma^{2}v\left(\tilde{\rho}+\tilde{p}\right)\right]$
$\displaystyle=$ $\displaystyle\dot{Q}$ (12)
$\displaystyle\frac{1}{r^{2}}\partial_{r}\left[r^{2}\gamma^{2}v^{2}\left(\tilde{\rho}+\tilde{p}\right)\right]$
$\displaystyle=$ $\displaystyle-\partial_{r}\tilde{p}$ (13)

where $r$ is the radial position with respect to the center of the Sun,
$\gamma=(1-v^{2})^{-1/2}$ is Lorentz factor associated to the radial bulk
velocity $v$ of the fluid, $\tilde{\rho}(r)$ and $\tilde{p}(r)$ are the
comoving density and pressure of the fluid, and $\dot{Q}(r)$ is the power per
unit volume injection from the Sun in the form of $\chi\bar{\chi}$ pairs.
Quantities with a tilde $\tilde{}$ on top of it are defined in the rest frame
of the fluid and those without it are defined in Sun’s frame.

For simplicity, we assume here that the $\chi$ particles are massless or
sufficiently light that they are adequately described as a radiation-dominated
fluid with $\tilde{\rho}=3\tilde{p}=a\tilde{T}^{4}$, with $a$ a constant.
Integrating the energy equation (12) gives us the comoving temperature
$\tilde{T}$ of the fluid in terms of the fluid velocity $v$

$\frac{4a}{3}\tilde{T}^{4}=\frac{\int_{0}^{r}\dot{Q}(r^{\prime})4\pi r^{\prime
2}dr^{\prime}}{4\pi\gamma^{2}vr^{2}}$ (14)

Substituting the above into the momentum equation (13) results in an equation
for the fluid velocity only

$\left(\frac{1/3-v^{2}}{1/3+v^{2}}\right)\frac{\partial\ln v}{\partial\ln
r}=f(r)-\frac{2(1-v^{2})}{1+3v^{2}}$ (15)

where we defined a source function

$f(r)\equiv\frac{\dot{Q}(r)r^{3}}{\int_{0}^{r}\dot{Q}(r^{\prime})r^{\prime
2}dr^{\prime}}$ (16)

The first term on the RHS of (15) is due to the energy injection (i.e. inertia
injection) from the Sun, while the second term stems from the pressure
gradient of the fluid. Since the source $\dot{Q}(r)$ enters only via the
dimensionless quantity $f(r)$, the resulting velocity profile $v(r)$ does not
depend on the normalization of $\dot{Q}$, but only on the radial variation of
the quantity $f(r)$. In particular, for the model (1) under consideration
$f(r)$ is independent of both $\epsilon$ and $\alpha_{D}$.

Figure 1: Dark solar wind profiles as a function of the radial distance $r$
from the Sun. The Lorentz factor $\gamma(r)$ solves the fluid velocity
equation (15) with the transonic boundary condition (17). The source function
$f(r)\equiv\dot{Q}(r)r^{3}/\int_{0}^{r}\dot{Q}(r^{\prime})r^{\prime
2}dr^{\prime}$ is computed with the millicharged-particle energy density
injection rate $\dot{Q}$ of Chang  _et al._ (2021) and the numerical solar
profiles of Bahcall  _et al._ (2005). Both $\gamma(r)$ and $f(r)$ are
independent of $\epsilon$ and $\alpha_{D}$. The comoving temperature
$\tilde{T}(r)\propto(\epsilon^{2}\alpha_{D})^{1/4}$ of the fluid is evaluated
from the integrated energy equation (14) for $\epsilon=10^{-15}$ and
$\alpha_{D}=1$. Also shown are the location of the sonic point,
$r=0.24r_{\odot}$, and the Earth’s orbit radius, $r=215r_{\odot}$.

In order to solve (15) we need to specify some boundary conditions. Since the
Sun produces $\chi$ particles with no net radial momenta to begin with, $v$
must vanish at the origin. As we go to larger radii, there are two possible
types of profiles, depending on whether or not $v$ goes above the speed of
sound $c_{s}=1/\sqrt{3}$ during the outflow. Subsonic solutions444This
includes hydrostatic solutions., where $v<c_{s}$ all the way, predict
$v\propto r^{-2}$ at large distances from the Sun which, in turn, implies a
non-zero comoving temperature $\tilde{T}$ at infinity through (14). However,
the latter boundary condition is unphysical. Unless we add new ingredients to
the model that provide pressure support on the fluid far from the Sun, e.g.
through the high cosmic abundance of some particles interacting with the
fluid, we expect $\tilde{T}$ to vanish at large $r$. This leaves us with the
remaining possibility, namely the transonic solution, in which case the fluid
velocity increases from subsonic speeds ($v<c_{s}$) at small $r$ through the
sonic point ($v=c_{s}$) to supersonic speeds ($v>c_{s}$) at large $r$. The
monotonically increasing velocity implies a monotonically decreasing
$\tilde{T}$, thus satisfying the vanishing $\tilde{T}$ boundary condition at
infinity.

The sonic point ($v=1/\sqrt{3}$) can be crossed smoothly only if it coincides
with the zero of the RHS in (15) for $v=1/\sqrt{3}$. The latter most likely
occurs at around the Sun’s core radius $r_{\rm core}$, where $\dot{Q}$ starts
to drop rapidly, and numerically we found it to lie at $r=r_{\rm sonic}\approx
0.24r_{\odot}$. Thus, we require

$v=\frac{1}{\sqrt{3}}\quad\text{at}\quad r=r_{\rm sonic}\approx 0.24r_{\odot}$
(17)

This completes the boundary conditions for numerically solving the velocity
equation (15). We plot the resulting transonic Lorentz factor profile
$\gamma(r)$, together with the source function $f(r)$, and the comoving fluid
temperature $\tilde{T}$ from (14) in Fig. 1. It shows that beyond the sonic
point $r\gtrsim r_{\rm sonic}$ the Lorentz factor $\gamma$ of the fluid flow
asymptotes toward the well-known fireball solution Piran  _et al._ (1993),
$\gamma\sim r/r_{\rm sonic}$, for an adiabatically expanding fluid. In fact,
as shown in Ref. Piran  _et al._ (1993), when $\gamma\gg 1$ the fireball
solution $\gamma\sim r/r_{\rm sonic}$ solves not only the the sourceless
$\dot{Q}=0$ steady-state fluid equations (12) and (13) but also the time-
dependent fluid equations. Thus, we expect $\gamma\sim r/r_{\rm sonic}$ to
hold robustly outside the Sun even if for some reason the fluid flow deviates
from the assumed steady-state solution inside the Sun.

The structure of the fluid equation (15) considered here and the singling out
of the transonic solution are mathematically analogous to that of Parker’s
solar wind Parker (1965) (see also Bondi accretion Bondi (1952)). However, the
physical mechanisms behind them are completely different. Parker’s solar wind
is isothermal, non-relativistic, and accelerated by an interplay between
pressure gradient and gravity. On the other hand, the dark solar wind is
adiabatic, relativistic, and accelerated by an interplay between the pressure
gradient and energy/inertia injection from the Sun.

### Properties near the Earth.—

As the flow expands to larger $r$ and accelerates to higher Lorentz factors
$\gamma$, the comoving temperature of the fluid $\tilde{T}$ cools down
adiabatically according to (14). The integral in (14) for $r\gtrsim r_{\rm
core}$ yields the luminosity $L_{\chi}$ of the $\chi$ particles produced in
the Sun, resulting in $(4a/3)\tilde{T}^{4}\approx\gamma^{-2}(L_{\chi}/4\pi
r^{2})$ for $v\approx 1$. Interestingly, in the highly-relativistic limit
$\gamma\gg 1$ expected at $r\gg r_{\rm sonic}$ the scalings with $r$ of the
average energy per particle
$\left<E\right>\sim\gamma\tilde{T}\approx\text{const.}$ and number density
$n\sim\gamma\tilde{T}^{3}\propto r^{-2}$ in the Sun’s frame are identical to
those in the free streaming case, i.e. it is as if these particles simply free
streamed from the surface at which the fireball approximation starts to hold
($r\sim r_{\rm sonic}$). The latter is understandable because the acceleration
of the fluid to relativistic bulk velocities manifests itself at the particle
level as the velocities of the particles becoming increasingly radial the
farther they are from the Sun (relativistic beaming).

All things considered, the average energy $\left<E\right>$ and number density
$n$ of the dark particles at $r\gg r_{\rm sonic}$ in the Sun’s frame are given
up to $O(1)$ factors by

$\displaystyle\left<E\right>$ $\displaystyle\sim\left(\frac{L_{\chi}}{r_{\rm
sonic}^{2}}\right)^{1/4}\sim 1\text{
eV}\left(\frac{L_{\chi}}{10^{-2}L_{\odot}}\right)^{1/4}$ (18) $\displaystyle
n$ $\displaystyle\sim\frac{n_{\rm FS}\left<E\right>_{\rm
FS}}{\left<E\right>}\sim 10^{3}n_{\rm
FS}\left(\frac{L_{\chi}}{10^{-2}L_{\odot}}\right)^{-1/4}$ (19)

By contrast, in the free-streaming case the average energy per particle is
given by the core temperature of the Sun $\left<E\right>_{\rm
FS}\sim\text{keV}$ and if these particles are massless energy conservation
then gives $n_{\rm FS}\sim L_{\chi}/r^{2}\left<E\right>_{\rm FS}$. When the
dark particle luminosity saturates the cooling limit, $L_{\chi}\sim
10^{-2}L_{\odot}$, the dark solar wind gives $\sim 10^{3}$ lower average
energy $\left<E\right>$ and $\sim 10^{3}$ higher number density $n_{\chi}$
compared to those in the free-streaming case. The results in the two cases
deviate even more for $L_{\chi}\ll 10^{-2}L_{\odot}$. For the dark fermion
dark photon model (1) considered in this Letter, and with the $O(1)$ factors
included, the results at $r\gg r_{\rm sonic}$ are

$\displaystyle\gamma$ $\displaystyle\approx 893\left(\frac{r}{1\text{
AU}}\right)$ (20) $\displaystyle\tilde{T}$ $\displaystyle\approx
0.14\,\textrm{meV}\left(\frac{\epsilon}{10^{-15}}\right)^{\frac{1}{2}}\left(\frac{\alpha_{D}}{1}\right)^{\frac{1}{4}}\left(\frac{1\text{
AU}}{r}\right)$ (21) $\displaystyle\langle E\rangle$ $\displaystyle\approx
4\gamma\tilde{T}\approx
0.5\,\textrm{eV}\left(\frac{\epsilon}{10^{-15}}\right)^{\frac{1}{2}}\left(\frac{\alpha_{D}}{1}\right)^{\frac{1}{4}}$
(22) $\displaystyle n$
$\displaystyle=\frac{5\zeta(3)}{\pi^{2}}\gamma\tilde{T}^{3}\approx\frac{2\times
10^{5}}{\text{cm}^{3}}\left(\frac{\epsilon}{10^{-15}}\right)^{\frac{3}{2}}\left(\frac{\alpha_{D}}{1}\right)^{\frac{3}{4}}\left(\frac{1\text{
AU}}{r}\right)^{2}$ (23)

where we have used $\tilde{\rho}=a\tilde{T}^{4}$ with $a=(2+4\times
7/8)(\pi^{2}/30)$ corresponding to dark photons and $\chi\bar{\chi}$ pairs in
obtaining the above results, and $\left<E\right>$ was found by averaging over
all species using their spectra as seen on Earth. Naively, the relevant
spectrum of these particles on Earth would be that of a thermal distribution
boosted by $\gamma_{\text{AU}}\approx 893$. The spectrum that is actually seen
in an Earth-based detector is however biased by the fact that the detector
samples dark particles in an anisotropic way, e.g. more from the side facing
the incoming dark flow than from the opposite side. Accounting for both these
effects yields a distribution similar to that of a $\gamma_{\text{AU}}^{-2}$
rescaled thermal distribution Henry  _et al._ (1968)

$\frac{dn}{dEd\Omega}=\frac{1}{\gamma_{\text{AU}}^{2}}\frac{g_{*}}{(2\pi)^{3}}\frac{E^{2}}{e^{E/T_{\rm
eff}}\pm 1}$ (24)

with the degrees of freedom of the particle of interest $g_{*}$ and a
direction-dependent effective temperature

$T_{\rm eff}(\theta)\approx\frac{\tilde{T}}{\gamma_{\text{AU}}(1-v_{\rm
AU}\cos\theta)}$ (25)

where $v_{\rm AU}=(1-1/\gamma_{\rm AU}^{2})^{1/2}$, $d\Omega$ is a
differential solid angle, and $\theta$ is the angle of the particle velocity
relative to the local direction of the dark solar wind. Due to relativistic
beaming the angular distribution has a strong support only within
$\theta\lesssim\gamma_{\rm AU}^{-1}$.

Figure 2: Viable parameter space and predictions. The blue region violates the
stellar cooling limit (4). The red region does not satisfy the thermalization
requirement (11). The dashed lines correspond to different combinations of
dark solar wind density $n$ (23) and average energy per particle
$\left<E\right>$ (22) near the Earth.

### Detection.—

We next briefly discuss detection prospects of dark solar wind on Earth. In
principle, the dark photons as well as the dark fermions in the dark plasma
can be detected. Since the prospects are futuristic, we only provide order-of-
magnitude estimates and defer a systematic study for future work.

We start by considering experiments proposed for detecting non-relativistic
dark photons making up the dark matter via absorption Hochberg  _et al._
(2017); Bunting  _et al._ (2017); Knapen  _et al._ (2022); Mitridate  _et al._
(2021). A subset of these experiments are also sensitive to the relativistic
dark photons that have energies within the sensitive ranges of the detectors.
However, the absorption rate for the relativistic dark photons are suppressed
compared to their non-relativistic dark matter counterparts due to inferior
number density as well as suppression of the form
$\left(\alpha_{D}^{1/2}\tilde{T}/\langle E\rangle\right)^{\beta}$, where
$\beta=4(2)$ for the transverse (longitudinal) relativistic dark photons. The
latter comes from taking the relativistic limits of the dark-photon absorption
rates in An  _et al._ (2015). Due to this suppression, our preliminary
estimate suggests that the sensitivity of recently proposed dark-photon-
absorption experiments to dark solar wind is still a few orders of magnitude
above the cooling bound.

The dark fermions in the dark solar wind can potentially be probed through
their elastic or inelastic scattering with electrons in dark matter direct
detection experiments. In elastic scattering, since the dark fermions have
energies much lower than the mass of the electron, the energy deposited in the
electron is suppressed. The maximum average energies for the dark fermions in
this scenario are in the $\sim$ eV range resulting in energy depositions
$\lessapprox$ meV. This is lower than the planned thresholds of near-term
experiments but could potentially be detected in the future. Inelastic
processes wherein a bound electron is kicked out of its shell can
kinematically permit the dark fermion to lose $\mathcal{O}\left(1\right)$ of
its kinetic energy i.e. almost $\sim$ eV. While these processes are
kinematically more favorable, the cross-section for such inelastic processes
is suppressed by the momentum transferred during the process, requiring
detectors with larger target masses to probe phenomenologically interesting
parts of parameter space. We discuss these possibilities in detail in future
work.

The dark plasma flow studied here has properties that are very distinct from
that of cold dark matter, permitting detection strategies beyond scattering in
direct detection experiments. The uni-directional, relativistic and strongly-
coupled nature of the flow may allow us to probe this parameter space through
experiments in the same spririt as the so-called direct deflection Berlin  _et
al._ (2020); Berlin and Schutz (2021). The idea is to somehow perturb the dark
plasma wind and measure its back-reaction in the form of SM electric or
magnetic fields downstream. The generically tiny dark Debye length of the dark
plasma in our scenario quickly erases any static dark electric field in the
dark plasma once it is removed from the perturber. However, dark electric
currents may persist long enough to be detected Baym  _et al._ (1996). It
would be interesting to quantify this nontrivial dynamics with the help of
numerical simulations.

### Discussion.—

We pointed out a new generic phenomenon, here referred to as dark solar wind,
that arises in a generic light dark sector with sufficiently weak SM
interactions to avoid stellar cooling limits and sufficiently strong, number-
changing self-interactions to self-thermalize upon emission. Unlike in the
free-streaming scenario, the solar emission in this regime is less energetic,
denser, and behaves like a relativistically expanding fluid. Since the
properties of this fluid are dictated by thermal-equilibrium and steady-state
hydrodynamics, they are not sensitive to the details of the underlying
microphysics. We considered dark fermions charged under a dark photon that
kinetically mixes with the SM photon as an example, spelled out a sufficient
condition for achieving thermalization in the dark sector, numerically solved
the hydrodynamic equations for the resulting fluid, and worked out their
properties on Earth or elsewhere in the solar system. There might exist
parameter space in other self-interacting dark sector models such as light
non-Abelian sectors (motivated by potential solutions to the cosmological
constant problem Graham  _et al._ (2018, 2019); Berghaus  _et al._ (2021); Ji
_et al._ (2022)) and other models considered in Boddy  _et al._ (2014); March-
Russell  _et al._ (2020) which also exhibits dark solar wind.

Though we assumed that the dark fermion $\chi$ is massless in our discussions,
they are still valid for non-zero but light enough dark fermion mass
$m_{\chi}$. The thermalization condition is unaffected as long as $m_{\chi}$
is less than the lowest relevant energy scale, namely the pre-thermalization
debye frequency $\omega_{D}^{\rm pre}$, while the fluid dynamics is unchanged
as long as the dark fermions remain relativistic, i.e. $m_{\chi}\ll\tilde{T}$,
up to the radius of interest. As we increase $m_{\chi}$ from zero, we cross
$\omega_{D}^{\rm pre}$ way before the comoving temperature near the Earth
$\tilde{T}(r=1\text{ AU})$. We discuss this regime where the thermalization
condition is parameterically different but the fluid dynamics and hence the
model predictions are unchanged in the Supplemental Material. We leave the
explorations of yet higher $m_{\chi}$ regime as well as the effects of a non-
zero dark photon mass for future work.

We assumed in our analysis that the transonic steady-state solution that
solves the stationary fluid equation (15) is stable and has been reached by
now. We reiterate here that regardless of whether this stationary solution is
reached and maintained, our predictions would still be $O(1)$ correct since
the asymptotic fireball solution $\gamma\sim r/r_{\rm sonic}$ accurately
solves the time-dependent fluid equation outside the Sun. That said, it would
be interesting to run time-dependent hydrodynamical simulations in order to
confirm that this is indeed true and also to study the potential time
variability of the dark fluid properties.

The same SM couplings that produce the dark particles in the Sun could also
lead to a relic cosmological abundance of these particles Chu  _et al._
(2012); Evans  _et al._ (2020); Dvorkin  _et al._ (2019, 2021); Fernandez  _et
al._ (2021). Compared to the solar emissions, these cosmological relics are
more model-dependent as they rely on the early universe cosmology and physics
at higher energies. Nevertheless, we speculate on the impacts of such relics
here. Firstly, limits on dark matter self-couplings from the Bullet cluster
observation do not necessarily preclude the possibility that the same
strongly-coupled elementary particles produced in the Sun also make up all the
dark matter. Such astrophysical bounds can be relaxed or avoided if, for
instance, the cosmological relics come in the form of clumpy gas Buckley and
DiFranzo (2018) or composite bound states Cline  _et al._ (2014); Gresham  _et
al._ (2017). Even if these relics only make up a small fraction of the dark
matter, they can still have significant astrophysical impacts at galactic
scales Fan  _et al._ (2013). Furthermore, at the solar system scale, their
interactions with the dark solar wind may lead to a spectrum of outcomes. On
one extreme the dark matter may affect the dark solar wind by changing the
boundary conditions of the fluid equations, while on the other extreme the
dark solar wind may blow away and replace the typically assumed $\sim
0.3\text{ GeV}/\text{cm}^{3}$ dark matter energy density (or whatever its
expected relic density is) near the Sun.

Analogous fluid-like emissions of dark particles may also appear in a wide
range of astrophysical systems other than the Sun and be potentially
detectable on Earth. The dark wind emanating from the Earth’s core, while
closer in distance to terrestrial detection, is inferior in number density and
average energy in comparison to the dark solar wind due to far smaller dark
luminosity. Other astrophysical objects are of course farther from us than the
Sun, but their cores being much hotter and larger may yield new opportunities
for probing dark sector parameter space. Furthermore, their emissions may
differ from the solar one not only quantitatively but also qualitatively. For
instance, supernovae are transient and inject momentum in addition to energy
into the dark fluid, which might lead to completely different fluid solutions.
This work is thus a first step toward exploring more general dark fluid
outflows originating from a wide range of astrophysical systems.

### Acknowledgments.—

We thank William DeRocco, Peizhi Du, Rouven Essig, Anson Hook, Gustavo
Marques-Tavares, Mukul Sholapurkar, and Ken Van Tilburg for useful
discussions. J.H.C is supported by the NSF grant PHY-1914731, the Maryland
Center for Fundamental Physics, and the JHU Joint Postdoc Fund. D.E.K and S.R
are supported in part by the U.S. National Science Foundation (NSF) under
Grant No. PHY-1818899. This work was supported by the U.S. Department of
Energy (DOE), Office of Science, National Quantum Information Science Research
Centers, Superconducting Quantum Materials and Systems Center (SQMS) under
contract No. DE-AC02-07CH11359. S.R is also supported by the DOE under a
QuantISED grant for MAGIS, and the Simons Investigator Award No. 827042. H.R
acknowledges the support from the Simons Investigator Award 824870, DOE Grant
DE-SC0012012, NSF Grant PHY2014215, DOE HEP QuantISED award no. 100495, and
the Gordon and Betty Moore Foundation Grant GBMF7946.

## References

  * Raffelt (1990) G. G. Raffelt, Phys. Rept. 198, 1 (1990).
  * Parker (1965) E. Parker, Space Science Reviews 4, 666 (1965).
  * Holdom (1986) B. Holdom, Phys. Lett. B 166, 196 (1986).
  * Alexander  _et al._ (2016) J. Alexander _et al._ (2016) arXiv:1608.08632 [hep-ph] .
  * Battaglieri  _et al._ (2017) M. Battaglieri _et al._ , in  _U.S. Cosmic Visions: New Ideas in Dark Matter_ (2017) arXiv:1707.04591 [hep-ph] .
  * Arvanitaki  _et al._ (2021) A. Arvanitaki, S. Dimopoulos, M. Galanis, D. Racco, O. Simon, and J. O. Thompson, JHEP 11, 106 (2021), arXiv:2108.04823 [hep-ph] .
  * Dvorkin  _et al._ (2019) C. Dvorkin, T. Lin, and K. Schutz, Phys. Rev. D 99, 115009 (2019), arXiv:1902.08623 [hep-ph] .
  * Mathur  _et al._ (2020) A. Mathur, S. Rajendran, and E. H. Tanin, Phys. Rev. D 102, 055015 (2020), arXiv:2004.12326 [hep-ph] .
  * Chang  _et al._ (2021) J. H. Chang, R. Essig, and A. Reinert, JHEP 03, 141 (2021), arXiv:1911.03389 [hep-ph] .
  * Bahcall  _et al._ (2005) J. N. Bahcall, A. M. Serenelli, and S. Basu, Astrophys. J. Lett. 621, L85 (2005), arXiv:astro-ph/0412440 .
  * Vogel and Redondo (2014) H. Vogel and J. Redondo, JCAP 02, 029 (2014), arXiv:1311.2600 [hep-ph] .
  * Davidson and Sarkar (2000) S. Davidson and S. Sarkar, JHEP 11, 012 (2000), arXiv:hep-ph/0009078 .
  * Harigaya and Mukaida (2014) K. Harigaya and K. Mukaida, JHEP 05, 006 (2014), arXiv:1312.3097 [hep-ph] .
  * Mukaida and Yamada (2016) K. Mukaida and M. Yamada, JCAP 02, 003 (2016), arXiv:1506.07661 [hep-ph] .
  * Arnold  _et al._ (2001) P. B. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 11, 057 (2001), arXiv:hep-ph/0109064 .
  * Kurkela and Moore (2011) A. Kurkela and G. D. Moore, JHEP 12, 044 (2011), arXiv:1107.5050 [hep-ph] .
  * Garny  _et al._ (2019) M. Garny, A. Palessandro, M. Sandora, and M. S. Sloth, JCAP 01, 021 (2019), arXiv:1810.01428 [hep-ph] .
  * Stellingwerf and Buff (1978) R. F. Stellingwerf and J. Buff, Astrophys. J. 221, 661 (1978).
  * Keto (2020) E. R. Keto, Monthly Notices of the Royal Astronomical Society 493, 2834 (2020).
  * Garlick (1979) A. R. Garlick, aap 73, 171 (1979).
  * Petterson  _et al._ (1980) J. A. Petterson, J. Silk, and J. P. Ostriker, mnras 191, 571 (1980).
  * Theuns and David (1992) T. Theuns and M. David, Astrophys. J. 384, 587 (1992).
  * Ray and Bhattacharjee (2002) A. K. Ray and J. K. Bhattacharjee, Phys. Rev. E 66, 066303 (2002), arXiv:astro-ph/0212558 [astro-ph] .
  * Piran  _et al._ (1993) T. Piran, A. Shemi, and R. Narayan, Mon. Not. Roy. Astron. Soc. 263, 861 (1993), arXiv:astro-ph/9301004 .
  * Bondi (1952) H. Bondi, mnras 112, 195 (1952).
  * Henry  _et al._ (1968) G. R. Henry, R. B. Feduniak, J. E. Silver, and M. A. Peterson, Phys. Rev. 176, 1451 (1968).
  * Hochberg  _et al._ (2017) Y. Hochberg, T. Lin, and K. M. Zurek, Phys. Rev. D 95, 023013 (2017), arXiv:1608.01994 [hep-ph] .
  * Bunting  _et al._ (2017) P. C. Bunting, G. Gratta, T. Melia, and S. Rajendran, Phys. Rev. D 95, 095001 (2017), arXiv:1701.06566 [hep-ph] .
  * Knapen  _et al._ (2022) S. Knapen, J. Kozaczuk, and T. Lin, Phys. Rev. D 105, 015014 (2022), arXiv:2104.12786 [hep-ph] .
  * Mitridate  _et al._ (2021) A. Mitridate, T. Trickle, Z. Zhang, and K. M. Zurek, JHEP 09, 123 (2021), arXiv:2106.12586 [hep-ph] .
  * An  _et al._ (2015) H. An, M. Pospelov, J. Pradler, and A. Ritz, Phys. Lett. B 747, 331 (2015), arXiv:1412.8378 [hep-ph] .
  * Berlin  _et al._ (2020) A. Berlin, R. T. D’Agnolo, S. A. R. Ellis, P. Schuster, and N. Toro, Phys. Rev. Lett. 124, 011801 (2020), arXiv:1908.06982 [hep-ph] .
  * Berlin and Schutz (2021) A. Berlin and K. Schutz, (2021), arXiv:2111.01796 [hep-ph] .
  * Baym  _et al._ (1996) G. Baym, D. Bodeker, and L. D. McLerran, Phys. Rev. D 53, 662 (1996), arXiv:hep-ph/9507429 .
  * Graham  _et al._ (2018) P. W. Graham, D. E. Kaplan, and S. Rajendran, Phys. Rev. D 97, 044003 (2018), arXiv:1709.01999 [hep-th] .
  * Graham  _et al._ (2019) P. W. Graham, D. E. Kaplan, and S. Rajendran, Phys. Rev. D 100, 015048 (2019), arXiv:1902.06793 [hep-ph] .
  * Berghaus  _et al._ (2021) K. V. Berghaus, P. W. Graham, D. E. Kaplan, G. D. Moore, and S. Rajendran, Phys. Rev. D 104, 083520 (2021), arXiv:2012.10549 [hep-ph] .
  * Ji  _et al._ (2022) L. Ji, D. E. Kaplan, S. Rajendran, and E. H. Tanin, Phys. Rev. D 105, 015025 (2022), arXiv:2109.05285 [hep-ph] .
  * Boddy  _et al._ (2014) K. K. Boddy, J. L. Feng, M. Kaplinghat, and T. M. P. Tait, Phys. Rev. D 89, 115017 (2014), arXiv:1402.3629 [hep-ph] .
  * March-Russell  _et al._ (2020) J. March-Russell, H. Tillim, and S. M. West, Phys. Rev. D 102, 083018 (2020), arXiv:2007.14688 [astro-ph.CO] .
  * Chu  _et al._ (2012) X. Chu, T. Hambye, and M. H. G. Tytgat, JCAP 05, 034 (2012), arXiv:1112.0493 [hep-ph] .
  * Evans  _et al._ (2020) J. A. Evans, C. Gaidau, and J. Shelton, JHEP 01, 032 (2020), arXiv:1909.04671 [hep-ph] .
  * Dvorkin  _et al._ (2021) C. Dvorkin, T. Lin, and K. Schutz, Phys. Rev. Lett. 127, 111301 (2021), arXiv:2011.08186 [astro-ph.CO] .
  * Fernandez  _et al._ (2021) N. Fernandez, Y. Kahn, and J. Shelton, (2021), arXiv:2111.13709 [hep-ph] .
  * Buckley and DiFranzo (2018) M. R. Buckley and A. DiFranzo, Phys. Rev. Lett. 120, 051102 (2018), arXiv:1707.03829 [hep-ph] .
  * Cline  _et al._ (2014) J. M. Cline, Z. Liu, G. D. Moore, and W. Xue, Phys. Rev. D 90, 015023 (2014), arXiv:1312.3325 [hep-ph] .
  * Gresham  _et al._ (2017) M. I. Gresham, H. K. Lou, and K. M. Zurek, Phys. Rev. D 96, 096012 (2017), arXiv:1707.02313 [hep-ph] .
  * Fan  _et al._ (2013) J. Fan, A. Katz, L. Randall, and M. Reece, Phys. Dark Univ. 2, 139 (2013), arXiv:1303.1521 [astro-ph.CO] .
  * Braaten and Segel (1993) E. Braaten and D. Segel, Phys. Rev. D 48, 1478 (1993), arXiv:hep-ph/9302213 .
  * Peigne and Smilga (2009) S. Peigne and A. V. Smilga, Phys. Usp. 52, 659 (2009), arXiv:0810.5702 [hep-ph] .
  * Dokshitzer and Kharzeev (2001) Y. L. Dokshitzer and D. E. Kharzeev, Phys. Lett. B 519, 199 (2001), arXiv:hep-ph/0106202 .

Supplemental Material: Dark Solar Wind

##  I Plasma Effects in the Sun

The dominant production channel of $\chi\bar{\chi}$ pairs in the Sun is
through transverse plasmon decays, with the number density production rate
$\dot{n}$ and energy density production rate $\dot{Q}$ given by Chang  _et
al._ (2021)

$\displaystyle\dot{n}$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}g_{T}f_{\textrm{BE}}(\omega_{T})\Gamma_{\gamma^{*}_{T}\rightarrow\chi\bar{\chi}}$
(S1) $\displaystyle\dot{Q}$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}g_{T}f_{\textrm{BE}}(\omega_{T})\omega_{T}\Gamma_{\gamma^{*}_{T}\rightarrow\chi\bar{\chi}}\,,$
(S2)

where $g_{T}=2$ is the degree of freedom for transverse plasmon,
$f_{\textrm{BE}}(E)$ is the Bose-Einstein distribution, $\omega_{T}(k)$ is the
energy of the transverse plasmon, and

$\displaystyle\Gamma_{\gamma^{*}_{T}\rightarrow\chi\bar{\chi}}=\frac{\epsilon^{2}\alpha_{D}}{3}\frac{\omega_{T}^{2}-k^{2}}{\omega_{T}}\,.$
(S3)

In the core of the Sun, the dispersion relation of the photon is modified by
plasma effects. These effects are captured by the real parts of polarization
tensors of the photon, $\text{Re}\Pi_{T,L}$, which come from the self-energy
calculation with background electrons. In this work, we take the classical
limit in Braaten and Segel (1993), which gives

$\displaystyle\text{Re}\Pi_{T}$ $\displaystyle=$
$\displaystyle\omega_{p}^{2}\left(1+\frac{k^{2}}{\omega_{T}^{2}}\frac{T}{m_{e}}\right)\,,$
(S4) $\displaystyle\text{Re}\Pi_{L}$ $\displaystyle=$
$\displaystyle\omega_{p}^{2}\left(1+3\frac{k^{2}}{\omega_{L}^{2}}\frac{T}{m_{e}}\right)\,,$
(S5)

for transverse and longitudinal plasmons with a 4-momentum
$(\omega_{T,L},\vec{k})$ in a plasma with a temperature $T$. Here,
$\omega_{p}$ is the plasma frequency,

$\displaystyle\omega_{p}^{2}=\frac{4\pi\alpha
n_{e}}{m_{e}}\left(1-\frac{5}{2}\frac{T}{m_{e}}\right)\,,$ (S6)

$m_{e}$ is the mass of the electron, and $n_{e}$ is the electron density in
the plasma. Then the dispersion relation is given by

$\displaystyle\omega_{T}^{2}$ $\displaystyle=$ $\displaystyle
k^{2}+\text{Re}\Pi_{T}\,,\quad 0\leq k<\infty\,,$ (S7)
$\displaystyle\omega_{L}^{2}$ $\displaystyle=$
$\displaystyle\text{Re}\Pi_{L}\,,\quad 0\leq
k<\omega_{p}\sqrt{1+\frac{3T}{m_{e}}}\,.$ (S8)

For the plasmon decay to dark fermions, the transverse mode dominates because
the decay rate is proportional to the thermal mass,
$\sqrt{\omega_{T,L}^{2}-k^{2}}$, and the longitudinal thermal mass drops
quickly for large $k$. Also note we set the transverse vertex renormalization
constant to 1 ($Z_{T}=1$) for our calculations.

##  II Dark Sector Thermalization for Massive Dark Fermions

Though we assumed that $\chi$ is massless in deriving the thermalization
requirement in the main text, it is still valid as long as the $\chi$ mass is
below the pre-thermalization dark debye frequency,
$m_{\chi}\lesssim\omega_{D}^{\rm pre}$. Here, we generalize the discussion to
cases with more massive but still ultra-relativistic dark fermions $\chi$,
while keeping the dark photons $\gamma_{D}$ massless. As described in the main
text, the relaxation of the initially underpopulated $\chi$ particles produced
in the Sun toward thermal equilibrium relies on plasma-induced radiative
processes to increase the number of particles. Such processes have been
studied extensively in the QED literature, see Arnold  _et al._ (2001); Peigne
and Smilga (2009) for representative reviews on radiative loses in an abelian
plasma and also Mukaida and Yamada (2016); Garny  _et al._ (2019) which
analyzed a setup similar to ours, namely the thermalization an underoccupied
abelian plasma. The discussion below is an adaption of the results of these
earlier works in our context.

As we increase $m_{\chi}$ from zero, different effects become relevant as it
crosses the following key scales in the setup:

  * •

Pre-thermalization dark debye frequency inside the Sun

$\omega_{D}^{\rm pre}\sim\sqrt{\frac{\alpha_{D}n_{\rm hard}}{E_{\rm
hard}}}\sim 7\times
10^{-5}\left(\frac{\epsilon}{10^{-15}}\right)\left(\frac{\alpha_{D}}{1}\right)\text{
eV}$ (S9)

  * •

Post-thermalization dark debye frequency inside the Sun

$\omega_{D}^{\rm th}\sim\alpha_{D}^{1/2}\tilde{T}\sim(\alpha_{D}^{2}n_{\rm
hard}E_{\rm hard})^{1/4}\sim
0.3\left(\frac{\epsilon}{10^{-15}}\right)^{1/2}\left(\frac{\alpha_{D}}{1}\right)^{3/4}\text{
eV}$ (S10)

Note that we always have $\omega_{D}^{\rm th}\gg\omega_{D}^{\rm pre}$ because
$n_{\rm hard}\ll\tilde{T}^{3}$ and $E_{\rm hard}\gg\tilde{T}$.

  * •

The core temperature of the Sun $T_{\odot}\sim\text{ keV}$ multiplied with
$\alpha_{D}^{1/2}$

$\alpha_{D}^{1/2}T_{\odot}\sim
10^{3}\left(\frac{\alpha_{D}}{1}\right)^{1/2}\text{ eV}$ (S11)

which is always above $\omega_{D}^{\rm th}$ since we always have
$T_{\odot}\gg\tilde{T}$.

These scales divide the massive dark fermion case into several regimes, each
with a different expression for $\Gamma_{2\rightarrow 3}$. We discuss each of
these regimes below. The $m_{\chi}$ dependence of $\Gamma_{2\rightarrow 3}$
and the associated thermalization requirement $\Gamma_{2\rightarrow 3}r_{\rm
core}\gtrsim 1$ are illustrated in FIG. S1. Note that, in any of these
regimes, if $m_{\chi}\gtrsim\tilde{T}(1\text{ AU})$ dark fermions can
thermalize inside the Sun but will fall out of equilibrium (as they become
non-relativistic in the fluid frame) before arriving on Earth, hence changing
the predicted dark solar wind properties on Earth.

Figure S1: Left: The pre-thermalization
$\chi\chi\rightarrow\chi\chi\gamma_{D}$ splitting rate $\Gamma_{2\rightarrow
3}$ (up to $O(1)$ factors) as a function of the dark fermion mass $m_{\chi}$.
This splitting process is efficient at thermalizing the dark sector if
$\Gamma_{3\rightarrow 2}r_{\rm core}\gtrsim 1$, with $r_{\rm core}\approx
0.2r_{\odot}$. Right: The viable parameter space for different values of
$m_{\chi}$. Finite-mass effects lead to parameterically smaller
$\Gamma_{2\rightarrow 3}$ which in turn makes the thermalization requirement
more restrictive.

  1. 1.

$m_{\chi}\lesssim\omega_{D}^{\rm pre}$
We begin by reviewing the $m_{\chi}\ll\omega_{D}^{\rm pre}$ case in order to
highlight the assumptions that led to the $\Gamma_{2\rightarrow
3}\sim\alpha_{D}^{3/2}\omega_{D}^{\rm pre}$ expression used in the main text.
The sparse high-energy dark fermions $\chi$ produced in the Sun lose their
energies dominantly though near-collinear
$\chi\chi\rightarrow\chi\chi\gamma_{D}$ splittings (and their variants). The
rate for such $2\rightarrow 3$ splittings is naively given by
$\alpha_{D}\Gamma_{2\rightarrow 2}^{\rm soft}$, where $\Gamma_{2\rightarrow
2}^{\rm soft}\sim\alpha_{D}^{2}/k_{\rm min}^{2}$ is the rate of elastic
fermion scatterings that is infrared-regulated by the softest dark photon
momentum scale, $k_{\rm min}\sim\omega_{D}^{\rm pre}$. However, the actual
rate is suppressed by the so called LPM effect, which is due to the finite
time $t_{\rm form}$ it takes to resolve the emitted dark photon (daughter)
from the emitting dark fermion (parent). Dark photons emitted within $t_{\rm
form}$ interfere destructively and so only one dark photon can be emitted
effectively during $t_{\rm form}$. The splitting rate can therefore be written
as

$\Gamma_{\rm 2\rightarrow 3}\sim\alpha_{D}\text{
min}\left(\Gamma_{2\rightarrow 2}^{\rm soft},t_{\rm
form}^{-1}\right),\quad\text{for}\quad m_{\chi}\lesssim\omega_{D}^{\rm pre}$
(S12)

The formation time $t_{\rm form}$ can be understood as the time it takes for
the daughter particle to acquire an $O(1)$ phase in the parent’s frame. This
is especially relevant in our case because the separation between the parent
and daughter particles relies on either their same-direction velocity
difference $1-v_{\chi}\sim(m_{\chi}/E_{\chi})^{2}$, which is suppressed for a
highly relativistic parent $\chi$, or their relative transverse velocity,
which is (initially) suppressed due to the collinearity of the emission. In an
abelian plasma, the parent-daughter relative transverse velocity is dominated
by the deflection of the parent fermion555By constrast, in a non-abelian
plasma, the parent-daughter relative transverse velocity is dominated by the
deflection of the daughter gauge boson, which can scatter with other gauge
bosons in the plasma owing to their non-abelian nature Kurkela and Moore
(2011)., whose transverse momentum grows as
$p_{\perp}\sim\sqrt{N}\omega_{D}^{\rm pre}$ as it undergoes $\Delta
p_{\perp}\sim\omega_{D}^{\rm pre}$ soft collisions $N\sim\Gamma_{2\rightarrow
2}^{\rm soft}t_{\rm form}$ times in the medium. Taking into account both
effects, the reciprocal formation time is given by Peigne and Smilga (2009);
Garny  _et al._ (2019)

$t_{\rm form}^{-1}\sim k\,\text{max}\left[\sqrt{\frac{\alpha_{D}^{2}n_{\rm
hard}}{E_{\rm hard}^{2}k}},\left(\frac{m_{\chi}}{E_{\rm
hard}}\right)^{2}\right]$ (S13)

The first factor inside the square brackets corresponds to $t_{\rm form}$
being determined by the transverse velocity difference and is the relevant one
in the main text. The second factor accounts for the same-direction velocity
difference, which does not become important until
$m_{\chi}\gtrsim\omega_{D}^{\rm th}$. Furthermore, as along as
$m_{\chi}\lesssim\alpha_{D}^{1/2}T_{\odot}$ the rate $\Gamma_{2\rightarrow 3}$
is set by $t_{\rm form}^{-1}$ which is maximized for $k\sim E_{\rm hard}$,
resulting in the following expression used in the main text

$\Gamma_{2\rightarrow 3}\sim\alpha_{D}t_{\rm
form}^{-1}\sim\alpha_{D}^{3/2}\omega_{D}^{\rm pre},\quad\text{for}\quad
m_{\chi}\lesssim\omega_{D}^{\rm pre}$ (S14)

  2. 2.

$\omega_{D}^{\rm pre}\lesssim m_{\chi}\lesssim\omega_{D}^{\rm th}$
Finite-mass effects begin to matter when $m_{\chi}$ crosses the pre-
thermalization Debye frequency $\omega_{D}^{\rm pre}$ shown in (8). When
$m_{\chi}$ is above this scale, the squared amplitude for the $\gamma_{D}$
bremsstrahlung of a scattered $\chi$ particle is effectively suppressed by a
factor of $(\omega_{D}^{\rm pre}/m_{\chi})^{2}$, and consequently the
$2\rightarrow 3$ rate is now modified as Peigne and Smilga (2009)

$\Gamma_{\rm 2\rightarrow 3}\sim\left[\alpha_{D}\left(\frac{\omega_{D}^{\rm
pre}}{m_{\chi}}\right)^{2}\right]t_{\rm
form}^{-1}\sim\left(\frac{\omega_{D}^{\rm
pre}}{m_{\chi}}\right)^{2}\alpha_{D}^{3/2}\omega_{D}^{\rm
pre},\quad\text{for}\quad\omega_{D}^{\rm pre}\lesssim
m_{\chi}\lesssim\omega_{D}^{\rm th}$ (S15)

This suppression can be traced back to suppression in the emission amplitude,
known as the dead cone effect Dokshitzer and Kharzeev (2001), for dark photon
emission angles $\theta\lesssim m_{\chi}/E_{\rm hard}$.

  3. 3.

$\omega_{D}^{\rm th}\lesssim m_{\chi}\lesssim\alpha_{D}^{1/2}T_{\odot}$
The next scale that $m_{\chi}$ crosses is $(\alpha_{D}^{2}n_{\rm hard}E_{\rm
hard})^{1/4}$ which is parameterically equivalent to the post-thermalization
Debye frequency. Since we are only concerned with the pre-thermalization state
of the dark plasma here, the appearance of $\omega_{D}^{\rm th}$ in this
discussion is just a mathematical coincidence. In this regime, the reciprocal
formation time $t_{\rm form}^{-1}$ is set by its finite-mass lower bound,
$t_{\rm form}^{-1}\gtrsim k(m_{\chi}/E_{\rm hard})^{2}$. This has the effect
of parametrically enhancing $\Gamma_{2\rightarrow 3}$ relative to the
expression applicable in the $\omega_{D}^{\rm pre}\lesssim
m_{\chi}\lesssim\omega_{D}^{\rm th}$ regime. The $\Gamma_{2\rightarrow 3}$ in
this case is still maximized at $k\sim E_{\rm hard}$ with the value

$\Gamma_{2\rightarrow 3}\sim\left[\alpha_{D}\left(\frac{\omega_{D}^{\rm
pre}}{m_{\chi}}\right)^{2}\right]t_{\rm
form}^{-1}\sim\,\left(\frac{m_{\chi}}{\omega_{D}^{\rm
th}}\right)^{2}\left(\frac{\omega_{D}^{\rm
pre}}{m_{\chi}}\right)^{2}\alpha_{D}^{3/2}\omega_{D}^{\rm
pre},\quad\text{for}\quad\omega_{D}^{\rm th}\lesssim
m_{\chi}\lesssim\alpha_{D}^{1/2}T_{\odot}$ (S16)

which is independent of $m_{\chi}$.

  4. 4.

$m_{\chi}\gtrsim\alpha_{D}^{1/2}T_{\odot}$
Then, at $m_{\chi}\gtrsim\alpha_{D}^{1/2}T_{\odot}$ the reciprocal formation
time $t_{\rm form}^{-1}$ goes above the soft $2\rightarrow 2$ scattering rate
$\Gamma_{2\rightarrow 2}^{\rm soft}$, which means the LPM effect stops
operating, and

$\Gamma_{2\rightarrow 3}\sim\left[\alpha_{D}\left(\frac{\omega_{D}^{\rm
pre}}{m_{\chi}}\right)^{2}\right]\Gamma_{2\rightarrow 2}^{\rm
soft}\sim\left(\frac{\alpha_{D}^{1/2}T_{\odot}}{m_{\chi}}\right)\left(\frac{\omega_{D}^{\rm
pre}}{m_{\chi}}\right)\alpha_{D}^{3/2}\omega_{D}^{\rm
pre},\quad\text{for}\quad m_{\chi}\gtrsim\alpha_{D}^{1/2}T_{\odot}$ (S17)

However, it turns out that we always have
$\alpha_{D}^{1/2}T_{\odot}\gtrsim\tilde{T}(1\text{ AU})$ in the viable
parameter space where the dark sector particles satisfy both the cooling limit
and self-thermalization condition inside the Sun. Therefore, if $m_{\chi}$ was
high enough to be in this $m_{\chi}\gtrsim\alpha_{D}^{1/2}T_{\odot}$ regime,
the steady-state fluid properties inside the Earth’s orbit would have changed
relative to those derived in the main text.


Corresponding author: Lekan Molu (e-mail: lekanmolu@microsoft.com).

# Comments on “Time-Varying Lyapunov Functions for Tracking Control of
Mechanical Systems With and Without Frictions”

OLALEKAN OGUNMOLU1  Microsoft Research, 300 Lafayette Street, New York, NY
10012, USA.

###### Abstract

In the article111Ren, W., Zhang, B, Li, H, and Yan L. _Time-Varying Lyapunov
Functions for Tracking Control of Mechanical Systems With and Without
Frictions_. IEEE Access, vol. 8. pp. 51510-51517. 2020., the authors
introduced a time-varying Lyapunov function for the stability analysis of
nonlinear systems whose motion is governed by standard Newton-Euler equations.
The authors established asymptotic stability with the choice of two symmetric
positive definite matrices restricted by certain eigenvalue bounds in the
control law. Exponential stability in the sense of Lyapunov using integrator
backstepping and Lyapunov redesign is established in this note using just one
matrix in the derived controller. We do not impose minimum eigenvalue bound
requirements on the symmetric positive definite matrix introduced in our
analysis to guarantee stability. Reducing the parameters needed in the control
law, our analysis improves the stability and convergence rates of tracking
errors reported in the article1.

###### Index Terms:

Nonlinear mechanical systems, time-varying Lyapunov functions, tracking
control.

=-15pt

The well-established method of Lyapunov analysis is a great tool for
synthesizing stable controllers for nonlinear systems such as robot
manipulators, space flexures, parallel robots, and magnetic bearings inter
alia. The kinetics of such systems are typically derived from (time-invariant)
Euler-Lagrange equations. When the underlying nonlinear system is time-
varying, however, stability is difficult to establish with traditional time-
invariant Lyapunov functions.

Considering frictionless mechanical systems, the authors of the article [1]
(i) proposed a time-varying Lyapunov function for nonlinear systems;  (ii)
established sufficient conditions for their practical stability; and  (iii)
under a zero friction derived an asymptotic convergence rate for the filtered
velocity and position tracking errors.  In this note, we will prove
exponential error convergence when friction is absent from the dynamics. In
addition, only one feedback matrix (here the symmetric positive definite
matrix $P$) is needed for the establishment of our results – as opposed to two
symmetric positive definite matrices ($P$ and $K$) reported by the authors of
[1].

The non-negative function used to prove asymptotic stability in the paper [1]
is given as

$\displaystyle
Q(t)=\dfrac{1}{2}q_{r}(t)^{T}M(q)q_{r}(t)+\dfrac{1}{2}\tilde{q}^{T}(t)P\tilde{q}(t)$
(1)

where all variables are as defined in the paper [1] i.e. $q$ is the joint
space variable, $q_{r}$ is the filtered velocity tracking error, defined as

$\displaystyle q_{r}(t)=\dot{\tilde{q}}(t)+{\lambda}\tilde{q}(t),$ (2)

for a ${\lambda}>0$; $\tilde{q}(t)$ is the position tracking error given as
$\tilde{q}(t)=q_{d}(t)-q(t)$, for a desired joint position, $q_{d}$. The
matrix $M(q)$ is the symmetric positive definite inertia matrix of the system.

The model of the system to be controlled is given by

$M(q)\ddot{q}+C(q,\dot{q})\dot{q}+G(q)+F(t,\dot{q})=\tau$ (3)

where $C(q,\dot{q})$ denotes the Coriolis and centrifugal forces, $G(q)$
models the gravity and non-conservative external forces, $F(t,\dot{q})$ is the
time-varying friction, and $\tau$ is the vector of actuator torques.

Observe: The matrix $(\dot{M}(q)-2C(q,\dot{q}))$ is skew-symmetric for a
suitable choice of $C(q,\dot{q})$. One may also write this property as
$\dot{M}(q)=C(q,\dot{q})+C^{T}(q,\dot{q})$.

######  Lemma 1 (Asymptotic stability without frictions)

For the nonlinear manipulator (3) with a zero time-varying friction i.e.
$F(t,\dot{q})=0$, every motion of the mechanical system is asymptotically
stable when $F(t,\dot{q})\equiv 0$, i.e. $q(t)\rightarrow q_{d}(t)$ as
$t\rightarrow\infty$ provided we choose a virtual input, $\eta(t)$, which is
backstepped as an integrator as follows

$\displaystyle\eta(t)$ $\displaystyle=\dot{q}(t),$ (4a)
$\displaystyle\triangleq\dot{q}_{d}(t)-\tilde{q}(t)-\dfrac{1}{2}M^{-1}(q)\dot{M}(q)\tilde{q}(t).$
(4b)

###### Proof 1

Define the tracking error dynamics for a desired joint space variable,
$q_{d}$, as

$\displaystyle\dot{\tilde{q}}(t)$ $\displaystyle=\eta(t)-\dot{q}_{d}(t).$ (5)

Now, consider the following Lyapunov function candidate

$\displaystyle V_{1}(t)=\frac{1}{2}\tilde{q}^{T}(t){M(q)}\tilde{q}(t).$ (6)

We find that

$\displaystyle\dot{V}_{1}(t)$
$\displaystyle=\tilde{q}^{T}(t){M}(q)\dot{\tilde{q}}(t)+\dfrac{1}{2}\tilde{q}^{T}(t)\dot{M}(q)\tilde{q}(t),$
(7)
$\displaystyle=\tilde{q}^{T}(t){M}(q)(\eta(t)-\dot{q}_{d}(t))+\dfrac{1}{2}\tilde{q}^{T}(t)\dot{M}(q)\tilde{q}(t),$
(8) $\displaystyle=-\tilde{q}^{T}(t){M}(q)\tilde{q}(t).$ (9)

Equation (9) implies that the system must be asymptotically stable [2].

######  Theorem 1 (Exponential stability of the system)

Consider the entire nonlinear system (3), suppose that there exists a
$P=P^{T}\,\succ 0$, then trajectories that emanate from the system are
exponentially stable if we choose the control law, $\tau$, according to

$\displaystyle\tau$
$\displaystyle={M}(q)\left(\ddot{q_{d}}+\lambda\dot{\tilde{q}}+q_{r}\right)+C(q,\dot{q})\left(\dot{q}+q_{r}\right)+G(q).$
(10)

###### Proof 2

Consider the Lyapunov function candidate

$\displaystyle V(t)=V_{1}(t)+\frac{1}{2}q_{r}^{T}(t)PM(q)q_{r}(t)$ (11)

so that the time derivative of $V(t)$ (dropping the r.h.s templated arguments
for notation terseness) is

$\displaystyle\dot{V}(t)=\dot{V}_{1}(t)+q_{r}^{T}PM\dot{q}_{r}+\frac{1}{2}q_{r}^{T}P\dot{M}q_{r},$
(12a)
$\displaystyle=\dot{V}_{1}(t)+q_{r}^{T}PM(\ddot{\tilde{q}}+\lambda\dot{\tilde{q}})+\frac{1}{2}q_{r}^{T}P\dot{M}q_{r},$
(12b)
$\displaystyle=\dot{V}_{1}(t)+q_{r}^{T}PM(\ddot{q}_{d}-\ddot{q}+\lambda\dot{\tilde{q}})+\frac{1}{2}q_{r}^{T}P\dot{M}q_{r},$
(12c)
$\displaystyle=\dot{V}_{1}(t)+q_{r}^{T}PM(\ddot{q}_{d}+\lambda\dot{\tilde{q}})$
$\displaystyle\qquad\qquad+q_{r}^{T}P\left[C(q,\dot{q})(\dot{q}+{q}_{r})+G(q)-\tau\right].$
(12d)

Where we have used the skew-symmetric property $\dot{M}-2C$ in arriving at
(12d). Let $\tau$ be as defined in (10), then we must have,

$\displaystyle\dot{V}(t)$
$\displaystyle=-\tilde{q}(t)M(q)\tilde{q}(t)-q_{r}^{T}(t)PM(q)q_{r}(t),$ (13a)
$\displaystyle\dot{V}(t)$ $\displaystyle\triangleq-2V(t).$ (13b)

A fortiori, exponential stability follows [2] as a result of (13b).

Through a careful choice of virtual input design, the convergence rate of
system tracking errors can be significantly improved upon. In addition, the
number of parameters introduced into the torque control equation is reduced
from two to one. This is particularly of importance when designing fast real-
time learning-based optimization for suitable controllers for nonlinear
systems for which equation (3) is applicable.

## References

  * [1] Ren, W., Zhang, B, Li, H, and Yan L. _Time-Varying Lyapunov Functions for Tracking Control of Mechanical Systems With and Without Frictions_. IEEE Access, vol. 8. pp. 51510-51517. 2020.
  * [2] Ioannou, Petros A., and Jing Sun. _Robust adaptive control._ Courier Corporation, 2012.