# A composite solution to the EDGES anomaly

Anubhav Mathur  Department of Physics and Astronomy, Johns Hopkins University,
3400 N. Charles St., Baltimore, MD 21218, USA  Surjeet Rajendran  Department
of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles St.,
Baltimore, MD 21218, USA  Harikrishnan Ramani  Stanford Institute for
Theoretical Physics, Stanford University, Stanford, CA 94305, USA

(September 2021)

###### Abstract

Subcomponent millicharged dark matter that cools baryons via Coulomb
interactions has been invoked to explain the EDGES anomaly. However, this
model is in severe tension with constraints from cosmology and stellar
emissions. In this work, we consider the consequences of these millicharged
particles existing in composite Q-balls. The relevant degrees of freedom at
high temperature are minuscule elementary charges, which fuse at low
temperatures to make up Q-balls of larger charge. These Q-balls serve as the
degrees of freedom relevant in cooling the baryons sufficiently to account for
the EDGES anomaly. In such a model, cosmology and stellar constraints (which
involve high-temperature processes) apply only to the feebly-interacting
elementary charges and not to the Q-balls. This salvages a large range of
parameter space for millicharged Q-balls that can explain the EDGES anomaly.
It also opens up new parameter space for direct detection, albeit at low
momentum transfers.

##  I Introduction

The Experiment to Detect the Global Epoch of Re-ionization Signature (EDGES)
experiment has reported a dip in the 21 cm spectrum corresponding to strong
absorption around $z=17$ Bowman  _et al._ (2018). This can be interpreted as a
3.8 $\sigma$ deviation from the $\Lambda$CDM prediction for the baryon
temperature Barkana (2018); Barkana  _et al._ (2018); Kovetz  _et al._ (2018).
Dark matter (DM) cooling of the baryonic fluid has been invoked as an
explanation for this excess Barkana (2018); Barkana  _et al._ (2018); Kovetz
_et al._ (2018). A DM model that maximizes the cross-section around cosmic
dawn is sub-component millicharged dark matter (mCDM), which has a larger
cross-section with Standard Model (SM) charges at the lowest relative
velocities. However, the millicharge parameter space is extremely constrained
due to limits from CMB and BBN, cooling of SN1987A and stars, and terrestrial
experiments Kovetz  _et al._ (2018). It has subsequently been shown that even
this limited parameter space results in over-production of mCDM through
freeze-in Berlin  _et al._ (2018); Creque-Sarbinowski  _et al._ (2019).

These difficulties have led to two other ways to solve the EDGES anomaly. The
first involves heating the CMB relative to baryons Feng and Holder (2018);
Ewall-Wice  _et al._ (2018); Pospelov  _et al._ (2018); Fialkov and Barkana
(2019), while the second involves mCDM which is tightly coupled to an
additional cold component that forms the dominant DM which does the bulk of
the cooling Liu  _et al._ (2019). In this paper, we point to a third
possibility. The mCDM explanations for the EDGES anomaly to date have treated
the millicharged particle (mCP) as elementary without internal structure. As a
result, the same mCP is the physical particle at all energies. In this work,
we explore the consequence of this mCP being a composite state of elementary
mCPs with much smaller mass and charge, glued together by a force that
confines at low temperatures. The elementary charges are the relevant degrees
of freedom at temperatures and energies much higher than cosmic dawn. As a
result, in our model, constraints from CMB, BBN, over-closure, stellar and SN
cooling as well as colliders all apply only to the elementary charges. We
demonstrate here that there is a drastic increase in new parameter space for
mCDM as long as it is in a composite state. Furthermore, we explore the unique
thermal history for the dark sector that involves confinement when $T_{D}$
(the dark temperature) falls below $\Lambda_{D}$ (the dark confining scale)
and deconfinement if the dark temperature increases subsequently. We point out
a novel dark phase where thermal contact with the SM results in a thermostatic
dark bath, i.e. the dark bath staying at the same temperature with the heat
dump from baryons exactly cancelled by Hubble cooling.

##  II Model

We consider mCP scalars $s$ that carry electric charge $\epsilon_{s}$ and mass
$m_{s}$111As we explain later, fermionic mCPs share many of the properties of
scalar composites discussed in this work via nucleosynthesis. However, the
Fermi temperature of the constituents exceeds the cold temperatures required
to explain the EDGES anomaly.. The scalar is assumed to have quartic and
effective sextic self-interactions. These self-interactions make it possible
for the scalars to form Q-balls i.e.  these Q-balls will be the composite
states of interest to us. We also assume that the mCPs are charged under a
dark $U(1)$ with a dark charge $g_{D}$. The interactions of the scalar field
$s$ are captured in the Lagrangian

$\mathcal{L}=\frac{m_{s}^{2}|s|^{2}}{2}-\lambda|s|^{4}+\frac{1}{\Lambda_{6}}|s|^{6}+D_{\mu}sD^{\mu}s^{*}+\epsilon
F^{\mu\nu}F^{\prime}_{\mu\nu}$ (1)

Here $D_{\mu}=\partial_{\mu}-ig_{D}A^{\prime}_{\mu}$, and $F$ and $F^{\prime}$
are the field strength tensors of the SM U(1) and the dark U(1) i.e.
$A^{\prime}$ respectively. This dark $U(1)$ accomplishes two goals: first, it
allows the electric charge of the mCP to be generated via kinetic mixing and
second, it provides Coulomb repulsion that prevents the Q-balls from getting
too big. The mass of the dark photon is unimportant as long as it is
sufficiently long ranged to allow the Q-balls to interact with baryons via
Coulomb scattering. We can take this mass to be zero, or low enough to evade
direct stellar constraints on the dark photon.

For simplicity we will assume that these elementary mCPs that confine all have
the same charge sign. For net charge neutrality, we envision an asymmetric
dark component with the opposite charge, just like the standard model. Under
the correct circumstances, the scalars fuse to form composite Q-balls or
Q-balls $\chi$ with “atomic number” Q such that they carry charge
$\epsilon_{\chi}=Q\epsilon_{s}$. The mass of these Q-balls is given by

$m_{\chi}=Q\sqrt{m_{s}^{2}-\frac{\lambda^{2}\Lambda_{6}^{2}}{2}}=Q\Lambda_{D}$
(2)

i.e. the difference between the total mass square of the constituents minus
the square of the binding energy. We assume
$m_{s}\gtrsim\frac{\lambda\Lambda_{6}}{2}\approx\Lambda_{D}$ and $\Lambda_{D}$
can be viewed as the binding energy per scalar parton.

The volume of these Q-balls is given by

$V_{\chi}=\frac{Q}{\sqrt{m_{s}^{2}-\frac{\lambda^{2}\Lambda_{6}^{2}}{2}}}\frac{2}{\lambda\Lambda_{6}^{2}}\approx\lambda\frac{Q}{\Lambda_{D}^{3}}$
(3)

Thus the volume of these Q-balls can be parametrically smaller than their
fermionic cousins, the nuclei $N$ which scale as
$V_{N}=\frac{Q}{\Lambda_{D}^{3}}$. Finally, the Q-balls are stable only for
large enough $Q$, where the surface tension term can be ignored, this gives
Grabowska  _et al._ (2018)

$Q_{\rm min}\approx\frac{1}{\lambda}$ (4)

###  II.1 Synopsis of Q-ball Evolution

We start tracking the dark bath at temperature $T_{D}\gtrsim\Lambda_{D}$ at
$z=1000$, with Q-balls broken into their constituent particles. At early
times, the SM baryon temperature $T_{b}\gg\Lambda_{D}$. The constituents have
low enough charge that they are not in thermal equilibrium. At late times,
they undergo a phase transition and rapidly form small Q-balls which are
assumed not to be in thermal equilibrium with the standard model Kusenko and
Shaposhnikov (1998). These small Q-balls will begin to fuse and form larger
Q-balls provided the inter-particle spacing between Q-balls is larger than the
size of Q-balls i.e. $n_{\rm Q-ball}^{-1}>V_{\rm Q-ball}$. This is satisfied
as long as the phase transition happens later than

$\displaystyle z_{\rm tr,max}$
$\displaystyle=\left(\frac{m_{\chi}}{f_{D}\rho_{\text{DM},0}}\frac{\Lambda_{D}^{3}}{\lambda
Q}\right)^{1/3}$ $\displaystyle\approx 24\left(\frac{\Lambda_{D}}{1\text{
K}}\right)^{4/3}\left(\frac{f_{D}}{0.4\%}\right)^{-1/3}\left(\frac{\lambda}{10^{-7}}\right)^{-1/3}$
(5)

It is important to note that the analogous calculation for fermions rather
than Q-balls (corresponding to $\lambda=1$, as discussed above) forces the
bound states to form much later than cosmic dawn, ruling out the entire
parameter space. Stated differently, depending on their temperature in the
early universe, throughout their evolution fermionic mCPs are either cold
enough that they constitute a degenerate gas (in which case the composite
states considered in this work are not the relevant degrees of freedom), or so
warm that they cannot cool the baryons sufficiently to explain EDGES.

We choose the initial dark temperature of the Q-balls at $z=1000$ to be high
enough that the phase transition happens after $z_{\rm tr,max}$. The
requirement that the first Q-balls formed are (i) small enough to be out of
thermal equilibrium with SM baryons at this redshift and (ii) stable at their
size sets

$\lambda>7\times
10^{-8}\left(\frac{m_{\chi}}{\text{MeV}}\right)^{-1/5}\left(\frac{\Lambda_{D}}{1\text{
K}}\right)^{1/5}\left(\frac{f_{D}}{0.4\%}\right)^{-1/5}$ (6)

As the Q-balls become bigger, they will interact with the SM baryons and start
extracting energy from the standard model, leading to heating of the Q-ball
sector. But the temperature of this sector cannot exceed $\Lambda_{D}$ since
this results in Q-ball fission, leading to loss of thermal contact with the
baryon bath, which leads to Hubble cooling. Once $T_{b}$ drops sufficiently
that it is unable to transfer enough heat to hinder Q-ball growth, there is
rapid fusion resulting in large Q-balls. The maximum size of the Q-ball in
this case is set by Coulomb repulsion—as the Q-ball becomes larger, the
repulsion from the dark $U(1)$ grows and it inhibits Q-ball growth beyond a
certain size. Parameters are chosen so that this phase of Q-ball formation
occurs around the redshifts of interest to the EDGES experiment. At this
stage, the Q-balls scatter with the SM baryons, cooling the SM baryons and
explaining the EDGES observations.

###  II.2 The Size of the Q-ball

It is necessary to limit the size of the Q-balls so that they can coherently
scatter with the baryons as well as provide enough heat capacity in the Q-ball
sector to cool the SM bath. The dark $U(1)$ provides the Coulomb repulsion
necessary to enforce this limit. Since $g_{D}\gg\epsilon_{s}$, we will ignore
the Coulomb repulsion from electromagnetism in this section.

How do the Q-balls form? We follow the prescription developed in Krnjaic and
Sigurdson (2015), to account for the Coulomb repulsion due to the dark $U(1)$.
The ratio of fusion rate to Hubble rate for $\\{Q_{1},Q_{2}\\}$ fusion is

$\frac{n\sigma
v}{H}=\frac{10^{32}f_{D}\lambda^{\frac{2}{3}}z^{\frac{3}{2}}}{Q^{\frac{5}{6}}}\left(\frac{10~{}{\rm
K}}{\Lambda_{D}}\right)^{2}\sqrt{\frac{T_{D}}{10~{}{\rm
K}}\left(\frac{10~{}{\rm K}}{m_{s}}\right)^{3}}\gg 1$

As the Q-balls grow bigger, there is an increased Coulomb barrier to fusion as
treated in Krnjaic and Sigurdson (2015). As the Q-balls become larger, this
cross-section is suppressed by the Coulomb barrier. This is captured by
$P_{G}(T)$ the Gamow factor Gamow (1963) which is the temperature dependent
factor that captures the effects of the Coulomb barrier. This factor is:

$\displaystyle P_{G}(T)=e^{-G_{E}}=e^{-\sqrt{\frac{2\mu
Q_{1}^{2}Q_{2}^{2}g_{D}^{4}}{T_{D}}}}$ (7)

Here $Q_{1}$ and $Q_{2}$ are the sizes of the two Q-balls respectively and
$\mu$ is their reduced mass. Thus fusion freeze-out depends critically on the
Gamow factor. From the Gamow factor, it is clear that large-Q-ball small-Q-
ball fusion will dominate over large-large type fusion due to weaker Coulomb
repulsion. Moreover, as the Q-balls grow in size, the number density of larger
Q-balls is lower than that of smaller Q-balls. Further the cross-section for a
smaller Q-ball to merge with a larger Q-ball is set by the geometric size of
the larger Q-ball. All of these factors imply that the growth of the Q-balls
in our case is dominated by the mergers of small Q-balls with larger Q-balls.

Let us now see how big these Q-balls can get i.e. estimate the freeze out of
the fusion process. As seen in Eqn. II.2, in the absence of the exponentially
suppressed Gamow factor, the rate of the fusion process is very rapid compared
to Hubble. The size of the Q-ball is then restricted purely by the exponential
suppression from the Gamow factor which forces the process to freeze out.

Taking $T_{D}\approx\mu\approx\Lambda_{D}$, for small-large fusion, the Gamow
exponent is

$G_{E}\approx Qg_{D}^{2}$ (8)

This places a bound on fusion growth,

$Q^{\rm lim}_{\rm Gamow}\approx g_{D}^{-2}$ (9)

This limit on the Q-ball size arising from the inhibition of their growth is
the same ball-park as the stability limit Grabowska  _et al._ (2018) $Q^{{\rm
lim}}_{\rm stab}=1/g_{D}^{2}$ that can be placed on their size due to Coulomb
repulsion. Q-ball freeze out occurs only due to the exponential dependence on
Q-ball size in the Gamow factor. Thus, Q-balls whose sizes are close to, but
smaller than the Gamow limit are rapidly formed. This implies that as the
universe expands and the temperatures drop, Q-balls will continue to grow
until the Gamow limit is reached.

It is also important to consider the heat that is released by the fusion
process as the Q-balls grow. Each fusion process releases roughly
$\sim\Lambda_{D}$ per unit $Q$. Thus in the roughly $\sim Q$ fusion processes
that occur to form a Q-ball of size $Q$, approximately $Q\Lambda_{D}$ energy
into $\sim Q$ particles is released. Thus the heat released in the fusion
cannot change the temperature by more than $\Lambda_{D}$ and thus does not
hinder fusion.

While we expect this mechanism to produce a range of Q-ball masses, the charge
to mass ratio of all these Q-balls is the same. As seen later in Sec. IV, for
most of the relevant parameter space, results depend only on the charge to
mass ratio, so it is justified to make a simplifying assumption that all
Q-balls are of the same mass $m_{\chi}$. Note that this analysis of the Q-ball
size is independent of the baryon temperature $T_{b}$. As we show in the
following section, $T_{b}$ is an important parameter in determining the number
of mCPs that are fused into Q-balls but it does not determine the maximum size
of a Q-ball.

###  II.3 Heat Transfer

In order to understand heat transfer with the SM bath, we start by deriving
the transfer cross-section for Q-balls to scatter with baryons. The
differential cross-section for a mCP with charge $\epsilon$ to scatter with
protons/electrons is Dvorkin  _et al._ (2014),

$\frac{d\sigma}{d\cos\theta}=\frac{2\pi\epsilon^{2}\alpha^{2}}{\mu^{2}v_{\rm
rel}^{4}(1-\cos\theta)^{2}}.$ (10)

with $\mu$ the reduced mass and $v_{\rm rel}$ the relative velocity.

The forward divergence is cut off by the Debye mass of the mediator. For the
SM photon, the Debye mass squared is given by

$\Pi_{A}=e^{2}\left(\frac{x_{e}n_{b}}{T_{\rm b}}\right)$ (11)

where $x_{e}\equiv n_{e}/n_{\rm H}$ is the free-electron fraction, determined
using Ali-Haïmoud and Hirata (2011); Lee and Ali-Haïmoud (2020). The Debye
mass is approximately $10^{-6}$ eV at $z=1000$ and $3\times 10^{-8}$ eV at
$z=10$. The Debye mass square of the dark photon is

$\Pi_{A^{\prime}}=g_{D}^{2}\left(\frac{n_{D}}{T_{\rm D}}\right)$ (12)

For the parameter space we are interested in, $g_{D}^{2}n_{D}\ll
e^{2}x_{E}n_{b}$ and $T_{D}\leq T_{b}$, such that
$\Pi_{A^{\prime}}\ll\Pi_{A}$. Hence we take only the SM photon Debye mass to
regulate the divergence. Finally, for elementary charges, $q_{\rm max}=2\mu
v_{\rm rel}$, such that the $\theta$ integral is taken between the limits
$\theta=\\{-1,\frac{2\epsilon\alpha\sqrt{\Pi_{A}}}{3T_{b}}\\}$. For Q-balls,
$q_{\rm max}\sim\text{Min}\left(R_{\rm Q-ball}^{-1},2\mu v_{\rm rel}\right)$,
such that $\theta_{\min}=1-\frac{q_{\rm max}^{2}}{2\mu^{2}v_{\rm rel}^{2}}$.
The thermal-averaged transfer cross-section in the $q^{2}_{\rm max}\gg\Pi_{A}$
limit is given by integrating Eqn. 10 over $\theta$, giving,

$\sigma_{T}=\frac{2\pi\epsilon^{2}\alpha^{2}\xi}{\mu^{2}v_{\rm rel}^{4}}$ (13)

with
$\xi=\ln\left({\frac{9T_{b}^{3}}{4\pi\epsilon^{2}\alpha^{3}x_{e}n_{b}}}\frac{q_{\rm
max}^{2}}{2\mu^{2}v_{\rm rel}^{2}}\right)$. In the region of interest, it is
safe to ignore the factor $\frac{q_{\rm max}^{2}}{2\mu^{2}v_{\rm rel}^{2}}$
since it is inside the log.

Next, we compare the rate for Q-balls scattering off baryons to the Hubble
rate:

$\frac{n_{b}\sigma_{T}v_{\text{rel}}}{H}\sim
10^{-18}\left(\frac{m_{\chi}}{\text{MeV}}\right)^{-1/2}\left(\frac{T_{D}}{10\text{
K}}\right)^{-3/2}\left(\frac{\epsilon_{\chi}}{10^{-14}}\right)^{2}z^{3/2}$
(14)

As a result, the smallest Q-balls with charges $\sim\epsilon_{s}$ which we
take to obey stellar-cooling constraints discussed next in Eqn. 20, are never
in thermal contact with the SM. We also see that larger Q-balls with charge
$\epsilon_{\chi}\gtrsim 10^{-7}$ can interact with the SM bath. At
temperatures around $T_{D}\approx\Lambda_{D}$, both small Q-balls and large
Q-balls can co-exist. Defining $\mathcal{F}_{\rm Q-ball}(z)$ as the mass
fraction of large Q-balls i.e. Q-balls with $Q\sim Q^{\rm lim}_{\rm Gamow}$,

$\displaystyle\dot{T}_{D}^{\rm ref}(\mathcal{F}_{\rm Q-ball})=$
$\displaystyle-2HT_{D}+\frac{2}{3}\frac{m_{\chi}x_{e}\rho_{b}}{\left(m_{\chi}+m_{b}\right)^{2}}\mathcal{F}_{\rm
Q-ball}\frac{\sigma_{0}}{u_{\chi,b}^{3}}$
$\displaystyle\times\left\\{\sqrt{\frac{2}{\pi}}\left(T_{b}-T_{D}\right)\right\\}$
(15)

Here $u_{\chi,b}=\sqrt{\frac{T_{D}}{m_{\chi}}+\frac{T_{b}}{m_{b}}}$ is the
average relative velocity due to thermal motion and
$\sigma_{0}=\sigma_{T}v_{\rm rel}^{4}$. We have verified that the bulk
relative velocity between the $\chi$ bath and SM fluids does not contribute
substantially to the thermal evolution of either fluid. When
$T_{D}\gtrsim\Lambda_{D}$, the relevant degrees of freedom are the elementary
charges, which have no thermal contact with the SM such that the dark fluid
cools due to Hubble expansion. When $T_{D}$ drops below $\Lambda_{D}$ there is
rapid Q-ball formation. These Q-balls can now interact with the SM and heat
up, but the temperature cannot exceed $\Lambda_{D}$; after all, thermal
contact with the SM would immediately be lost. Consequently if the second term
in Eqn. 15 dominates for $\mathcal{F}_{\rm Q-ball}\rightarrow 1$, then
$\mathcal{F}_{\rm Q-ball}$ adjusts to smaller values so as to keep
$\dot{T}_{D}=0$. Thus, we set [Note: eqn below changed since initial dark temp
is greater than Lambda, not zero.]

$\dot{T}_{D}=\begin{cases}\dot{T}_{D}^{\rm ref}(\mathcal{F}_{\rm
Q-ball}=0)\quad&T_{D}>\Lambda_{D}\\\ \max\left(0,\dot{T}_{D}^{\rm
ref}(\mathcal{F}_{\rm Q-ball}=0)\right)&T_{D}\leq\Lambda_{D}\end{cases}$ (16)

In the regime where $\dot{T}_{D}=0$, we can solve for the $z$ dependent
fraction in Q-balls $\mathcal{F}_{\rm Q-ball}$ by setting Eqn. 15 to 0. We
find for $T_{D}\leq\Lambda_{D}$,

$\displaystyle\mathcal{F}_{\rm Q-ball}=$
$\displaystyle\textrm{Min}\left(1,2HT_{D}\times\left[\frac{2}{3}\frac{m_{\chi}x_{e}\rho_{b}}{\left(m_{\chi}+m_{b}\right)^{2}}\frac{\sigma_{0}}{u_{\chi,b}^{3}}\right.\right.$
$\displaystyle\left.\left.\times\left\\{\sqrt{\frac{2}{\pi}}\left(T_{b}-T_{D}\right)\right\\}\right]^{-1}\right)$
(17)

We can see that even after the elementary charges have cooled to $\Lambda_{D}$
and formed Q-balls, in the limit where $T_{b}\gg T_{D}$ and when interactions
are strong enough, the quantity in square brackets is much larger than Hubble
cooling and hence $\mathcal{F}_{\rm Q-ball}\rightarrow 0$. This happens
because in this limit, large Q-balls that form immediately break up into
smaller ones that are not in thermal contact with the SM. As the disparity
between $T_{b}$ and $T_{D}$ shrinks, $\mathcal{F}_{\rm Q-ball}\rightarrow 1$.

The time evolution of the baryon temperature obeys

$\displaystyle\dot{T}_{\rm b}=$
$\displaystyle-2HT_{b}+\frac{2}{3}\frac{m_{b}x_{e}\rho_{D}}{\left(m_{\chi}+m_{b}\right)^{2}}\frac{\mathcal{F}_{\rm
Q-ball}f_{D}}{1+f_{\rm He}+x_{e}}\frac{\sigma_{0}}{u_{\chi,b}^{3}}$
$\displaystyle\times\left\\{\sqrt{\frac{2}{\pi}}\left(T_{D}-T_{b}\right)\right\\}+\Gamma_{C}(T_{\rm
CMB}-T_{b})$ (18)

where $f_{\rm He}\equiv n_{\rm He}/n_{\rm H}$ is the helium fraction and
$\Gamma_{C}$ is the Compton scattering rate.

While the initial temperatures of the CMB and baryon gas are set by
observation, the initial dark temperature has to be chosen such that the
Q-ball phase transition occurs no earlier than $z_{\rm tr,max}$. As described
in Section II.2, we find this value by incorporating Eqn. 6 into Eqn. 5, to
ensure that the Q-balls are stable and out of thermal equilibrium with the
baryonic gas when they first form. The initial dark temperature must also
exceed $2\Lambda_{D}$ so that we always start off with Q-balls completely
broken apart into their constituents. Combining these requirements yields the
initial conditions

$\displaystyle T_{b}(z=1000)$ $\displaystyle=T_{\rm CMB}(z=1000)\approx T_{\rm
CMB}^{0}\times 1000$ $\displaystyle T_{\rm CMB}^{0}$
$\displaystyle=2.725~{}\textrm{K}$ $\displaystyle T_{D}(z=1000)$
$\displaystyle=\max\left(2\Lambda_{D},\left(\frac{1000}{z_{\rm
tr,max}}\right)^{2}\Lambda_{D}\right)$
$\displaystyle\approx\Lambda_{D}\max\left(2,1100\left(\frac{m_{\chi}}{\text{MeV}}\right)^{-2/15}\right.$
$\displaystyle\hphantom{\Lambda_{D}\max\left(2,1100\right.}\left.\left(\frac{\Lambda_{D}}{1\text{
K}}\right)^{-38/15}\left(\frac{f_{D}}{0.4\%}\right)^{8/15}\right)$

##  III Existing Limits

As alluded to in the introduction, the constraints on composite mCPs can be
quite different from elementary mCPs of the same charge. We elucidate further
below.

Stellar bounds: For $\Lambda_{D}\ll 1\,\textrm{keV}$, the relevant degrees of
freedom in the interior of stars and supernovae are the elementary mCPs, and
their charge is restricted to $\epsilon_{s}<10^{-14}$ for small enough
$m_{s}$. The Q-balls are never produced in stellar environments. However the
limit on the elementary charges translates to a limit on Q-ball charge:

$\epsilon_{\chi}<10^{-14}\frac{m_{\chi}}{\Lambda_{D}}.$ (20)

BBN and CMB  $\mathbf{N_{\rm\bf eff}:}$ As we have seen in the previous
section, when there is significant thermal contact with baryons and
$T_{b}\gg\Lambda_{D}$, $\mathcal{F}_{\rm Q-ball}\rightarrow 0$ and the
relevant degrees of freedom are the elementary charges before recombination.
Thermal equilibrium with the SM is reached only if Muñoz and Loeb (2018);
Creque-Sarbinowski  _et al._ (2019) $\epsilon_{s}\gtrsim
10^{-8}\left(\frac{m_{\chi}}{10~{}{\rm K}}\right)^{\frac{1}{2}}$. This is more
restrictive than stellar constraints only when $m_{s}\approx\Lambda_{D}\leq
1~{}\mu\textrm{eV}$. Dark photons arising from bremsstrahlung and mesons from
dark fusion are produced at the temperature of the dark bath and hence do not
contribute appreciably to $N_{\rm eff}$ either.

CMB power spectrum: The effect of mCP scattering on protons was investigated
in Kovetz  _et al._ (2018), and constraints from Planck 2015 data effectively
ruled out mCPs as a solution to EDGES for $f_{D}>0.4\%$. It is interesting to
note that since these limits only depend on the charge to mass ratio
$\frac{\epsilon_{\chi}}{m_{\chi}}$, they apply equally to Q-balls as well as
elementary charges. However, it was found in Kovetz  _et al._ (2018) that no
limits exist for $f_{D}\leq 0.4\%$, so we restrict ourselves to smaller
fractions.

##  IV Results

We now display results obtained by numerically solving the coupled
differential equations for time evolution. We consider a benchmark Q-ball mass
$m_{\chi}=1$ MeV, and charge $\epsilon_{\chi}=4\times 10^{-6}$ and
$f_{D}=0.4\%$. We start by tracking $\mathcal{F}_{\rm Q-ball}(z)$ for
different $\Lambda_{D}$ in Fig. 1. For large $z$, the Q-balls are broken apart
into their constituents so $\mathcal{F}_{\rm Q-ball}=0$. For lower $z$,
cooling due to Hubble expansion results in a phase transition and the
subsequent formation of large Q-balls, causing an increase in
$\mathcal{F}_{\rm Q-ball}$. We see that for greater $\Lambda_{D}$, the
increase in the fraction of the dark bath in Q-balls happens earlier as it is
easier for Q-balls to form.

Figure 1: The evolution of the fraction of the dark millicharged bath in large
Q-balls

as a function of redshift is shown for different $\Lambda_{D}$, the dark
confining scale. Smaller $\Lambda_{D}$ leads to smaller Q-ball fractions.

In Fig. 2, the time evolution of the baryonic temperature $T_{b}$ and the dark
temperature $T_{D}$ are shown for different choices of $\Lambda_{D}$, the dark
confining scale. The CMB temperature $T_{\rm CMB}$ and the baryon temperature
$T_{b}$ in the absence of interacting DM are shown in black for reference.

The solid and dashed colored lines track the dark and baryonic temperatures
for different $\Lambda_{D}$. We see that as per Eqn. 5, models with smaller
$\Lambda_{D}$ are forced to reach the phase transition temperature at later
times, and correspondingly must have a greater dark temperature at $z=1000$,
as in Eqn. LABEL:init-cond. For the choice of parameters in the figure, this
condition is restrictive only when $\Lambda_{D}\lesssim 10\text{ K}$, in which
case the initial dark temperature is too large for there to be much cooling of
SM baryons. The requirement that the Q-balls are initially broken apart
prevents models with $\Lambda_{D}\gtrsim 1000\text{ K}$ from being effective
at cooling, through a similar increase in the initial dark temperature.

Figure 2: Temperature evolution of the baryonic and DM bath are plotted as a
function of redshift $z$. The CMB temperature and the baryon temperature
without DM are plotted in black. The solid lines track the evolution of the
dark temperature $T_{D}$ for different $\Lambda_{D}$, the dark confining
temperature. The dashed lines track the baryon temperature $T_{b}$ for
different $\Lambda_{D}$ with the same color code as $T_{D}$. The error bar
marks the baryonic temperature at $z=17$ as measured by the EDGES
collaboration. For $\Lambda_{D}\gtrsim 10\text{ K}$ the initial dark
temperature decreases with lower $\Lambda_{D}$, but for low enough
$\Lambda_{D}$ the requirement from Eqn. LABEL:init-cond becomes restrictive so
the initial dark temperature begins to increase with lower $\Lambda_{D}$.

Figure 3: The contours that explain the EDGES anomaly in the Q-ball charge
$\epsilon_{\chi}$ vs Q-ball mass $m_{\chi}$ plane are shown for different
choices of the confining scale $\Lambda_{D}$ for mCP bath fractions of
$f_{D}=0.04\%$ (left) and $f_{D}=0.4\%$ (right). Also shown are stellar
cooling constraints from Eqn. 20. The elementary charge solution from Muñoz
and Loeb (2018) is shown in black.

We next discuss the contours that explain EDGES in the $\epsilon_{\chi}$ vs
$m_{\chi}$ plane and compare it to the parameter space derived for elementary
charges in Muñoz and Loeb (2018). Given a DM fraction, elementary charges that
explain EDGES obey $\epsilon_{\rm elem}\propto m_{\rm elem}$ (as seen with the
black curve). This happens due to the following reason: for a fixed DM
fraction, a drop in $T_{b}$, $\Delta T_{b}$ is associated with an increase in
dark temperature $\Delta T_{D}\propto m_{\rm elem}\times\Delta T_{b}$, i.e.
larger elementary masses $m_{\rm elem}$ undergo larger temperature gain
because of equipartition. Another way to see this is that the total energy
gained is equal to $n_{\rm elem}\times\Delta T_{D}$ and the number density is
inversely proportional to $m_{\rm elem}$, and hence $T_{D}$ is directly
proportional to $m_{\rm elem}$. Starting with an initially-cold dark bath
$T_{D}\ll T_{b}$, the proportionality factor ensures that $T_{D}\propto m_{\rm
elem}$ throughout. This in turn implies that the elementary charges’ thermal
velocity is independent of the mCP mass. Finally, the heat transfer is
proportional to the transfer cross-section given in Eqn. 13, which is
dependent only on the charge to mass ratio since the velocity is mass-
independent. Thus, this behavior applies to very small masses. It was also
pointed out in Muñoz and Loeb (2018) that for a choice of DM fraction, there
is also a maximum mass due to the same equipartition arguments, $m_{\rm
elem}\leq\mu_{b}f_{\rm elem}\Omega_{c}/\Omega_{b}$. The elementary charge
required to explain EDGES obeys Muñoz and Loeb (2018),

$\epsilon_{\rm elem}\approx 6\times 10^{-7}\left(\frac{m_{\rm elem}}{\rm
MeV}\right)\left(\frac{10^{-2}}{f_{\rm elem}}\right)^{\frac{3}{4}}.$ (21)

It is important to note that the entirety of the elementary charge solution is
ruled out Creque-Sarbinowski  _et al._ (2019).

Next we discuss the contours for Q-balls with different $\Lambda_{D}$, shown
in Fig. 3. In each case, we mark out the unphysical region where the
constituent elementary charges are ruled out by stellar constraints from Eqn.
20. For the same reason as explained for the elementary charge solution, we
observe a linear relationship $\epsilon_{\chi}\propto m_{\chi}$ for Q-balls in
the intermediate mass range (most clearly visible for $\Lambda_{D}=10\text{
K}$). However, at lower mass there is deviation from this behavior as the
requirement from Eqn. 5 becomes restrictive and forces a large initial dark
temperature. This means that a greater millicharge is required for cooling
than would be expected from the $\epsilon_{\chi}\propto m_{\chi}$ relation.
Once again, there exists a cut-off mass that is now $\Lambda_{D}$ dependent.
For smaller $\Lambda_{D}$, the mCP bath stays elementary for longer, i.e.
$\mathcal{F}_{\rm Q-ball}=0$ for longer. To compensate, a smaller Q-ball mass
$m_{\chi}$ is required to increase heat capacity, so as to reach temperatures
below $\Lambda_{D}$ sufficiently soon. As a corollary, larger $\Lambda_{D}$
results in an enhanced range in mass where the EDGES solution is viable.
However, larger $\Lambda_{D}$ translates to stricter stellar constraints and
for large enough $\Lambda_{D}$, the charge required to explain EDGES is ruled
out by Eqn. 20. The Q-balls in most of the parameter space shown in Fig. 3 do
not survive galaxy formation. The parameter space for which Q-balls do not
break up in the galaxy is given in Eqn. 23 and can be recast as,

$m_{\chi}^{\rm gal}\lesssim
86~{}\textrm{eV}\left(\frac{\Lambda_{D}}{1~{}\textrm{K}}\right)$ (22)

Thus even for $\Lambda_{D}\approx 10~{}\textrm{K}$, the Q-balls resize
themselves to masses below $1~{}\textrm{keV}$, making prospects for direct
detection tricky.

##  V Conclusion

Making mCDM inherently composite is a simple nuance with parallels in SM
baryons. In this work, we have considered this possibility and explored its
myriad consequences with specific emphasis on explaining the EDGES anomaly.
The DM degrees of freedom are Q-balls at temperatures below the confining
scale and elementary charges at temperatures above it. For an appropriately
chosen confining scale $\Lambda_{D}$, the elementary charges are the degrees
of freedom during BBN, CMB and in the interior of stars. The elementary
charges are chosen to be feeble enough to evade all these constraints.
However, at temperatures below the confining scale, these rapidly fuse into
Q-balls increasing in size till they reach a size determined by stability
considerations due to repulsion. These Q-balls now have large enough charges
that coherently scatter with baryons at temperatures around $z=17$, relevant
for physics during the dark ages, without suffering from the strict stellar
and cosmology constraints that apply to elementary mCPs. Thus, we find a large
unconstrained parameter space for mCP Q-balls for $f_{D}\leq 0.4\%$, that
explains the EDGES anomaly. In the next few years, this signal will also be
accessible to a slew of experiments sensitive to the global 21 cm signal such
as SARAS2 Singh  _et al._ (2018), LEDA Price  _et al._ (2018), SCI-HI/PRIZM
Voytek  _et al._ (2014), HYPERION Presley  _et al._ (2015) and CTP Nhan  _et
al._ (2017). We also find that there is a novel dark phase, where the dark
bath can exist as an admixture of elementary charges that do not interact with
baryons and composite Q-balls that do, with the fraction in each adjusting so
as to balance the heat transfer from baryons with cooling due to Hubble
expansion. This keeps the dark bath at a constant temperature until the
baryons become cool enough that Hubble cooling dominates heat transfer from
baryons. While this phase was an intriguing curiosity in this work, in the
early universe this can have interesting consequences to thermal freeze-out of
mCPs with a confining force. Finally, the avoidance of stellar and cosmology
constraints due to the composite nature of the DM Q-balls provides a vastly
larger parameter space that is unconstrained compared to elementary mCPs that
do not confine. It is interesting to ask if these Q-balls can be probed in
terrestrial experiments. This task is made more difficult by the fact that
galaxy formation has the potential to destabilize the Q-balls. Galaxy
formation results in DM gaining virial velocities velocities $v_{\rm
vir}\approx 10^{-3}$. Self-interactions are large enough to break up the
Q-balls once more if the kinetic energy exceeds the confining scale. Thus the
Q-balls stay intact till today only if

$m_{\chi}v_{\rm vir}^{2}\lesssim\Lambda_{D}.$ (23)

Thus for large enough Q-ball masses $m_{\chi}$, there is significant fission
in galaxies, the Q-balls are resized into smaller ones that obey Eqn. (23)
which are present in the galaxy today. These smaller Q-balls should
nevertheless be present in the galaxy today since the dark photon sets the
range for self-interactions Lasenby (2020) and cuts off long-range galactic
processes such as evacuation from the galactic disk Chuzhoy and Kolb (2009);
Dunsky  _et al._ (2019) and retention in galactic magnetic fields Harnik  _et
al._ (2020), and prevents the mCP from being blown away by the solar wind
Dunsky  _et al._ (2019); Emken  _et al._ (2019).

Independent of its implications for the EDGES anomaly, this parameter space
increases the scope of direct detection experiments sensitive to masses lower
than $1$ MeV, albeit at momentum transfers smaller than $R_{\rm Q-ball}^{-1}$
to retain coherence. Experiments such as SENSEI Barak  _et al._ (2020), DAMIC
Aguilar-Arevalo  _et al._ (2020), super-CDMS Agnese  _et al._ (2018), and even
future proposals Griffin  _et al._ (2018); Essig  _et al._ (2019); Bunting
_et al._ (2017); Chen  _et al._ (2020) are not sensitive to momentum transfers
$q\leq\Lambda_{D}\approx\textrm{meV}$. Instead, manipulation with electric and
magnetic fields Berlin  _et al._ (2020) is a promising detection strategy. For
large enough Q-ball charge, terrestrial accumulation and subsequent detection
Pospelov and Ramani (2020) might be a viable avenue.

## Acknowledgments

We thank Asher Berlin, Diego Redigolo, Gordan Krnjaic and Hongwan Liu for
useful discussions. SR is supported in part by the NSF under grant
PHY-1818899. SR is also supported by the DoE under a QuantISED grant for MAGIS
and the SQMS quantum center. HR is supported in part by NSF Grant PHY-1720397
and the Gordon and Betty Moore Foundation Grant GBMF7946.

## References

  * Bowman  _et al._ (2018) Judd D. Bowman, Alan E. E. Rogers, Raul A. Monsalve, Thomas J. Mozdzen, and Nivedita Mahesh, “An absorption profile centred at 78 megahertz in the sky-averaged spectrum,” Nature 555, 67–70 (2018), arXiv:1810.05912 [astro-ph.CO] .
  * Barkana (2018) Rennan Barkana, “Possible interaction between baryons and dark-matter particles revealed by the first stars,” Nature 555, 71–74 (2018), arXiv:1803.06698 [astro-ph.CO] .
  * Barkana  _et al._ (2018) Rennan Barkana, Nadav Joseph Outmezguine, Diego Redigolo, and Tomer Volansky, “Strong constraints on light dark matter interpretation of the EDGES signal,” Phys. Rev. D 98, 103005 (2018), arXiv:1803.03091 [hep-ph] .
  * Kovetz  _et al._ (2018) Ely D. Kovetz, Vivian Poulin, Vera Gluscevic, Kimberly K. Boddy, Rennan Barkana, and Marc Kamionkowski, “Tighter limits on dark matter explanations of the anomalous EDGES 21 cm signal,” Phys. Rev. D 98, 103529 (2018), arXiv:1807.11482 [astro-ph.CO] .
  * Berlin  _et al._ (2018) Asher Berlin, Dan Hooper, Gordan Krnjaic, and Samuel D. McDermott, “Severely Constraining Dark Matter Interpretations of the 21-cm Anomaly,” Phys. Rev. Lett. 121, 011102 (2018), arXiv:1803.02804 [hep-ph] .
  * Creque-Sarbinowski  _et al._ (2019) Cyril Creque-Sarbinowski, Lingyuan Ji, Ely D. Kovetz, and Marc Kamionkowski, “Direct millicharged dark matter cannot explain the EDGES signal,” Phys. Rev. D 100, 023528 (2019), arXiv:1903.09154 [astro-ph.CO] .
  * Feng and Holder (2018) Chang Feng and Gilbert Holder, “Enhanced global signal of neutral hydrogen due to excess radiation at cosmic dawn,” Astrophys. J. Lett. 858, L17 (2018), arXiv:1802.07432 [astro-ph.CO] .
  * Ewall-Wice  _et al._ (2018) A. Ewall-Wice, T. C. Chang, J. Lazio, O. Dore, M. Seiffert, and R. A. Monsalve, “Modeling the Radio Background from the First Black Holes at Cosmic Dawn: Implications for the 21 cm Absorption Amplitude,” Astrophys. J. 868, 63 (2018), arXiv:1803.01815 [astro-ph.CO] .
  * Pospelov  _et al._ (2018) Maxim Pospelov, Josef Pradler, Joshua T. Ruderman, and Alfredo Urbano, “Room for New Physics in the Rayleigh-Jeans Tail of the Cosmic Microwave Background,” Phys. Rev. Lett. 121, 031103 (2018), arXiv:1803.07048 [hep-ph] .
  * Fialkov and Barkana (2019) Anastasia Fialkov and Rennan Barkana, “Signature of Excess Radio Background in the 21-cm Global Signal and Power Spectrum,” Mon. Not. Roy. Astron. Soc. 486, 1763–1773 (2019), arXiv:1902.02438 [astro-ph.CO] .
  * Liu  _et al._ (2019) Hongwan Liu, Nadav Joseph Outmezguine, Diego Redigolo, and Tomer Volansky, “Reviving Millicharged Dark Matter for 21-cm Cosmology,” Phys. Rev. D 100, 123011 (2019), arXiv:1908.06986 [hep-ph] .
  * Grabowska  _et al._ (2018) Dorota M. Grabowska, Tom Melia, and Surjeet Rajendran, “Detecting Dark Blobs,” Phys. Rev. D 98, 115020 (2018), arXiv:1807.03788 [hep-ph] .
  * Kusenko and Shaposhnikov (1998) Alexander Kusenko and Mikhail Shaposhnikov, “Supersymmetric Q-balls as dark matter,” Phys. Lett. B 418, 46–54 (1998), arXiv:9709492 [hep-ph] .
  * Krnjaic and Sigurdson (2015) Gordan Krnjaic and Kris Sigurdson, “Big Bang Darkleosynthesis,” Phys. Lett. B 751, 464–468 (2015), arXiv:1406.1171 [hep-ph] .
  * Gamow (1963) George Gamow,  _The Quantum Theory of the Atomic Nucleus_ (US Atomic Energy Commission, Division of Technical Information Extension, 1963).
  * Dvorkin  _et al._ (2014) Cora Dvorkin, Kfir Blum, and Marc Kamionkowski, “Constraining Dark Matter-Baryon Scattering with Linear Cosmology,” Phys. Rev. D 89, 023519 (2014), arXiv:1311.2937 [astro-ph.CO] .
  * Ali-Haïmoud and Hirata (2011) Yacine Ali-Haïmoud and Christopher M. Hirata, “Hyrec: A fast and highly accurate primordial hydrogen and helium recombination code,” Phys. Rev. D 83, 043513 (2011), arXiv:1011.3758 [astro-ph.CO] .
  * Lee and Ali-Haïmoud (2020) Nanoom Lee and Yacine Ali-Haïmoud, “Hyrec-2: A highly accurate sub-millisecond recombination code,” Phys. Rev. D 102, 083517 (2020), arXiv:2007.14114 [astro-ph.CO] .
  * Muñoz and Loeb (2018) Julian B. Muñoz and Abraham Loeb, “A small amount of mini-charged dark matter could cool the baryons in the early Universe,” Nature 557, 684 (2018), arXiv:1802.10094 [astro-ph.CO] .
  * Singh  _et al._ (2018) Saurabh Singh _et al._ , “SARAS 2 constraints on global 21-cm signals from the Epoch of Reionization,” Astrophys. J. 858, 54 (2018), arXiv:1711.11281 [astro-ph.CO] .
  * Price  _et al._ (2018) DC Price, LJ Greenhill, Anastasia Fialkov, GIANNI Bernardi, H Garsden, BR Barsdell, J Kocz, MM Anderson, SA Bourke, J Craig, _et al._ , “Design and characterization of the large-aperture experiment to detect the dark age (leda) radiometer systems,” Monthly Notices of the Royal Astronomical Society 478, 4193–4213 (2018).
  * Voytek  _et al._ (2014) Tabitha C. Voytek, Aravind Natarajan, José Miguel Jáuregui García, Jeffrey B. Peterson, and Omar López-Cruz, “Probing the Dark Ages at $z\sim$ 20: The SCI-HI 21 cm All-sky Spectrum Experiment,” Astrophys. J. Lett. 782, L9 (2014), arXiv:1311.0014 [astro-ph.CO] .
  * Presley  _et al._ (2015) Morgan E. Presley, Adrian Liu, and Aaron R. Parsons, “Measuring the Cosmological 21 cm Monopole with an Interferometer,” Astrophys. J. 809, 18 (2015), arXiv:1501.01633 [astro-ph.CO] .
  * Nhan  _et al._ (2017) Bang D. Nhan, Richard F. Bradley, and Jack O. Burns, “A polarimetric approach for constraining the dynamic foreground spectrum for cosmological global 21-cm measurements,” Astrophys. J. 836, 90 (2017), arXiv:1611.06062 [astro-ph.IM] .
  * Lasenby (2020) Robert Lasenby, “Long range dark matter self-interactions and plasma instabilities,” JCAP 11, 034 (2020), arXiv:2007.00667 [hep-ph] .
  * Chuzhoy and Kolb (2009) Leonid Chuzhoy and Edward W. Kolb, “Reopening the window on charged dark matter,” JCAP 07, 014 (2009), arXiv:0809.0436 [astro-ph] .
  * Dunsky  _et al._ (2019) David Dunsky, Lawrence J. Hall, and Keisuke Harigaya, “CHAMP Cosmic Rays,” JCAP 07, 015 (2019), arXiv:1812.11116 [astro-ph.HE] .
  * Harnik  _et al._ (2020) Roni Harnik, Ryan Plestid, Maxim Pospelov, and Harikrishnan Ramani, “Millicharged Cosmic Rays and Low Recoil Detectors,” (2020), arXiv:2010.11190 [hep-ph] .
  * Emken  _et al._ (2019) Timon Emken, Rouven Essig, Chris Kouvaris, and Mukul Sholapurkar, “Direct Detection of Strongly Interacting Sub-GeV Dark Matter via Electron Recoils,” JCAP 09, 070 (2019), arXiv:1905.06348 [hep-ph] .
  * Barak  _et al._ (2020) Liron Barak _et al._ (SENSEI), “SENSEI: Direct-Detection Results on sub-GeV Dark Matter from a New Skipper-CCD,” Phys. Rev. Lett. 125, 171802 (2020), arXiv:2004.11378 [astro-ph.CO] .
  * Aguilar-Arevalo  _et al._ (2020) A. Aguilar-Arevalo _et al._ (DAMIC), “Results on low-mass weakly interacting massive particles from a 11 kg-day target exposure of DAMIC at SNOLAB,” Phys. Rev. Lett. 125, 241803 (2020), arXiv:2007.15622 [astro-ph.CO] .
  * Agnese  _et al._ (2018) R. Agnese _et al._ (SuperCDMS), “Results from the Super Cryogenic Dark Matter Search Experiment at Soudan,” Phys. Rev. Lett. 120, 061802 (2018), arXiv:1708.08869 [hep-ex] .
  * Griffin  _et al._ (2018) Sinead Griffin, Simon Knapen, Tongyan Lin, and Kathryn M. Zurek, “Directional Detection of Light Dark Matter with Polar Materials,” Phys. Rev. D 98, 115034 (2018), arXiv:1807.10291 [hep-ph] .
  * Essig  _et al._ (2019) Rouven Essig, Jesús Pérez-Ríos, Harikrishnan Ramani, and Oren Slone, “Direct Detection of Spin-(In)dependent Nuclear Scattering of Sub-GeV Dark Matter Using Molecular Excitations,” Phys. Rev. Research. 1, 033105 (2019), arXiv:1907.07682 [hep-ph] .
  * Bunting  _et al._ (2017) Philip C. Bunting, Giorgio Gratta, Tom Melia, and Surjeet Rajendran, “Magnetic Bubble Chambers and Sub-GeV Dark Matter Direct Detection,” Phys. Rev. D 95, 095001 (2017), arXiv:1701.06566 [hep-ph] .
  * Chen  _et al._ (2020) Hao Chen, Rupak Mahapatra, Glenn Agnolet, Michael Nippe, Minjie Lu, Philip C. Bunting, Tom Melia, Surjeet Rajendran, Giorgio Gratta, and Jeffrey Long, “Quantum Detection using Magnetic Avalanches in Single-Molecule Magnets,” (2020), arXiv:2002.09409 [physics.ins-det] .
  * Berlin  _et al._ (2020) Asher Berlin, Raffaele Tito D’Agnolo, Sebastian A. R. Ellis, Philip Schuster, and Natalia Toro, “Directly Deflecting Particle Dark Matter,” Phys. Rev. Lett. 124, 011801 (2020), arXiv:1908.06982 [hep-ph] .
  * Pospelov and Ramani (2020) Maxim Pospelov and Harikrishnan Ramani, “Earth-bound Milli-charge Relics,” (2020), arXiv:2012.03957 [hep-ph] .


sec:tracking_prob

# Development of Linear Battery Model for Path Planning with Mixed Integer
Linear Programming: Simulated and Experimental Validation

Drew Scott1, Satyanarayana G. Manyam2, David W. Casbeer3, Manish Kumar4,
Isaac E. Weintraub5, Michael J. Rothenberger6 1Department of Mechanical and
Materials Engineering, University of Cincinnati<EMAIL_ADDRESS>Scientist, Infoscitex corp., a DCS company, Dayton OH, 45431
<EMAIL_ADDRESS>Area Lead, Cooperative & Intelligent Control,
Control Science Center, Air Force Research Laboratory, WPAFB, OH, 45433
<EMAIL_ADDRESS>Professor in Department of Mechanical and Materials
Engineering, University of Cincinnati<EMAIL_ADDRESS>Electronics
Engineer, Control Science Center, Air Force Research Laboratory, WPAFB, OH,
45433<EMAIL_ADDRESS>Mechanical Engineer, Flight Systems
Integration Branch, Air Force Research Laboratory, WPAFB, OH,
45433Distribution Statement A: Approved for Public Release; Distribution is
Unlimited. PA# AFRL-2022-4598This work has been supported in part by AFOSR
LRIR No. 21RQ-COR084.We acknowledge the support from Bryan Brown and Austin
Wessels at the UAV MASTER Lab at University of CincinnatiThis material is
based on research sponsored by the Air Force Research Laboratory and Strategic
Council for Higher Education under agreement FA8650-19-2-9300. The U.S.
Government is authorized to reproduce and distribute reprints for Governmental
purposes notwithstanding any copyright notation thereon. The views and
conclusions contained herein are those of the authors and should not be
interpreted as necessarily representing the official policies or endorsements,
either expressed or implied by the Strategic Council for Higher Education and
the Air Force Research Laboratory (AFRL) or the U.S. Government.

###### Abstract

Mixed Integer Linear Programs (MILPs) are often used in the path planning of
both ground and aerial vehicles. Such a formulation of the path planning
problem requires a linear objective function and constraints, limiting the
fidelity of the the tracking of vehicle states. One such parameter is the
state of charge of the battery used to power the vehicle. Accurate battery
state estimation requires nonlinear differential equations to be solved. This
state estimation is important in path planning to ensure flyable paths,
however when using MILPs to formulate the path planning problem these
nonlinear equations cannot be implemented. Poor accuracy in battery estimation
during the path planning runs the risk of the planned path being feasible by
the estimation model but in reality will deplete the battery to a critical
level. To the end of higher accuracy battery estimation within a MILP, we
present here a simple linear battery model which predicts the change in state-
of-charge (SOC) of a battery given a power draw and duration. This model
accounts for changes in battery voltage due to applied electrical load and
changes in battery SOC. The battery model is presented and then tested against
alternate battery models in numerical and in experimental tests. Further, the
effect the proposed linear model has over a simpler SOC estimation on the
time-to-solve a resource constrained shortest path problem is also evaluated,
implemented in two different algorithms. It is seen that the linear model
performs well in battery state estimation while remaining implementable in a
Linear Program or MILP, with little affect on the time-to-solve.

##  I INTRODUCTION

Path planning and trajectory planning are problems which show up in many
domains, including vehicle routing, Unmanned Aerial Vehicle (UAV) path
planning, and general robotics. These are frequently solved in a discrete
manner, where the environment to be traversed is represented as a set of
discrete nodes and edges for the agent to travel along. In applications where
a battery is used to power the vehicle, the tracking of battery state-of-
charge (SOC) along the path is important to ensure planned paths are feasible
on the real system. If the battery capacity is limited relative to the
expected returned path distance, there lies the risk of planning a path which
exceeds the vehicle’s battery limits but is calculated to be feasible with the
battery model used in the planning. This danger is present in the case of
small UAVs (sUAV), hybrid-fuel UAVs which combine battery power with another
power source, such as a fuel cell of combustion engine, and electric or
hybrid-fuel ground vehicles. Tracking battery SOC along a path in a Mixed
Integer Linear Program (MILP) requires a linear function for the battery,
however battery SOC in reality changes in a nonlinear manner. To track with a
higher fidelity nonlinear battery model, the problem ceases to be linear and
thus must be solved with more costly Mixed Integer Nonlinear Programming
techniques. However, the cost of increased problem complexity may not be
outweighed by the improved battery tracking, depending on the exact problem at
hand. Multi-agent routing, in a MILP form, is often very computationally
expensive and quickly becomes infeasible to solve as the problem sizes
increase. If transformed to a nonlinear programming problem, to the end of
battery tracking accuracy, the time-to-solve scaling is even harsher. We
present here a simple linear battery model which accounts for voltage changes
due to power draw and drop in battery SOC changes. We compare this to several
alternate battery models in both numerical tests and experimental tests in
terms of accuracy of state of charge tracking. Further, this model is
implemented in a Resource Constrained Shortest Path Problem (RCSPP) and the
time-to-solve is compared to an alternate resource constraint where a constant
battery voltage is assumed.

We are motivated by a prior work in [1] where noise-aware path planning is
done for a hybrid-fuel UAV which combines a generator and battery pack as
power sources. Here a simple equation is used to update the battery SOC where
a constant voltage is assumed, so that the change in SOC along an edge is
constant regardless of power draw and current SOC. This is a simplification
made to keep the path planning problem easier, where in reality the change in
battery SOC for a given flight leg depends on many factors, two major ones
being electrical load and the current battery SOC.

The primary contributions of this paper are: i) Numerical and experimental
testing of a linear battery model, on both single-cell and multi-cell
batteries, developed to the end of use in a MILP ii) analysis of the affect
the linear model has on the time-to-solve a RCSPP variant as compared to a
constant voltage assumption. The paper is organized as follows. The battery
tracking problem of concern is presented in Section II. Section III presents
the linear model, along with the alternate models. Section IV presents the
results of a numerical test comparing the linear model with the alternate
models on an 18650 battery cell. Section V presents the results of
experimental tests on a Lithium Polymer (LiPo) battery. An evaluation of time-
to-solve a RCSPP with and without the linear battery model presented in
Section VI. This RCSPP is solved with two algorithms: i) Branch-and-bound ii)
Labeling Algorithm.

###  I-A Prior Work

To the best of the authors’ knowledge, a linear battery model specifically for
use in path planning with MILPs has not been presented and tested in the
existing literature. However, there exist numerous models for battery systems.
These are nonlinear and seek to model the battery behavior as accurately as
possible. There is also much work done on real-time estimation of SOC and
battery health.

Battery models are divided into 3 main types: i) Electro-chemical models; ii)
Equivalent-circuit models; iii) Data-driven models. Electro-chemical models
seek to model the behavior of the battery at a chemical level, whereas the
equivalent-circuit models represent the battery as a circuit system, often a
resistor-capacitor (RC) circuit, with behavior that mimics that of the real
battery. Data-driven models develop a numerical model of the battery system
based on gathered data. Battery state estimation techniques either seek to
determine the state-of-energy, being the remaining energy that can be pulled
from the battery, state-of-charge estimation, power capability, or state-of-
health (SOH). A comprehensive survey of the above estimation and model types
is given in [2].

Tracking battery state in path planning problems has been studied before,
primarily in the case of Electric Vehicles (EVs) and sUAVs where the only
power source is a battery. In [3], a model is presented to check if a given
path is feasible with regard to battery SOC. Here, a path is first produced
and then checked for battery feasibility after the fact in an online fashion.
In [4], the path planning of a UAV is done while accounting for battery SOH.
UAV path planning is considered in [5] where an algorithm is proposed which
includes a sUAV performance model and battery model within the path planning,
which was solved using a particle swarm optimization algorithm.

There is also existing work on Vehicle Routing Problem (VRP) variants applied
to electric vehicles (EV). The Electric VRP (EVRP) is that of finding minimum
global cost routes for a fleet of electric vehicles, often with additional
constraints such as recharging events. An extensive survey of the problems and
solution approaches is given in [6]. There are several studies on hybrid
vehicle routing [7, 8, 9, 10] where a constant charge/discharge is used for
the battery, ignoring the effect battery SOC has on the change in battery SOC
for the same edge. Another study [1] involves a battery/generator hybrid UAV,
where the battery is discharged and charged throughout the flight, where the
battery is also tracked using a constant charge/discharge for a given edge.
These problems, all formulated as MILPS, can be improved to have greater
accuracy in battery modeling using the model presented in this paper without
the significant increase in computational cost when using nonlinear
constraints.

##  II Battery Tracking Problem

Figure 1: Graph with Power Consumption Parameters

As discussed above, the primary challenge in tracking the battery SOC within a
MILP path planning problem arises due to the nonlinear behavior of batteries.
Depending on the chemical principles used for modeling, this nonlinear
behavior can be more pronounced. The path planning, when done on a graph,
returns a path, optimal by some metric, as a series of edges. An example graph
is given in Figure 1, in a case where edges are defined by their power
requirement and time to travel along. Graphs may not be parameterized in this
exact manner, but by metrics which track energy or power usage in some
alternate manner. This path planning is often formulated as a MILP, where the
battery SOC is tracked by a linear equation that is a function of the graph
parameters and decision variables. An example of such is given in Equation
(1), a simplified equation from [1].

$\displaystyle b_{j}\leq
b_{i}+\Delta(b_{i},P_{ij},t_{ij})+M(1-x_{ij})\quad\forall(i,j)\in E$ (1)

where $b_{i}$ is the SOC of the battery at node $i$, $x_{ij}$ is a binary
decision variable for usage of the edge from node $i$ to node $j$, and $E$ is
the set of all edges in the problem. The value $M$ is a large constant such
that if edge $(i,j)$ is not used in the solution then the constraint is
trivially satisfied. $\Delta(b_{i},P_{ij},t_{ij})$ is the function which
determines the change in battery SOC from node $i$ to node $j$ given the power
draw $P$ and time duration $t$ for the given edge. The change in SOC $\Delta$
will be nonlinear. The aim of the linear model is to estimate the $\Delta$
function while still being implementable in a MILP as a linear constraint. The
above equation is for a shortest path variant, whereas a similar constraint
for a VRP or Traveling Salesman Problem variant would require adaptation to
the tour format rather than path format.

##  III Battery Update Equations

The change in SOC with respect to time is found with Equation (2)

$\displaystyle\frac{dS}{dt}=-\frac{I}{C_{m}}=-\frac{I(S,P)}{C_{m}}$ (2)

where $S$ is battery SOC normalized by the maximum charge in Coulombs $C_{m}$
such that at full SOC $S=1$, $t$ is time, and $I$ is the current through the
battery which is positive if the battery is being drained. Assuming there is a
known power draw $P$, the current can be calculated using Equation (3).

$\displaystyle I(S,P)=\frac{P}{V(S,P)}$ (3)

where $V(S,P)$ is a function of battery SOC and power draw.

Substituting yields the change in SOC with respect to time, given in Equation
(4).

$\displaystyle\frac{dS}{dt}=-\frac{P}{V(S,P)C_{m}}$ (4)

Integrating Equation (4) gives a function for $S(t)$, given in Equation (5).

$\displaystyle
S(t)=S_{0}-\int_{t_{0}}^{t}\\!\frac{P}{V(S(\tau),P)C_{m}}\,d\tau$ (5)

where a positive power value $P$ is that which drains the battery. This can be
solved numerically with methods such as Runge-Kutta techniques or similar. The
solution obtained from the numerical methods can be compared to the linear
model estimation, which uses a single time-step in calculation of $S(t)$. For
comparison against the linear model, three alternate battery models are used.
The first is a simple Ohmic drop model and the second is a 1st order RC model,
both from [11]. These are models (2) and (7) respectively as referred to in
the cited study. The Ohmic drop model was chosen due to its simplicity,
comparable to that of the linear model presented here. The 1st order RC model
was chosen as its accuracy was comparable to more complex models presented in
[11], where the higher complexity models saw diminishing returns in accuracy
improvements. The third alternate model used is one which assumes constant
voltage. Here, this is taken to be the nominal voltage of the battery types
testing in Section IV and V.

###  III-A Simple Ohmic Drop Model

The open-circuit voltage (OCV) of a battery can be obtained experimentally and
utilized by a look-up table to find OCV for a given SOC. However, when a load
is applied to the battery, the actual voltage supplied by the battery drops.
This effect, called Ohmic drop, increases with increasing load on the battery.
For higher loads, the battery is supplying a lower voltage and thus, for the
same power draw, a larger current. Therefore, the change in SOC depends on
both the current SOC and the power. The simple Ohmic drop model is defined as:

$\displaystyle V(S,I)=OCV(S)-IR_{0}$ (6) $\displaystyle P=V(S,I)I$ (7)

where $R_{0}$ is the internal resistance of the battery. The current $I$ is
assumed to be a function of power and voltage, as in our case there is a
constant power draw through which the current $I$ is defined, whereas in most
cases the current is assumed to be known directly. Substituting (7) into (6)
gives the following where voltage is now a function of SOC and power draw:

$\displaystyle V(S,P)=\frac{OCV(S)+\sqrt{OCV(S)^{2}-4PR_{0}}}{2}$ (8)

Equation (8) is the simple Ohmic drop model when constant power is used. Here,
$V(S,P)$ will always have a real, positive solution so long as the current is
not so large to cause a negative voltage. Thus, the cases of complex solutions
are ignored.

###  III-B 1st Order RC Model

The first-order model is defined in Equations (9) and (10).

$\displaystyle V_{k}=OCV(SOC_{k})-I_{k}R_{0}-U_{k}$ (9) $\displaystyle
U_{k+1}=\exp(-\tfrac{\Delta t}{\tau})U_{k}+R_{1}[1-\exp(-\tfrac{\Delta
t}{\tau})]\frac{P_{k}}{V_{k}}$ (10)

where $\tau$ is the time constant of the RC circuit which models the battery
behavior, $R_{0}$ is the internal resistance of the battery, $k$ is the
numerical step, $R_{1}$ and $\tau$ are the resistance and time constant of the
simulated circuit, and $U_{k}$ is a secondary variable tracking the hysteresis
effects over time.

###  III-C Nominal Voltage Only

The fourth estimation model used is one which uses the nominal battery voltage
exclusively, and determines change in SOC based off this nominal voltage,
power draw, and time duration of the load. This type of model is used often in
path planning, as discussed in Section I-A. This estimation model is solved
using a single time step as a fair comparison to the linear model, as both can
be implemented in a MILP in this manner. The update equation is given in
Equation (11), where $V_{nom}$ is the nominal battery voltage.

$\displaystyle S(t)=S_{0}-\frac{P(t-t_{0})}{C_{m}V_{nom}}$ (11)

###  III-D Linear Model

The linear model to implement Equation (5) into a MILP constraint is now
presented. A linear fit is made utilizing the simple Ohmic drop in Equation
(6). However, since the change in SOC depends on $\frac{1}{V}$, a linear fit
will be made of the inverse of the voltage given in Equation (6). The linear
model which accounts for the power draw and current SOC is defined in Equation
(12).

$\displaystyle S(t)=S_{0}-P(AS_{0}+BP+C)(t-t_{0})/C_{m}$ (12)

where (A,B,C) are the fit parameters of the linear model. This is a linear fit
of $\frac{1}{V(S,P)}$ where $V(S,P)$ is calculated with Equation (8). Note
that over the time period $(t_{0},t)$, the voltage will change as the battery
is discharged. This is addressed in the RC model and simple Ohmic drop model
via the Runge-Kutta techniques to solve Equation (5). As the linear model is
to be implemented in a MILP, such techniques cannot be used and instead a
single time-step is used, where the voltage, which changes over the time
period, is taken to be constant as the voltage at time $t_{0}$.

Figure 2: OCV Curve - 4S LiPo and 18650 Cell Figure 3: Numerical Test - 18650
cell Figure 4: Experimental Test - 4S LiPo

##  IV Numerical Tests

Here we present the results of a numerical simulation for a single 18650
battery cell, where the constant-power flight legs of the path-planning
problem are simulated. In the path-planning problem, the battery SOC is to be
tracked along a path, where constant power is drawn for some known duration.
The path, in this form, is a queue of these flight legs or edges, and thus can
be represented as a series of constant power draws from the battery for some
duration. An 18650 lithium-ion (Li-Ion) battery with 2500 mAh capacity is
tested, with the OCV-SOC look up table taken from [12]. Plots of the OCV-
discharge curves are given in Figure 2. In each set of tests the 4 estimation
models from Section III are used: i) simple Ohmic drop model in Equation (6)
ii) 1st Order RC model iii) Linear Model presented iv) Nominal Voltage only

The results of the simulation for the 18650 Li-Ion cell are given in Figure 3.
Here, the power pulses are applied such that the battery is nearly depleted by
the end of the simulation, so that the models are tested over nearly the
entire SOC of the battery. It can be seen that all models perform similarly,
with a few percent-SOC difference in the final estimations. Because the power
applied is constant, and thus the current changing very slowly only as a
result of change in battery voltage as the battery is depleted, the transient
effects of the battery are minimal in this application. Cases were the power
changes more frequently may show improved performance of the 1st order RC
model relative to the others as the transient effects will have a greater
influence.

##  V Experimental Tests

The results of using the models to predict battery SOC as compared to real
data from a 4S LiPo battery are presented in this section. A 4S LiPo was used,
with 5000 mAh capacity, and a 50C maximum discharge rate. Measurements for
power, voltage, and current were used with an 1585 RC Benchmark test stand.
The battery load was applied via a motor and propeller, with the motor
controlled by an Electronic Stability Control (ESC) Pulse-Width Modulation
(PWM) setting. In these tests, a constant ESC frequency was applied to achieve
approximately the desired power draw. As seen in Figure 4, the power draw was
not constant, including both the noise and slow reduction in power draw over
time. The drop in power is due to the constant ESC setting rather than power
draw. As the battery discharged and heated up throughout the test, the voltage
decreases, which then reduces power draw for the same ESC PWM setting. From
the battery data, a constant power draw was calculated for each ”pulse”. This
constant power draw, along with time duration of that specific pulse, was
given to each model to predict the battery discharge along that pulse.

The results of the model predictions and the actual battery SOC are given in
Figure 4. It can be seen that the models compare similarly among themselves,
as was the case in the simulation of the 18650 battery in Section IV. Between
the linear model and the actual battery SOC there was a 0.84% SOC difference.

The linear model, due to the limitations of being implemented in a MILP, only
uses a single time step to calculate change in state of charge. When compared
to the prediction model which also uses a single time step but only uses the
nominal voltage of the battery to predict change in state of charge, the
linear model is more accurate while remaining implementable in a MILP for path
planning or mission planning. The other models tested here, described in
Section III, used Runge-Kutta with 100 time-steps to solve Equation (5)
whereas the linear model uses a single time step. However, despite the drastic
difference in time-steps, the linear model is still able to predict battery
SOC with high accuracy for the path planning battery prediction problem in
both the numerical test and on the experimental data.

##  VI MILP Path Planning Test

By design, the linear model is implementable in a MILP for path planning,
however increased complexity in MILP constraints can have a significant effect
on computation time. This effect can also vary between algorithms used to
solve the problem. We implement the linear battery model in a RCSPP [13] and
compare the time to solve with the same problems using nominal battery voltage
only. The MILP formulation is given in Equations (13)-(20).

$\displaystyle\min_{x}\sum_{i}\sum_{j}(C_{ij}x_{ij})$ (13)
$\displaystyle\sum_{j}x_{Sj}=1$ (14) $\displaystyle\sum_{j}x_{jF}=1$ (15)
$\displaystyle\sum_{j}x_{ij}-\sum_{j}x_{ji}=0$ $\displaystyle\forall i\in
N\setminus\\{S,F\\}$ (16) $\displaystyle 0\leq b_{j}\leq B_{m}$
$\displaystyle\forall j\in N\setminus{S}$ (17) $\displaystyle b_{S}=B_{0}$
(18)

where $x_{ij}$ is a binary decision variable which is $1$ if edge from node
$i$to node $j$ is used in the solution and 0 otherwise, $b_{i}$ is the battery
SOC at node $i$, and $S$ and $F$ are the start and goal nodes respectively.
The constraints to track battery state are given as follows:

$\displaystyle b_{j}\leq b_{i}-P_{ij}t_{ij}(A\hskip 1.0ptb_{i}+B\hskip
1.0ptP_{ij}+C)$ $\displaystyle\qquad\quad+M(1-x_{ij})\quad$
$\displaystyle\forall(i,j)\in E$ (19) $\displaystyle b_{j}\leq
b_{i}-\frac{P_{ij}t_{ij}}{V_{nom}}+M(1-x_{ij})\quad$
$\displaystyle\forall(i,j)\in E$ (20)

where $P_{ij}$ is the power draw, $t_{ij}$ is the time to travel along the
edge, $(A,B,C)$ are the linear model parameters, and $V_{nom}$ is the nominal
voltage of the battery. The linear battery model is used as the resource
constraint as implemented with Equation (19). The alternate resource
constraint, using the nominal voltage model, is given in (20).

The effect on running time when using the linear model over the nominal
voltage only model is evaluated when solving with 2 different algorithms: i)
Branch-and-bound ii) Labeling algorithm. The MILP solver in IBM CPLEX v12.9 is
used to solve with brand-and-bound and the labeling algorithm presented in [1]
is modified for this RCSPP problem. Problems are made from 5 nodes up to 100
nodes, with 30 instances generated for each problem size. Nodes are randomly
placed in a 2D space and each node is connected to the 4 closest nodes, by
euclidean distance. Edge costs are set as the Euclidean distance between
nodes. The battery parameters for the problem are defined to simulate the 4S
LiPo battery described in Section V. The averaged results for time-to-solve
with the branch-and-bound is given in Figure 5 and for the labeling algorithm
in Figure 6.

Figure 5: Time-to-Solve with Branch-and-Bound - Battery Model Comparison
Figure 6: Time-to-Solve with Labeling Algorithm- Battery Model Comparison

The time to solve is increased in both the branch-and-bound and the labeling
algorithm when using the linear model. However, this increase in computation
time is more pronounced in the branch-and-bound as compared to the labeling
algorithm. This is expected, as the constraint complexity can have a
significant effect on time-to-solve when using a branch-and-bound method. The
labeling algorithm, which utilizes dynamic programming, has the computation
time primarily affected by the number of feasible, undominated paths and sub-
paths that exist within the problem. While the linear model increases the
battery constraint complexity, it has little effect on the number of feasible,
undominated paths within the problem. Therefore, the battery constraint is
expected to have little effect on the computational time when using the
labeling algorithm. The spikes seen in the time-to-solve with the labeling
algorithm are expected, where infeasible and loosely constrained problems
require many more iterations to complete compared to the average number of
iterations, for the same problem size. A presentation of labeling algorithms
and their use in constrained shortest path problem variants is given in [14].

##  VII Conclusions and Future Work

The linear battery model was presented to the end of accurate battery SOC
estimation for MILP path planning problems concerned with EVs and battery-
powered UAVs. The model was applied to an 18650 battery cell and a 4S LiPo
battery pack where a current was applied to simulate the series of steady-
state flights along graph edges in the path planning problem. Beyond the
accuracy testing, changes to the time-to-solve a RCSPP were also studied. It
was seen that the linear model provides greater accuracy in SOC estimation
over assuming a constant nominal voltage with little effect on computational
time, while remaining a linear constraint able to be applied in a MILP. Future
work includes further testing on the linear model under different
circumstances, and integration with a power or energy consumption model to
accurately predict power usage along an edge. An accurate power model is
important for path planning for the same reason as the battery model, that
being the need to avoid paths which are feasible by the battery and power
models but in actuality are not feasible as the battery will drain to a
critical level before reaching the goal.

## References

  * [1] D. Scott, S. G. Manyam, D. W. Casbeer, M. Kumar, M. J. Rothenberger, and I. E. Weintraub, “Power management for noise aware path planning of hybrid UAVs,” in _American Control Conference (ACC)_ , Atlanta, GA, USA, 2022, pp. 4280–4285.
  * [2] Y. Wang, J. Tian, Z. Sun, L. Wang, R. Xu, M. Li, and Z. Chen, “A comprehensive review of battery modeling and state estimation approaches for advanced battery management systems,” _Renewable and Sustainable Energy Reviews_ , vol. 131, p. 110015, 2020.
  * [3] C. Di Franco and G. Buttazzo, “Coverage path planning for UAVs photogrammetry with energy and resolution constraints,” _Journal of Intelligent & Robotic Systems_, vol. 83, no. 3, pp. 445–462, 2016.
  * [4] R. Schacht-Rodríguez, J.-C. Ponsart, C.-D. García-Beltrán, C.-M. Astorga-Zaragoza, D. Theilliol, and Y. Zhang, “Path planning generation algorithm for a class of UAV multirotor based on state of health of lithium polymer battery,” _Journal of Intelligent & Robotic Systems_, vol. 91, no. 1, pp. 115–131, 2018.
  * [5] A. R. Hovenburg, F. A. de Alcantara Andrade, R. Hann, C. D. Rodin, T. A. Johansen, and R. Storvold, “Long-range path planning using an aircraft performance model for battery-powered sUAS equipped with icing protection system,” _IEEE Journal on Miniaturization for Air and Space Systems_ , vol. 1, no. 2, pp. 76–89, 2020.
  * [6] T. Erdelić and T. Carić, “A survey on the electric vehicle routing problem: variants and solution approaches,” _Journal of Advanced Transportation_ , vol. 2019, pp. 1–48, 05 2019.
  * [7] A. Verma, “Electric vehicle routing problem with time windows, recharging stations and battery swapping stations,” _EURO Journal on Transportation and Logistics_ , vol. 7, no. 4, pp. 415–451, 2018.
  * [8] C. Doppstadt, A. Koberstein, and D. Vigo, “The hybrid electric vehicle–traveling salesman problem,” _European Journal of Operational Research_ , vol. 253, no. 3, pp. 825–842, 2016.
  * [9] F. Y. Vincent, A. P. Redi, Y. A. Hidayat, and O. J. Wibowo, “A simulated annealing heuristic for the hybrid vehicle routing problem,” _Applied Soft Computing_ , vol. 53, pp. 119–132, 2017.
  * [10] G. Hiermann, R. F. Hartl, J. Puchinger, and T. Vidal, “Routing a mix of conventional, plug-in hybrid, and electric vehicles,” _European Journal of Operational Research_ , vol. 272, no. 1, pp. 235–248, 2019.
  * [11] X. Hu, S. Li, and H. Peng, “A comparative study of equivalent circuit models for Li-ion batteries,” _Journal of Power Sources_ , vol. 198, pp. 359–367, 2012.
  * [12] F. Elmahdi, L. Ismail, and M. Noureddine, “Fitting the OCV-SOC relationship of a battery lithium-ion using genetic algorithm method,” in _E3S Web of Conferences_ , vol. 234. EDP Sciences, 2021, p. 00097.
  * [13] H. C. Joksch, “The shortest route problem with constraints,” _Journal of Mathematical analysis and applications_ , vol. 14, no. 2, pp. 191–197, 1966.
  * [14] J. Desrosiers, Y. Dumas, M. M. Solomon, and F. Soumis, “Time constrained routing and scheduling,” _Handbooks in operations research and management science_ , vol. 8, pp. 35–139, 1995.