aainstitutetext: Institute for Particle Physics Phenomenology, Department of
Physics, Durham University,
Durham DH1 3LE, U.K.

# Cosmic microwave background constraints on extended dark matter objects

Djuna Croon a Sergio Sevillano Muñoz<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Primordially formed extended dark objects would accrete baryonic matter and
impact the ionisation history of the Universe. Insisting on consistency with
the anisotropies of the cosmic microwave background, we derive constraints on
the dark matter fraction for various classes of objects, of different sizes.
We introduce a novel scaling technique to speed up numerical calculations and
release our calculation framework in the form of a Mathematica notebook.
Conservatively, we focus on spherical accretion and collisional ionisation. We
find strong constraints limiting the dark matter fraction to subpercent level
for objects of up to $10^{4}$ AU in size.

††preprint: IPPP/24/11

##  1 Introduction

Despite the fact that dark matter (DM) comprises $\sim 85\%$ of the total
amount of matter in our Universe, very little is known about its mass: the
potential mass range of DM candidates spans over 80 orders of magnitude. At
the high end of this spectrum, macroscopic DM candidates, with masses ranging
from comparable to large asteroids to large stars, offer an alternative to
particle- and wave- like DM. Such celestial objects, formed in the early
universe, are so sizeable that they can in principle be detected via their
gravitational effects. As such, important constraints have been set using
gravitational (micro-)lensing techniques (e.g. Alcock  _et al._ (1998a);
Niikura  _et al._ (2019a, b); Green and Kavanagh (2021); Croon  _et al._
(2020a); Bai  _et al._ (2020); Croon  _et al._ (2020b); Crispim Romão and
Croon (2024)), and gravitational waves (e.g. Bird  _et al._ (2016); Ali-
Haïmoud  _et al._ (2017); Kavanagh  _et al._ (2018); Bertone  _et al._ (2020);
Abbott  _et al._ (2019); Chen and Huang (2020); Franciolini (2021); Croon  _et
al._ (2023)) on these objects as DM candidates.

Formed in the early Universe, celestial DM objects are also expected to
accrete throughout their lifetime. Thus, one can study the interaction between
accreting compact objects formed in the early Universe and the cosmic
microwave background (CMB). As the objects accrete matter, they emit high-
energy radiation, which can influence the ionisation and thermal state of the
surrounding medium. This process of energy injection can lead to alterations
in the ionisation fraction, describing the ratio of free electrons to total
hydrogen atoms in the universe. Changes in the ionisation fraction prior to
and during the epoch of recombination have the potential to leave discernible
imprints on the CMB.

Measurements of the CMB by satellites such as Planck have yielded precise data
on its temperature fluctuations and polarisation, and therefore the ionisation
fraction. This has enabled stringent constraints to be placed on the abundance
of primordial black holes (PBHs), popular examples of primordial compact
objects, in the very high mass range $\rm 10\,M_{\odot}<M_{\rm
PBH}<10^{4}\,M_{\odot}$ (see e.g. Carr (1981); Ricotti  _et al._ (2008); Ali-
Haïmoud and Kamionkowski (2017); Poulin  _et al._ (2017)). Ref. Bai  _et al._
(2020) generalised the analysis of PBH to objects with a uniform density
profile comprising 100% of DM. In this work it was found that even relatively
dilute objects of several hundreds of AU could cause observable levels of
ionisation in the CMB.

In this work, we generalise the accretion formalism further to account for
extended dark objects (EDOs) of arbitrary radius and density profiles. In
contrast to Bai  _et al._ (2020), which focused on an analytical calculation
of the ionisation fraction, our computations are primarily numerical which
allows them to be easily adapted to any mass profile. We introduce a scaling
technique that significantly speeds up the numerical calculation. The primary
purpose of this work is to find the constraints on the fraction of DM that the
extended dark objects can comprise. We find stringent constraints for objects
over $10^{6}\,\rm R_{\odot}$ or $10^{4}\rm\,AU$ in the mass range above
$10^{3}\,\rm M_{\odot}$, complementary to those found in microlensing studies
Croon  _et al._ (2020a); Bai  _et al._ (2020); Croon  _et al._ (2020b) and
dominating the constraints found by Graham and Ramani (2023) from the
dynamical heating of stars in ultrafaint dwarfs in the high mass range.

The structure of this paper can be conveniently inferred from the Table of
Contents.

##  2 Accretion

The gravitational influence of extended dark matter objects can lead to the
accretion of matter, leading to a possible alteration of the ionisation
history of the universe. To calculate the actual impact, first we need to
obtain the density profile of accreted matter in the dark matter object, as
well as how this translates into its temperature, ionisation fraction and
electron/proton number. In this section, we will describe how this is
calculated for a generic set of EDO mass functions, closely following the
steps in Refs.Ali-Haïmoud and Kamionkowski (2017); Bai  _et al._ (2020).

###  2.1 Baryonic matter

As a good approximation, we can assume that at the relevant redshifts the
universe is made up of non-relativistic electrons, protons and hydrogen atoms,
which we refer to in this work as baryonic matter. In this way, we can define
the density of baryonic matter as

$\rho=m_{e}n_{e}+m_{p}n_{p}+m_{H}n_{H},$ (1)

where $n_{e}(r,t)$, $n_{p}(r,t)$ and $n_{H}(r,t)$ are the number densities of
electrons, protons and hydrogen atoms, respectively, and $m_{e}\simeq
0.511~{}\text{MeV}$ and $m_{p}\approx m_{H}\simeq 0.938~{}\text{GeV}$ their
corresponding masses. Assuming local charge neutrality implies

$n_{e}=n_{p},$ (2)

which allows us to quantify each parameter through the fraction of ionized
hydrogen atoms, given by

$x_{e}=\frac{n_{e}}{n_{e}+n_{H}}.$ (3)

Therefore, considering that $m_{e}\ll m_{p}$, we obtain

$\displaystyle n_{e}(r,t)=$ $\displaystyle
n_{p}(r,t)\approx\frac{1}{m_{p}}x_{e}(r,t)\rho(r,t),$ (4) $\displaystyle
n_{H}(r,t)=$ $\displaystyle\frac{1}{m_{p}}(1-x_{e}(r,t))\rho(r,t).$ (5)

Having defined these parameters, now we can address the evolution and
thermodynamics of the initially homogeneous baryonic fluid.

####  2.1.1 Evolution and thermodynamics of baryonic matter

The dynamic evolution of the baryon fluid around an EDO111We will calculate
the accretion profile for a single, isolated extended dark matter object. This
approximation will be justified later in section 3, when studying these
objects in a cosmological context. will be dictated by its mass density
$\rho(r,t)$, fluid velocity $\vec{v}(\vec{r},t)$, and internal energy density
$\mathcal{E}(r,t)$ via the set of Navier-Stokes equations

$\displaystyle\dot{\rho}$ $\displaystyle+\frac{1}{r^{2}}(r^{2}\rho
v)^{\prime}=0,$ (6) $\displaystyle\rho$ $\displaystyle\dot{v}+\rho
vv^{\prime}+P^{\prime}=\rho g,$ (7) $\displaystyle\rho$
$\displaystyle\dot{(\mathcal{E}/\rho)}+\rho
v(\mathcal{E}/\rho)^{\prime}+P\frac{1}{r^{2}}(r^{2}v)^{\prime}=\dot{q},$ (8)

which ensures the conservation of mass, energy, and momentum. Herein, we have
used $\dot{}=\partial/\partial t$, ${}^{\prime}=\partial/\partial r$, and have
imposed rotational symmetry such that $\vec{v}(\vec{r},t)=v(r,t)\vec{r}/r$. In
these equations, some terms are still to be defined and specified, which we
proceed to do now. Firstly, in the second equation, the $g(r,t)$ term
corresponds to the radial component of the potential energy sources acting on
the baryonic fluid. In particular, two terms are sourcing this component, the
main one is given by the gravitational pull exerted by the EDO,

$\vec{g}_{\rm grav}(\vec{r})=g_{\rm
grav}(r)\frac{\vec{r}}{r},\qquad\text{with}\qquad g_{\rm
grav}(r)=-G\frac{M(r)}{r^{2}},$ (9)

where $G$ is Newton’s gravitational constant and $M(r)$ is the mass
distribution of the extended dark matter object. In addition to the
gravitational pull, there is also a drag term arising from the interactions
between the accreting and the CMB radiation, but as in Bai  _et al._ (2020) it
can be ignored for the systems we will be studying in this work.

In terms of the internal energy density, $\mathcal{E}(r,t)$, and pressure,
$P(r,t)$, of the baryonic matter, we can use Fermi-Dirac statistics for
electrons and Maxwell-Boltzman statistics for protons and hydrogen, leading to

$\displaystyle\mathcal{E}$
$\displaystyle=\frac{3}{2}T\left[n_{e}f_{\mathcal{E}}(T/m_{e})+n_{p}+n_{H}\right]$
(10) $\displaystyle P$
$\displaystyle=T\left[n_{e}f_{P}(T/m_{e})+n_{p}+n_{H}\right],$ (11)

where $T(r,t)$ is set to be the common temperature of the baryon fluid in
local equilibrium. The function $f_{\mathcal{E}}(X)$ varies from
$f_{\mathcal{E}}(X\ll 1)\sim 1$ to $f_{\mathcal{E}}(X\gg 1)\sim
7\pi^{4}/(270\zeta(3))\simeq 2.10$, and $f_{P}(X)$ from $f_{P}(X\ll 1)\sim 1$
to $f_{P}(X\gg 1)\sim 7\pi^{4}/(540\zeta(3))\simeq 1.05$. These two functions
introduce the relativistic behaviour of the electron, which changes its
thermal properties when $T\gg m_{e}$. Notice that, for a complete analysis, we
would also have to include similar terms for the proton and hydrogen atom, but
throughout this work we will just consider temperatures much smaller than the
mass of the proton, for which the equivalent functions are set to unity. Using
the definitions in Eqs. (4) and (5) we can simplify these two parameters as
follows

$\displaystyle\mathcal{E}$
$\displaystyle=\frac{3}{2}\frac{1}{m_{p}}T\rho\left[1+x_{e}f_{\mathcal{E}}(T/m_{e})\right]$
(12) $\displaystyle P$
$\displaystyle=\frac{1}{m_{p}}T\rho\left[1+x_{e}f_{P}(T/m_{e})\right].$ (13)

The last thing left to define in Eqs. (6-8) is the $\dot{q}(r,t)$ term in the
third equation, which corresponds to the rate of heating (for $\dot{q}>0$) or
cooling (for $\dot{q}<0$) per unit volume. This is given by

$\dot{q}=\frac{4x_{E}\sigma_{\rm T}\rho_{\rm cmb}}{m_{p}m_{e}}(T_{\rm
cmb}-T)\rho,$ (14)

where $\sigma_{\rm T}\simeq 6.65\times 10^{-25}~{}\text{cm}^{2}$ is the
Thomson scattering cross section, $\rho_{\rm cmb}\simeq(1.97\times
10^{-15}~{}\text{eV}^{4})(1+z)^{4}$ is the CMB energy density, and $T_{\rm
cmb}\simeq(2.34\times 10^{-4}~{}\text{eV})(1+z)$ is the CMB temperature at
redshift $z$.

Using the system of equations in Eqs. (6-8) completely defines the evolution
of baryonic matter around the EDO. However, in practice, to handle these
equations we still need to take a set of approximations that will make it
possible to perform this calculation. These are the following:

  1. 1.

We will focus on obtaining a static approximation, for which we can neglect
the time evolution of all the parameters defined so far. In this way, we will
just calculate the solutions for the functions $\rho(r)$, $T(r)$, $P(r)$,
$n_{e}(r)$, $x_{e}(r)$, $\mathcal{E}(r)$ and $v(r)$.

  2. 2.

Taking the hydrostatic approximation, we can take the baryon flow to vanish,
such that $v(r)=0$. This is a valid approximation as long as the gravitational
interaction is not too intense, placing limits on EDOs with a radius of the
same order as their Schwarchild radius.

  3. 3.

We will set $x_{e}(r)$ to be constant when solving for the Navier-Stokes
equations, and will add the effects of the ionisation of the baryon fluid as a
correction to their solution.

  4. 4.

Taking the ratio between the internal energy, $\mathcal{E}(r)$, and pressure,
$P(r)$, we obtain the adiabatic index, $\gamma(r)$,

$\frac{\mathcal{E}(r)}{P(r)}=\gamma(r)-1,$ (15)

which for $x_{e}\ll 1$, gives $\gamma(r)\approx 5/3$, while for $x_{e}\approx
1$ we have $\gamma(r)\approx 5/3$ if $T(r)\ll m_{e}$ and $\gamma(r)\approx
13/9$ if $T(r)\gg m_{e}$. Therefore, we will approximate this function to be
always constant except for $T(r)=m_{e}$, such that

$\gamma(r)=\begin{cases}\frac{5}{3}&\quad T(r)\leq m_{e}\\\ \frac{13}{9}&\quad
T(r)>m_{e}.\end{cases}$ (16)

Doing so, allows us to relate $P(r)$ and the baryon density, $\rho(r)$ as
follows

$P(r)=K\rho(r)^{\gamma},$ (17)

where $K$ is a constant of mass dimension $4-4\gamma$. Inserting this relation
into Eq. (13) leads to

$T(r)=Km_{p}\frac{\rho(r)^{\gamma-1}}{1+x_{e}f_{P}},$ (18)

where in order for $T(r)$ to be continous at $T(r_{\rm rel})=m_{e}$ we will
have different, but constant, values for $K,\gamma$ and $f_{P}$ in each side
of that boundary.

Considering this set of approximations, the Navier-Stokes equations simplify
to

$\frac{GM(r)}{r^{2}}+\gamma K\rho(r)^{\gamma-2}\frac{d\rho(r)}{dr}=0.$ (19)

###  2.2 EDO Mass profiles

The only thing left to define is the mass function, $M(r)$, which will vary
depending on the formation mechanism for the EDO. In this work, we will
consider the following mass profiles:

  * •

Uniform spheres (of constant density): Although this is a toy model, this mass
profile describes quark nuggets, either from standard Witten (1984) or axion
Zhitnitsky (2003) QCD. The mass profile is given by

$M(r)=M\begin{cases}\left(\frac{r}{R}\right)^{3}&\qquad r\leq R\\\ 1&\qquad
R<r.\end{cases}$ (20)

  * •

Boson stars: Obtained by the gravitational collapse of a scalar field Schunck
and Mielke (2003). They do not have an analytical mass function, but they can
be obtained by solving Schrodinger’s equation for a scalar field under the
influence of its gravitational potential Bar  _et al._ (2018). As the solution
converges at infinity, we will define the radius where it encloses $90\%$ of
its mass.

  * •

Navarro-Frenk-White (NFW) subhalos: Empirically obtained density distribution
that closely matches the behaviour of cold dark matter clustering Navarro  _et
al._ (1996). The mass profile is given by

$M(r)=\int_{0}^{r}{\rm d}\hat{r}\,4\pi\hat{r}^{2}\rho_{\rm NFW}(\hat{r}),$
(21)

where

$\rho_{\rm
NFW}(\hat{r})=\frac{\rho_{0}}{\frac{\hat{r}}{R_{s}}\left(1+\frac{\hat{r}}{R_{s}}\right)^{2}},$
(22)

is the NFW density profile, with $R_{s}$ and $\rho_{0}$ being parameters that
define the total subhalo mass and $r_{90}$, the radius that encloses $90\%$ of
its mass. Because the mass of this profile technically diverges, we choose to
cut it off at $100R_{s}$.

  * •

Ultra-compact mini-halos: Although it is widely known that primordial density
fluctuations seed the creation of primordial black hole, this fact depends on
the gravitational pull of the over-density. Thus, for small enough over-
densities that do not form a pbh, it is still possible to seed the growth of
an ultra-compact mini-halo Ricotti and Gould (2009). The mass profile of these
objects is expected to be Bertschinger (1985)

$M(r)=M\begin{cases}\left(\frac{r}{R}\right)^{3/4}&\qquad r\leq R\\\ 1&\qquad
R<r.\end{cases}$ (23)

Figure 1: Mass profiles of the EDOs described in text.

We show the different mass profiles in Figure 1. With this, we are all set to
solve Eq. (19) and see the impact that EDOs can have on the ionisation history
of the universe. In the next subsection, we will proceed to obtain the
density, temperature, ionisation fraction and electron/proton number density
of the accreted baryonic matter for a generic EDO.

###  2.3 Accretion profiles of baryonic matter

As introduced above, we solve our system of equations numerically, which,
although it is less efficient than the analytical route from Ref. Bai  _et
al._ (2020), allows extending our results to any mass function. Moreover, in
Appendix A we present how this computation can be optimised so that solutions
for a specific EDO with a given total mass and radius can be easily
generalised to a different one.

Solving for the density profile of the accreted baryons around the extended
dark matter objects will first require solving Eq. (19) both inside and
outside the EDO, matching both solutions at the boundary. To do so, we also
need to specify the boundary conditions at infinity, defined as

$\displaystyle\lim_{r\to\infty}\rho(r)\equiv\,$ $\displaystyle\rho_{\infty}$
$\displaystyle\lim_{r\to\infty}T(r)\equiv\,$ $\displaystyle T_{\infty}$
$\displaystyle\lim_{r\to\infty}x_{e}(r)\equiv\,$ $\displaystyle\bar{x}_{e},$
(24)

which will be defined by the homogeneous evolution of the Universe, as we will
discuss in section 3. Using Eqs. (13) and (16), and the fact that
$T_{\infty}\ll m_{e}$, we further obtain

$\displaystyle
P_{\infty}=\frac{1}{m_{p}}T_{\infty}\rho_{\infty}\left[1+\bar{x}_{e}\right],\qquad\text{and}\qquad\gamma_{\infty}=5/3.$
(25)

From this, one can find a specific form for $\rho(r)$, which is extended to
the temperature and electron/proton number density using

$\displaystyle T(r)=$ $\displaystyle
Km_{p}\frac{\rho(r)^{\gamma-1}}{1+\bar{x}_{e}}$ (26) $\displaystyle n_{e}(r)=$
$\displaystyle\frac{1}{m_{p}}\bar{x}_{e}\rho(r),$ (27)

where we have defined $K=T_{\infty}/\rho_{\infty}^{\gamma_{\infty}-1}$ from
the boundary condition at infinity. However, this solution is incomplete until
we include the thermodynamic corrections that lead to the cooling due to
interactions with the CMB radiation, the ionisation of the fluid, and the
relativistic transition of electrons. We will now explain how to account for
each of these corrections to our solutions.

####  2.3.1 CMB cooling

One of the first assumptions we made in the previous subsection was setting
the the baryon fluid heating rate to zero (i.e., $\dot{q}=0$). Here we will
allow for different values of this function, and check its impact on the
temperature obtained through Eq. (26).

Starting from the static version for Eqs. (6) and (8),

$\displaystyle\frac{1}{r^{2}}(r^{2}\rho v)^{\prime}$ $\displaystyle=0,$ (28)
$\displaystyle\rho v(\mathcal{E}/\rho)^{\prime}+P\frac{1}{r^{2}}$
$\displaystyle(r^{2}v)^{\prime}=\dot{q},$ (29)

and using Eqs. (13) and (15), we obtain

$v\rho^{\gamma_{\infty}-1}\frac{{\rm d}}{{\rm
d}r}\left(\frac{T}{\rho^{\gamma_{\infty}-1}}\right)=\frac{8\bar{x}_{e}\sigma_{\rm
T}\rho_{\rm cmb}}{3m_{e}(1+\bar{x}_{e})}(T_{\rm cmb}-T),$ (30)

where we have also set the values at infinity for $\gamma(r)$ and $x_{e}(r)$.
To analyze the different limits of this differential equation, let us first
introduce two quantities that derive from the sound speed of the baryonic
fluid at infinity,

$c_{\infty}=\sqrt{\gamma_{\infty}P_{\infty}/\rho_{\infty}}.$ (31)

Defining the Bondi radius Bondi (1952) as the characteristic distance where
the EDO influences the baryon fluid, we obtain

$R_{\rm B}=GM/c_{\infty}^{2},$ (32)

which is related to the escape velocity of the EDO. Similarly, the Bondi time
can be defined as the time it takes a perturbation of the baryon fluid to
cover the Bondi radius, leading to

$t_{\rm B}=GM/c_{\infty}^{3}.$ (33)

With this, we can now introduce these characteristic scales in Eq. (30) to
define the impact of cooling on the EDO formation. Multiplying $t_{\rm B}$ on
both sides, we get

$\frac{v}{c_{\infty}}\rho^{\gamma_{\infty}-1}R_{\rm{B}}\frac{{\rm d}}{{\rm
d}r}\left(\frac{T}{\rho^{\gamma_{\infty}-1}}\right)=\Theta(T_{\rm cmb}-T),$
(34)

where following Ref. Bai  _et al._ (2020) we define the dimensionless cooling
factor as

$\Theta=\frac{8\bar{x}_{e}\sigma_{\rm T}\rho_{\rm
cmb}}{3m_{e}(1+\bar{x}_{e})}t_{\rm B}.$ (35)

For the region of interest, outside of the EDO, we expect both $R_{\rm B}{\rm
d}/{{\rm d}r}(T/\rho^{\gamma_{\infty}-1})/(T/\rho^{\gamma_{\infty}-1})$ and
$v/c_{\infty}$ to be $\mathcal{O}(1)$. Thus, we can see two extreme cases for
this differential equation: for $\Theta\ll 1$, the right-hand side vanishes,
recovering the simplified temperature profile described in Eq. (26). On the
contrary, when $\Theta\gg 1$, the only fixed point for Eq. (34) is at
$T(r)=T_{\rm cmb}$, cooling down the EDO due to the interactions with the CMB.
By continuity, there must be a threshold for which the cooling effect due to
the CMB only partially affects the temperature of the accreted matter. In Ref.
Ali-Haïmoud and Kamionkowski (2017), the authors obtain an estimate for this
boundary of $r_{\rm cool}=\Theta^{-2/3}R_{\rm B}$ using the analytic solution
for PBHs. However, we can arrive at the same conclusion by studying Eq. (34)
itself. Focusing on the cooling factor, we can see that the aforementioned two
extreme cases can be also expressed in terms of a cooling time (defined as the
inverse of the cooling rate), given by

$t_{\rm cool}=\frac{3m_{e}(1+\bar{x}_{e})}{8\bar{x}_{e}\sigma_{\rm T}\rho_{\rm
cmb}}.$ (36)

For $t_{\rm cool}\gg t_{\rm B}$, accreted matter has enough time to escape the
EDO before the CMB cools it down, leading to a negligible effect on its
temperature profile. However, when $t_{\rm cool}\ll t_{\rm B}$, the accreted
matter cools down before it crosses the Bondi radius, which effectively
shrinks the radius at which the fluid can escape. This new effective radius,
$r_{\rm cool}$, can be obtained following the same arguments as for the Bondi
radius, such that

$r_{\rm cool}{\equiv}\frac{GM}{c_{\rm cool}^{2}}=\frac{GMt_{\rm
cool}^{2}}{r_{\rm
cool}^{2}}=\frac{GM}{c_{\infty}^{2}}c_{\infty}^{2}\frac{t_{\rm
cool}^{2}}{r_{\rm cool}^{2}}=\frac{R_{\rm B}^{3}}{r_{\rm
cool}^{2}\Theta^{2}},$ (37)

where $c_{\rm cool}=r_{\rm cool}/t_{\rm cool}$ is the adiabatic sound speed of
the cooled fluid. Eq. (37) can be solved to give $r_{\rm
cool}=\Theta^{-2/3}R_{\rm B}$. Therefore, when adding the cooling effect, the
matter temperature will be given by

$T(r)=\begin{cases}T_{\rm cmb}\qquad&r\geq r_{\rm cool}\\\ T_{\rm
cmb}\left(\frac{\rho(r)}{\rho(r_{\rm cool})}\right)^{\gamma-1}\qquad&r<r_{\rm
cool},\end{cases}$ (38)

where $\rho(r)$ is obtained by solving Eq. (19). In our analysis, we will
incorporate this effect by shifting $R_{\rm B}\to r_{\rm cool}$ for values of
$\Theta>1$, as this is when the effect becomes relevant.

####  2.3.2 Ionisation of matter

Independently of the ionisation fraction at infinity, $\bar{x}_{e}$, the
rising energy of the accreted matter inside the EDO can provide enough energy
for the electrons in the hydrogen atoms to escape via scattering processes
such as $H+H\to H+e^{-}+p+\gamma$, what is known as collisional ionisation.
The minimum temperature at which this effect comes into play depends on the
properties of the fluid at infinity; however, following Ref. Bai  _et al._
(2020), we will take it to be $T_{\rm ion}\simeq 1.5\times
10^{4}~{}\text{K}\approx 1.3~{}\text{eV}$ Nobili  _et al._ (1991).

Therefore, the temperature of the fluid will evolve following Eq. (26) until
it reaches the first threshold, $r_{\rm start}$, given by $T(r_{\rm
start})=T_{\rm ion}=1.3~{}\text{eV}$. At that point, the accreted matter
collisionally ionizes until $r_{\rm end}$, defined by $x_{e}(r_{\rm end})=1$.
Inside this region of space, the temperature of the fluid can be taken to
remain constant as this process only redistributes the internal energy of the
gas, while the evolution of the ionisation fraction $x_{e}(r)$ can be
calculated using the first law of thermodynamics, leading to

$\frac{{\rm d}}{{\rm d}r}\left(\frac{3}{2}[1+x_{e}(r)]T(r)-[1-x_{e}(r)]E_{\rm
I}\right)=-[1+x_{e}(r)]T(r)\rho(r)\frac{{\rm d}\left(1/\rho(r)\right)}{{\rm
d}r},$ (39)

where $E_{\rm I}\simeq 13.6~{}\text{eV}$ is the binding energy of the neutral
hydrogen. Solving this equation gives

$x_{e}(r)=(1+\bar{x}_{e})\left(\frac{\rho(r)}{\rho(r_{\rm
start})}\right)^{\left(\frac{3}{2}+\frac{E_{\rm I}}{T_{\rm
ion}}\right)^{-1}}-1,$ (40)

where the $\rho(r_{\rm start})$ denominator ensures the matching with boundary
condition $x_{e}(r_{\rm start})=\bar{x}_{e}$. Once $x_{e}(r)$ reaches
$x_{e}(r_{\rm end})=1$, the temperature of the accreted matter unfreezes and
evolves similarly to Eq. (26), ensuring continuity. Therefore, we can account
for collisional ionisation as follows:

$\displaystyle T(r)=$ $\displaystyle\begin{cases}T_{\rm cmb}\qquad&r\geq
r_{\rm cool}\\\ T_{\rm cmb}\left(\frac{\rho(r)}{\rho(r_{\rm
cool})}\right)^{\gamma-1}\qquad&r_{\rm cool}\geq r>r_{\rm start}\\\ T_{\rm
ion}\qquad&r_{\rm start}\geq r\geq r_{\rm end}\\\ T_{\rm
ion}\left(\frac{\rho(r)}{\rho(r_{\rm end})}\right)^{\gamma-1}\qquad&r_{\rm
end}<r\\\ \end{cases}$ (41) $\displaystyle x_{e}(r)=$
$\displaystyle\begin{cases}\bar{x}_{e}\qquad&r>r_{\rm start}\\\
(1+\bar{x}_{e})\left(\frac{\rho(r)}{\rho(r_{\rm
start})}\right)^{\left(\frac{3}{2}+\frac{E_{\rm I}}{T_{\rm
ion}}\right)^{-1}}-1\qquad&r_{\rm start}\geq r\geq r_{\rm end}\\\
1\qquad&r_{\rm end}>r.\end{cases}$ (42)

We expect the photons produced through collisional ionisation also to ionize
nearby neutral atoms when their energy is higher than the binding energy,
$E_{\rm I}$. This is called photoionisation  Osterbrock and Ferland (2006),
and is a second source of ionisation that can take place inside EDOs.
Realistically, both processes would occur simultaneously, but it is common
practice to calculate both cases separately, obtaining higher and lower bounds
for the system. Since the most conservative results are given by collisional
ionisation, we will just focus on that case in our analysis.

####  2.3.3 Relativistic electrons

Depending on the specific EDO properties, it is possible to heat the accreting
matter enough so that electrons become relativistic, for which we will take
$T(r)\geq 2m_{e}/3$. As introduced in section 2.1.1, the relativistic nature
of electrons has an impact on the pressure and the internal energy density of
the baryon fluid, controlled by the functions $f_{\mathcal{E}}(T/m_{e})$ and
$f_{P}(T/m_{e})$. However, in terms of the differential equation for the
density profile in Eq. (19), this difference arises through $\gamma$, which,
as defined in Eq. (16), evaluates to $5/3$ for non-relativistic temperatures,
and $13/9$ otherwise.

Moreover, as we can see in Eq. (26), the temperature relation to density also
depends on the adiabatic factor $\gamma$. However, this dependence perfectly
cancels the one in the differential equation from Eq. (19), implying that the
evolution of the temperature is not affected by relativistic effects.
Therefore, the only thing that we need to address is the energy density, for
which we can take two different routes: use Eq. (19) to solve the evolution of
$\rho(r)$ for the relativistic region numerically, or use the already
calculated value for $T(r)$, and rescale it to $\rho(r)$ using the correct
value for $\gamma=13/9$ inside the relativistic radius, $r_{\rm rel}$. Taking
the second route, the relativistic effects are incorporated as follows:

$\displaystyle T(r)=$ $\displaystyle\begin{cases}T_{\rm cmb}\qquad&r\geq
r_{\rm cool}\\\ T_{\rm cmb}\left(\frac{\rho_{\rm sol}(r)}{\rho_{\rm
sol}(r_{\rm cool})}\right)^{\gamma-1}\qquad&r_{\rm cool}\geq r>r_{\rm
start}\\\ T_{\rm ion}\qquad&r_{\rm start}\geq r\geq r_{\rm end}\\\ T_{\rm
ion}\left(\frac{\rho_{\rm sol}(r)}{\rho_{\rm sol}(r_{\rm
end})}\right)^{\gamma-1}\qquad&r_{\rm end}<r\\\ \end{cases}$ (43)
$\displaystyle x_{e}(r)=$
$\displaystyle\begin{cases}\bar{x}_{e}\qquad&r>r_{\rm start}\\\
(1+\bar{x}_{e})\left(\frac{\rho_{\rm sol}(r)}{\rho_{\rm sol}(r_{\rm
start})}\right)^{\left(\frac{3}{2}+\frac{E_{\rm I}}{T_{\rm
ion}}\right)^{-1}}-1\qquad&r_{\rm start}\geq r\geq r_{\rm end}\\\
1\qquad&r_{\rm end}>r\end{cases}$ (44) $\displaystyle\rho(r)=$
$\displaystyle\begin{cases}\rho_{\rm sol}(r)&r>r_{\rm rel}\\\ \rho_{\rm
sol}(r_{\rm rel})\left(\frac{\rho_{\rm sol}(r)}{\rho_{\rm sol}(r_{\rm
rel})}\right)^{3/2}&r_{\rm rel}\geq r,\end{cases}$ (45)

where we have called $\rho_{\rm sol}(r)$ to the solution obtained by solving
Eq. (19) for $\gamma=5/3$. Having this, the number of electrons is given by
computing Eq. (4) with the obtained functions, such that

$n_{e}(r)=\frac{1}{m_{p}}x_{e}(r)\rho(r).$ (46)

In Figure 2 we show some examples of the density, temperature, ionisation
fraction, and number of electrons for constant-density EDOs with various radii
and masses. Additionally, in Figure 3, we show a comparison of the effect that
different mass functions have on the accreted matter.

Figure 2: Solutions for the electron/proton number density (top-left), energy
density (top-right), temperature (bottom-left) and ionisation fraction
(bottom-right) for a uniform EDO of mass ${M=10^{5}~{}M_{\odot}}$ but
different radii. The dashed solution represents an EDO at $z=800$, while plain
corresponds to $z=200$, as given in Eq. (52). We can see the effect that
collisional ionisation has on temperature, and how the CMB cooling delays the
heating of the accreted matter for the $z=800$ solution.  Figure 3:
Temperature (left) and energy density (right) of accreted matter into EDOs
with same total mass ($M=10^{5}M_{\odot}$) but different mass functions. As
expected, all profiles agree up to the EDO radius, $R=10^{7}\,{\rm m}$, after
which each function leads to a different behaviour for the accreted matter.

Taking the series of approximations and correcting for thermodynamic effects
we have thus obtained the static behaviour of a baryonic fluid around dark
matter compact objects. In particular, we expect the heating and ionisation of
the accreted matter to have an impact on the thermal history of the universe.
In the next section, we will study EDOs in this context, placing constraints
on their formation.

##  3 Ionisation of the cosmic microwave background

Before addressing the effect that EDOs have on the ionisation history of the
universe, we first need to define the boundary conditions at infinity for the
EDOs. Here we will use the Peebles Case B Peebles (1968); Poulin  _et al._
(2015) model of recombination, also called the three-level atom model, to
obtain the background evolution of baryon temperature ($T_{\rm M}$) and
ionisation fraction ($x_{e}$). This is described by the set of equations222We
recommend Ref. Baumann (2022) for a pedagogical introduction to these
equations.

$\displaystyle\frac{{\rm d}T_{\rm M}}{{\rm d}z}$
$\displaystyle=\frac{1}{(1+z)}\left[2T_{\rm M}+\frac{8\pi^{2}\sigma_{\rm
T}{T_{\rm cmb}}^{4}}{45H(z)m_{e}}\frac{x_{e}}{1+x_{e}}(T_{\rm M}-T_{\rm
cmb})\right],$ (47) $\displaystyle\frac{{\rm d}x_{e}}{{\rm d}z}$
$\displaystyle=C_{r}(z)\frac{\alpha_{\rm B}(T_{\rm
M})}{H(z)(1+z)}\left[nx_{e}^{2}+\left(\frac{m_{e}T_{\rm
M}}{2\pi}\right)^{3/2}e^{-\frac{E_{\rm I}}{T_{\rm M}}}(1-x_{e})\right].$ (48)

All the necessary definitions are described in what follows. $H(z)$ is the
Hubble parameter, where we will take Planck’s value for $H_{0}=73.2$ km/Mpc/s.
$\sigma_{\rm T}\simeq 6.65\times 10^{-25}\,{\rm cm}^{2}$ is the Thomson’s
constant, $E_{\rm I}=13.2$ eV is the binding energy for the ground state of
the hydrogen atom, $T_{\rm cmb}\simeq 2.35\times 10^{-4}\,(1+z)\,{\rm eV}$ is
the already mentioned CMB radiation temperature and $n(z)=\rho(z)/m_{\rm p}$
is the baryon density number. Finally, $C_{r}(z)$ is the effective branching
ratio (or the so-called Peebles factor), and it describes the probability that
excited electrons in a hydrogen atom will reach the ground state. This factor
is given by

$C_{r}(z)=\frac{1+K_{H}\Lambda_{H}n(1-x_{e})}{1+K_{H}(\Lambda_{H}+\beta_{H}(T_{\rm
M}))n(1-x_{e})}$ (49)

where $K_{H}={\lambda_{\rm Ly}}^{3}(8\pi H(z))$ with $\lambda_{\rm Ly}=121.5$
nm being wavelength of a Lyman-$\alpha$ photon,
$\Lambda_{H}=8.227~{}\text{s}^{-1}$ is is the two-photon decay rate of the
hydrogen 2S state, and

$\displaystyle\beta_{H}(T_{\rm M})$ $\displaystyle=\alpha_{\rm B}(T_{\rm
M})\left(\frac{m_{e}T_{\rm M}}{2\pi}\right)^{3/2}e^{3E_{\rm I}/4T_{\rm M}},$
(50) $\displaystyle\alpha_{\rm B}(T_{\rm M})$
$\displaystyle=9.8\frac{\alpha^{2}}{m_{e}^{2}}\left(\frac{E_{\rm I}}{T_{\rm
M}}\right)\log\left(\frac{E_{\rm I}}{T_{\rm M}}\right)$ (51)

are characteristic parameters for the ionisation and recombination rates,
respectively, with $\alpha=1/137$ being the fine-structure constant.

Therefore, we will use the solution for this system of equations as our
boundary conditions at infinity of the EDOs, such that

$\displaystyle\rho_{\infty}$ $\displaystyle\equiv m_{p}n(z),$ $\displaystyle
T_{\infty}\equiv$ $\displaystyle~{}T_{\rm M}(z),$
$\displaystyle\bar{x}_{e}\equiv x_{e}(z).$ (52)

However, if an EDO is placed in this background, we expect it to backreact
onto the baryon fluid, modifying the recombination history as given by Eqs.
(47-48). For this, we will calculate the luminosity and the energy deposition
rate of these objects into the background, allowing us to obtain the
modifications to the ionisation fraction of the universe.

###  3.1 Luminosity and energy deposition rate

For a spherically symmetric extended dark matter object such as the ones
described in section 2, in-falling electrons and protons may scatter in such a
way that can lead to bremsstrahlung radiation via processes such as $ee\to
ee\gamma$ or $ep\to ep\gamma$. The emissivity of this radiation, defined as
the emission power per volume per frequency per steradian, into photons of
frequency $\nu$ at a distance $r$ of the EDO is given by

$j_{\nu}(r)=\frac{8}{3}\left(\frac{2\pi
m_{e}}{3T(r)}\right)^{1/2}\frac{\alpha^{3}}{m_{e}^{2}}g_{ff}(\nu,T(r))e^{-2\pi\nu/T(r)}n_{e}(r)n_{p}(r),$
(53)

where $g_{ff}(\nu,T(r))$ is the free-free Gaunt factor that adds in the
quantum corrections, and $T(r)$ and $n_{e}(r)=n_{p}(r)$ are given by Eqs.(43)
and (46), respectively.

Integrating the frequency in the emissivity we obtain the radiation power
density, defined as the radiation per power per volume, which can be expressed
as Stepney and Guilbert (1983)

$\mathcal{L}(r)=n_{e}(r)n_{p}(r)\alpha\sigma_{\rm
T}T(r)\mathcal{J}(T(r)/m_{e}),$ (54)

where $\mathcal{J}(T(r)/m_{e})$ is a dimensionless function encoding the
information from the free-free Gaunt function. Considering both the formulas
for $e-p$ Svensson (1982) and $e-e$ Nozawa  _et al._ (2009) free-free
emissivities, Ref. Ali-Haïmoud and Kamionkowski (2017) provides a percent
accuracy level approximation given by

$\mathcal{J}(X)\approx\begin{cases}\frac{4\sqrt{2}}{\pi\sqrt{\pi}}X^{-1/2}(1+5.5X^{1.25}),\qquad&X<1,\\\
\frac{27}{2\pi}\left[\log\right(2Xe^{-\gamma_{\rm
E}}+0.08\left)+\frac{4}{3}\right],\qquad&X>1,\end{cases}$ (55)

where $\gamma_{\rm E}=0.557$ is the Euler Gamma constant. Finally, by
integrating over volume the radiation power density we obtain the total
luminosity. However, using Eq. (54) to calculate the luminosity will give us
an overestimate, since we are only interested in the luminosity difference due
to the addition of the EDO. Therefore, we will use the following definition of
the EDO accretion luminosity,

$L=\int_{0}^{\infty}{\rm d}r\,4\pi
r^{2}\left[\mathcal{L}(r)-\mathcal{L}(\infty)\right],$ (56)

which vanishes in the absence of an extended dark matter object. Notice that
our static solution for the accreted matter into an EDO is consistent only up
to $r=R_{B}$, placing an upper cut-off to this integral. Finally, the power
density at redshift $z$ is given by

$P(z)=L(z)n_{\rm EDO}(z),$ (57)

where the redshift dependence comes from the boundary conditions at infinity
as defined in Eq. (52), and $n_{\rm EDO}(z)$ is the number density of extended
dark objects at a given redshift, defined as $n_{\rm EDO}(z)=f_{\rm
DM}\,\rho_{\rm DM}(z)/M$, where $\rho_{\rm DM}(z)$ is the density of dark
matter and $f_{\rm DM}$ the fraction of dark matter covered by EDOs of same
total mass $M$.

###  3.2 Depositing EDOs in our universe: relative distances and velocities

So far, the solutions for extended dark matter objects we studied were
simplified so that they were isolated and in the same frame of reference as
the thermal bath. In this subsection, we will address these approximations,
getting a more realistic solution for their effects on the CMB.

To make sure that the isolation approximation is valid, we need to compare the
Bondi radius for the EDOs, defined in Eq. (32), with their average separation
as a dark matter candidate, $d_{\rm s}$. Using the dark matter density, we
find

$d_{s}(z)\sim\left(\frac{M}{\rho_{\rm DM}(z)f_{\rm DM}}\right)^{1/3}.$ (58)

Therefore, the condition $d_{s}(z)\gg R_{\rm B}$ translates into

$M\ll\sqrt{\frac{(\gamma_{\infty}T_{\infty}(z))^{3}}{(2Gm_{p})^{3}\rho_{\rm
DM}(z)f_{\rm DM}}}\overset{z\to 50}{\approx}\frac{3\times
10^{4}M_{\odot}}{\sqrt{f_{\rm DM}}},$ (59)

which places an upper limit in our analysis, depending on the dark matter
fraction corresponding to EDOs.

Taking into account relative velocities is more complicated. As one might
expect, an EDO with a high relative speed compared to the speed of sound of
baryon fluid, $v_{\rm rel}>c_{s}$, will be more evasive, leading to a smaller
total accretion of matter and therefore emitted radiation. To deal with this
scenario we will take the approach from Refs.Ali-Haïmoud and Kamionkowski
(2017); Bondi (1952); Bondi and Hoyle (1944), where these effects are
considered by shifting the speed of sound at infinity by
$c_{\infty}\to\sqrt{c_{\infty}^{2}+v_{\rm rel}^{2}}$, decreasing the Bondi
radius for the extended dark matter object. We thus expect $v_{\rm rel}$ to be
a stochastic variable, and take it to be controlled by a three-dimensional
Gaussian linear distribution with a standard deviation $\langle
v_{s}^{2}\rangle^{1/2}$.333For large relative speeds, it can be argued that
the spherical accretion approximation should break down. Instead, disk-like
accretion of baryonic matter should dominate above a certain threshold for
$v_{\rm rel}$ Poulin  _et al._ (2017); Serpico  _et al._ (2020). Including
this effect would lead to tighter bounds on EDOs. In this work, we remain
conservative and focus on this simplified system.,444Additionally, we do not
consider the effect of outflows. See Bosch-Ramon and Bellomo (2020) for a
discussion about the appropriateness of accretion approximations in light of
potential outflows, which may weaken the constraints, particularly for
supersonic relative velocities. Following Ref. Ali-Haïmoud and Kamionkowski
(2017), the redshift dependence of the standard deviation is given by Dvorkin
_et al._ (2014)

$\langle v_{s}^{2}\rangle^{1/2}=\min[1,z/10^{3}]\times 30\text{km/s}.$ (60)

For a realistic estimation of the energy deposition rate of EDOs into the CMB
radiation, we will have to calculate the weighted power density

$\langle P(z)\rangle=\langle L(z)\rangle n_{\rm EDO}(z),$ (61)

where the brackets include the effect of the relative velocities when applied
on a function $\mathcal{O}$ as follows

$\langle\mathcal{O}\rangle=\frac{4\pi}{(2\pi\langle
v_{s}^{2}\rangle/3)^{3/2}}\int_{0}^{\infty}{\rm d}v_{\rm rel}\,v_{\rm
rel}^{2}e^{-\frac{v_{\rm rel}^{2}}{2\langle
v_{s}^{2}\rangle/3}}\mathcal{O}|_{c_{\infty}\to\sqrt{c_{\infty}^{2}+v_{\rm
rel}^{2}}}.$ (62)

Figure 4: Relative-velocity averaged luminosity of a uniform EDO for different
masses (left), normalized over the Eddington luminosity, given by
$L_{\rm{Edd}}=4\pi GMm_{p}/\sigma_{\rm T}$, defined as the luminosity at which
radiation pressure due Thomson scattering balances the accretion of hydrogen
due to gravity. On the right, we plot the total averaged luminosity for the
different EDO profiles with the same total mass $M=10^{2}M_{\odot}$, also
normalized over the Eddington luminosity. Although the accretion profiles from
Figure 3 differed greatly inside the EDO radius, we can see that changing the
mass function is to very good approximation equivalent to shifting the uniform
sphere solution, as represented by the dashed and dot-dashed lines in the
left-hand side plot.

Numerical calculations of this expression for the total luminosity are
computationally expensive. In Appendix A we show how this process can be made
more efficient by conveniently rescaling the radial components of the accreted
matter profiles.

In Figure 4 we show the relative velocity integrated luminosity for EDOs with
different radii and mass functions. As we can see, the total luminosity from
any mass function can be obtained to very good approximation by just shifting
the uniform sphere solution. With this, we are ready to quantify how the
accretion of matter into extended dark matter objects changes the ionisation
fraction of the background evolution.

###  3.3 Effect on CMB anisotropies

After recombination at redshift $z\sim 1100$, it is assumed that most
electrons and protons are combined into hydrogen, leading to an ionisation
fraction of $x_{e}\sim 10^{-4}-10^{-3}$. This background neutrality allows the
CMB radiation to move freely through space until today, except for the
ignition of stars at $z\sim\mathcal{O}(10)$. Therefore, strong constraints can
be placed using the state of the universe through the so-called dark ages, as
any deviation from neutrality would be a possible source for the absorption of
photons, decreasing the visibility of CMB anisotropies. This is where the
existence of EDOs comes into play, as their emitted radiation does not just
overlay the CMB, but also heats the matter around it, ionizing the neutral
hydrogen and locally causing $x_{e}\to 1$. To place constraints, we need to
find how this local change in the ionisation fraction affects the global
evolution of the universe as studied in the Peebles case B recombination case
in Eqs. (47) and (48). As derived in Ref. Ali-Haïmoud and Kamionkowski (2017),
the fraction of the power density that is deposited into the background,
called the energy deposition rate, $\dot{\rho}_{\rm dep}(t)$, is given by

$a^{-7}\frac{{\rm d}}{{\rm d}t}\left(a^{7}\dot{\rho}_{\rm dep}\right)\approx
0.1\,n\,\sigma_{\rm T}(\langle P\rangle-\dot{\rho}_{\rm dep}),$ (63)

where the $0.1$ prefactor comes from the efficiency of the Compton scattering.
We see that for $0.1\,n\,\sigma_{\rm T}a^{7}>1$ all the energy produced by the
EDO gets deposited in the neutral background, until $z\approx 200$, after
which the energy injection decays as $a^{-7}$.

Figure 5:  Modification of the ionisation fraction due to the presence of
uniform-sphere EDOs. On the left, the ionisation fraction is presented against
the $\Lambda$CDM evolution (blue) obtained solving Eqs. (47-48), while on the
right we show the net modification for each set of parameters.

Once we obtain the total energy deposited by EDOs at any redshift into the
background, we are ready to see how they will impact the average temperature
and ionisation fraction of the universe. For this, we just need to modify the
Peeble’s equations from Eqs. (47) and (48) as follows

$\displaystyle\frac{{\rm d}T_{\rm M}}{{\rm d}z}$
$\displaystyle=\frac{1}{(1+z)}\left[2T_{\rm M}+\frac{8\pi^{2}\sigma_{\rm
T}{T_{\rm cmb}}^{4}}{45H(z)m_{e}}\frac{x_{e}}{1+x_{e}}(T_{\rm M}-T_{\rm
cmb})\right]-\frac{2}{3\,n}\frac{1+2x_{e}}{3H(z)(1+z)}\dot{\rho}_{\rm dep},$
(64) $\displaystyle\frac{{\rm d}x_{e}}{{\rm d}z}$
$\displaystyle=C_{r}(z)\frac{\alpha_{\rm B}(T_{\rm
M})}{H(z)(1+z)}\left[nx_{e}^{2}+\left(\frac{m_{e}T_{\rm
M}}{2\pi}\right)^{3/2}e^{-\frac{E_{\rm I}}{T_{\rm
M}}}(1-x_{e})\right]-\frac{1-x_{e}}{3H(z)(1+z)}\frac{\dot{\rho}_{\rm
dep}}{E_{\rm I}\,n}.$ (65)

In Figure 5, we show how the ionisation fraction is modified by the addition
of EDOs with different masses, radii, and dark matter fractions. To constrain
these objects against their impact on the CMB anisotropies, we would have to
implement the modified $x_{e}(z)$ in a Boltzman code and solve for the CMB
power spectra, as it was done for PBHs in Refs. Ali-Haïmoud and Kamionkowski
(2017); Aloni  _et al._ (2017). However, following Ref. Bai  _et al._ (2020),
we note that we may also compare our EDO results to the existing results from
PBHs without explicitly repeating this analysis. In this way, we will impose
the constrain that $\Delta x_{e}/x_{e}<\mathcal{O}(1)$ at all times, testing
particularly for different masses, radii, and dark matter fractions satisfying
$\Delta x_{e}(z=50)<10^{-4}$.555It is reasonable to expect that with better
observations this threshold will be constrained to smaller values for $\Delta
x_{e}(z=50)$. Since the corrections to the recombination equations from Eqs.
(64-65) are subdominant, we can generally predict the impact that a change in
$\Delta x_{e}(z=50)$ would have on the EDO constrains by just conveniently
shifting $f_{\rm DM}$, as they are related by $\delta\Delta x_{e}/\Delta
x_{e}=\delta f_{\rm DM}/f_{\rm DM}$.

We present our results in Figure 6. On the plot on the left, we express the
EDO radius which leads to large $\Delta x_{e}$ as a function of the mass,
while keeping $f_{\rm DM}$ constant. We can see that the contour for $f_{\rm
DM}=1$ perfectly agrees with the result from Ref. Bai  _et al._ (2020), while
we extend the calculations also to present constraints for $f_{\rm
DM}=\\{10^{-1},10^{-2},10^{-3}\\}$, finding a logarithmic shift of the
contours to higher masses for the same radii. On the right, we keep the radius
constant, constraining EDOs with a given mass constituting some fraction of
the total dark matter, going down to $f_{\rm DM}=10^{-5}$. These results
complete the impact that extended dark objects may have on the CMB
anisotropies, which can now be added to the rest of constraints via different
tests, as will be discussed in the next section.

Figure 6: Constraints on uniform density EDOs for different dark matter
fractions as a function of mass vs radius (right) and for different radii as a
function of mass vs dark matter fraction (left). The gray regions mark the
limitations of our analysis. Figure 7: Complementarity of constraints on the
maximum allowed dark matter fraction as a function of mass, for boson stars
and NFW subhaloes of particular sizes. Shown here are constraints from the
non-observation of microlensing signals Croon  _et al._ (2020a, b) in the
EROS-2 Alcock  _et al._ (1998b), OGLE-IV Niikura  _et al._ (2019b) (note that
this does not include the latest results based on a combination of the OGLE-
III and OGLE-IV datasets Mroz  _et al._ (2024)), and Subaru-HSC Niikura  _et
al._ (2019c) datasets; gravitational wave signals from the inspiral regime of
binary mergers assuming binary formation in the early Universe (Croon  _et
al._ (2023)), and collisional ionisation of the CMB due to accretion (this
work). The gray regions show the limitations of our analysis.

##  4 Discussion

Extended dark matter objects such as boson stars, subhaloes and miniclusters
are popular dark matter candidates, associated with a variety of creation
mechanisms and potential mass ranges. So far, gravitational microlensing and
gravitational waves from binary mergers have been used to constrain a wide
range of the EDO parameter space $\lesssim 10\,M_{\odot}$. At large masses,
the ionisation of the CMB due to the accretion of baryonic matter, which for
PBHs implies the leading constraint, can be used to constrain the dark matter
fraction of these extended objects. In Ref. Bai  _et al._ (2020), the authors
calculated the bounds for EDOs with uniform-sphere mass functions making up
the $100\%$ of dark matter. In this work we extend this analysis in two
important ways: we study different mass functions, and allow for these objects
to be a subcomponent of DM. Specifically, besides uniform spheres (20) we
focus on mini-haloes (23), boson stars, and NFW subhaloes (21). This allows us
to place large-mass constraints on EDOs of these classes. We explain how our
analysis can be easily extended to account for other dark objects.

The impact of different mass profiles is rather subtle, and primarily affects
the interpretation of constraints. As described in Figure 4, in the scenarios
studied in this work changing the mass function is equivalent to rescaling the
radius for the uniform sphere. However, the rescaling of the constraints
between mass functions is model-dependent and non-trivial, so one first needs
to calculate the averaged luminosity for each case to find the ratio between
them. Once this is obtained, the constraints from Figure 6 can be mapped into
the new rescaled radii. To ease up this process, we provide the Mathematica
notebook used to numerically generate the constraints for the uniform
sphere,666The repository can be found in the following link:
https://gitlab.com/SergioSevillano/edo-accretion. which works on any mass
function as long as it can be expressed as $M(r/R)$ where $R$ is defined such
that $M(1)=M_{\rm tot}$ (see Appendix A).

The constraints we find are complementary to those found through the non-
observation of gravitational microlensing events and binary mergers. In Figure
7 we show these complementary constraints for boson stars and NFW subhaloes
together with the CMB constraints derived in this work. As the gravitational
wave constraints were derived assuming the inspiral regime, in which the
waveform of these object binaries would be known, the objects in the plot have
too low of a compactness ($C\lesssim 10^{-3}$) to be observed by the LVK
collaboration, as the inspiral regime would have ended at frequencies outside
of its observational window. The apparent open window for stellar mass objects
with stellar radius can be closed through gravitational heating of dwarf
galaxies Graham and Ramani (2023).

In summary, this study extends the exploration of DM by introducing
constraints on large mass EDOs, setting upper limits on their allowed mass
fraction for a large range of masses and sizes. We have worked out constraints
on boson stars, NFW subhaloes, and miniclusters in detail, and provided a
Mathematica package to generalise our results to arbitrary other EDO profiles.
Our results complement existing constraints derived from the non-observation
of gravitational microlensing, binary mergers, and dynamical heating of dwarf
galaxies.

## Acknowledgements

We thank Joachim Kopp, Vivian Poulin, and Harikrishnan Ramani for useful
discussions. DC and SSM are supported by the STFC under Grant No.
ST/T001011/1.

##  Appendix A Improving the efficiency of the integral over relative
velocities

In section 2 we showed how to obtain the temperature, density, and ionisation
fraction for the accreted matter around an EDO, and in section 3 the impact
that they have on recombination history. Allowing for generic mass functions
requires solving this system numerically, which is computationally expensive,
especially when averaging the Luminosity for relative speeds of the EDOS. This
requires performing the Gaussian integral over $v_{\rm rel}$ from Eq. (66),

$\langle\mathcal{L}\rangle=\frac{4\pi}{(2\pi\langle
v_{s}^{2}\rangle/3)^{3/2}}\int_{0}^{\infty}{\rm d}v_{\rm rel}\,v_{\rm
rel}^{2}e^{-\frac{v_{\rm rel}^{2}}{2\langle
v_{s}^{2}\rangle/3}}\mathcal{L}|_{c_{\infty}\to\sqrt{c_{\infty}^{2}+v_{\rm
rel}^{2}}},$ (66)

where, a priori, for every evaluated value of the integrated, we would need to
follow all the steps from section 2 with a shift in the sound speed at
infinity by $c_{\infty}^{2}\to c_{\infty}^{2}+v_{\rm rel}^{2}$. However, in
this appendix, we will show how one can recycle EDO solutions for this problem
so that this integral can be solved more efficiently.

This shift in the sound speed leads to a shift of the Bondi radius, $R_{\rm
B}\to GM/(c_{\infty}^{2}+v_{\rm rel}^{2})$, which only impacts the density
equation,

$\frac{GM(r)}{r^{2}}+\gamma K\rho(r)^{\gamma-2}\frac{d\rho(r)}{dr}=0,$ (67)

through the constant factor $K$, defined outside of the ionizing region as

$K=\frac{GM}{\gamma_{\infty}R_{\rm B}\rho_{\infty}^{\gamma_{\infty}-1}}.$ (68)

Figure 8: Example for the rescaling technique to save computing time when
evalutating different relative velocities. The top two plots show the system
with a rescaled radial component, $r^{\prime}=rR_{\rm B}$, while the plots on
the bottom correspond to the real solutions, obtained by inversing the
rescaling for $r^{\prime}$.

Therefore, in what follows, we will show that Eq. 70 is independent of the
Bondi radius by conveniently rescaling most quantities. For this, we will have
to split the calculation into two regions

  * •

Outside the EDO ($r>R$): This region is defined as where the mass function
becomes a constant ($M(r>R)=M)$. Thus, inserting $K$ into Eq. (70) and
rescaling $r\to R_{\rm B}r^{\prime}$ and
$\rho(r)\to\rho_{\infty}\rho^{\prime}(r^{\prime})$ gives

$\frac{1}{r^{\prime
2}}+\rho^{\prime}(r^{\prime})^{\gamma-2}\frac{d\rho^{\prime}(r^{\prime})}{dr^{\prime}}=0,$
(69)

which can be solved for $r^{\prime}\leq 1$, and matched at the boundary
($r^{\prime}=R/R_{\rm B}$) with the interior solution.

  * •

Inside the EDO ($r\leq R$): In this case, we will need to take into account
the radial dependence of the mass function. The rescaling is possible as long
as the mass function can be expressed as a function of $r/R$, which is
generally the case. In this way, rescaling $r\to\hat{r}R$,
$\rho(r)\to\rho_{\infty}\hat{\rho}(\hat{r})R^{\frac{1}{\gamma_{\infty}-1}}$
and $M(\hat{r})\to\widetilde{M}(\hat{r})M_{\rm tot}$, we obtain

$\frac{\widetilde{M}(\hat{r})}{\hat{r}^{2}}+\hat{\rho}(\hat{r})^{\gamma-2}\frac{d\hat{\rho}(\hat{r})}{d\hat{r}}=0,$
(70)

which again does not explicitly depend on $R_{\rm B}$.

Therefore, the idea is to solve these equations without specifying $R_{\rm
B}$, and then rescale and match the solutions at the boundary for each
evaluated $v_{\rm rel}$ in the integral. Moreover, since the temperature of
the accreted matter and its density are related via

$T(r)=T_{\infty}\left(\frac{\rho(r)}{\rho_{\infty}}\right)^{\gamma-1},$ (71)

and this does not depend on $R_{\rm B}$ either, we find that the rescaled
quantities corresponding to the ionisation start and end, $r^{\prime}_{\rm
start}=r_{\rm start}/R_{\rm B}$ and $r^{\prime}_{\rm end}=r_{\rm end}/R_{\rm
B}$, respectively, will not depend on $R_{\rm B}$ either, so they can also be
rescaled after being computed for the general case. However, given that the
relation between $T(r)$ and $\rho(r)$ depends on their values at infinity, we
will have to compute the whole system for different redshifts.

In Figure 8 we show an example of this rescaling for different relative
velocities (upper plots); by rescaling $r\to r/R_{\rm B}$, the solutions
perfectly agree up to the beginning of the different interior regions.
However, this is only because in the rescaled units the boundary to the
interior region is $R/R_{\rm B}$, which depends on the relative velocity. All
lines correspond to EDOs with the same mass and radii, as can be seen in the
un-rescaled plots at the bottom.

## References

  * Alcock  _et al._ (1998a) C. Alcock, R. A. Allsman, D. Alves, R. Ansari, É. Aubourg, T. S. Axelrod, P. Bareyre, J.-P. Beaulieu, A. C. Becker, D. P. Bennett, S. Brehin, F. Cavalier, S. Char, K. H. Cook, R. Ferlet, J. Fernandez, K. C. Freeman, K. Griest, P. Grison, M. Gros, C. Gry, J. Guibert, M. Lachièze-Rey, B. Laurent, M. J. Lehner, É. Lesquoy, C. Magneville, S. L. Marshall, É. Maurice, A. Milsztajn, D. Minniti, M. Moniez, O. Moreau, L. Moscoso, N. Palanque-Delabrouille, B. A. Peterson, M. R. Pratt, L. Prévôt, F. Queinnec, P. J. Quinn, C. Renault, J. Rich, M. Spiro, C. W. Stubbs, W. Sutherland, A. Tomaney, T. Vandehei, A. Vidal-Madjar, L. Vigroux, and S. Zylberajch, The Astrophysical Journal 499, L9 (1998a).
  * Niikura  _et al._ (2019a) H. Niikura _et al._ , Nature Astron. 3, 524 (2019a), arXiv:1701.02151 [astro-ph.CO] .
  * Niikura  _et al._ (2019b) H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki, Phys. Rev. D 99, 083503 (2019b), arXiv:1901.07120 [astro-ph.CO] .
  * Green and Kavanagh (2021) A. M. Green and B. J. Kavanagh, J. Phys. G 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Croon  _et al._ (2020a) D. Croon, D. McKeen, N. Raj, and Z. Wang, Phys. Rev. D 102, 083021 (2020a), arXiv:2007.12697 [astro-ph.CO] .
  * Bai  _et al._ (2020) Y. Bai, A. J. Long, and S. Lu, JCAP 09, 044 (2020), arXiv:2003.13182 [astro-ph.CO] .
  * Croon  _et al._ (2020b) D. Croon, D. McKeen, and N. Raj, Phys. Rev. D 101, 083013 (2020b), arXiv:2002.08962 [astro-ph.CO] .
  * Crispim Romão and Croon (2024) M. Crispim Romão and D. Croon, (2024), arXiv:2402.00107 [astro-ph.CO] .
  * Bird  _et al._ (2016) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Ali-Haïmoud  _et al._ (2017) Y. Ali-Haïmoud, E. D. Kovetz, and M. Kamionkowski, Phys. Rev. D 96, 123523 (2017), arXiv:1709.06576 [astro-ph.CO] .
  * Kavanagh  _et al._ (2018) B. J. Kavanagh, D. Gaggero, and G. Bertone, Phys. Rev. D 98, 023536 (2018), arXiv:1805.09034 [astro-ph.CO] .
  * Bertone  _et al._ (2020) G. Bertone _et al._ , SciPost Phys. Core 3, 007 (2020), arXiv:1907.10610 [astro-ph.CO] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 123, 161102 (2019), arXiv:1904.08976 [astro-ph.CO] .
  * Chen and Huang (2020) Z.-C. Chen and Q.-G. Huang, JCAP 08, 039 (2020), arXiv:1904.02396 [astro-ph.CO] .
  * Franciolini (2021) G. Franciolini,  _Primordial Black Holes: from Theory to Gravitational Wave Observations_ , Ph.D. thesis, Geneva U., Dept. Theor. Phys. (2021), arXiv:2110.06815 [astro-ph.CO] .
  * Croon  _et al._ (2023) D. Croon, S. Ipek, and D. McKeen, Phys. Rev. D 107, 063012 (2023), arXiv:2205.15396 [astro-ph.CO] .
  * Carr (1981) B. Carr, Monthly Notices of the Royal Astronomical Society 194, 639 (1981).
  * Ricotti  _et al._ (2008) M. Ricotti, J. P. Ostriker, and K. J. Mack, Astrophys. J. 680, 829 (2008), arXiv:0709.0524 [astro-ph] .
  * Ali-Haïmoud and Kamionkowski (2017) Y. Ali-Haïmoud and M. Kamionkowski, Phys. Rev. D 95, 043534 (2017), arXiv:1612.05644 [astro-ph.CO] .
  * Poulin  _et al._ (2017) V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, Phys. Rev. D 96, 083524 (2017), arXiv:1707.04206 [astro-ph.CO] .
  * Graham and Ramani (2023) P. W. Graham and H. Ramani, (2023), arXiv:2311.07654 [hep-ph] .
  * Witten (1984) E. Witten, Phys. Rev. D 30, 272 (1984).
  * Zhitnitsky (2003) A. R. Zhitnitsky, JCAP 10, 010 (2003), arXiv:hep-ph/0202161 .
  * Schunck and Mielke (2003) F. E. Schunck and E. W. Mielke, Class. Quant. Grav. 20, R301 (2003), arXiv:0801.0307 [astro-ph] .
  * Bar  _et al._ (2018) N. Bar, D. Blas, K. Blum, and S. Sibiryakov, Phys. Rev. D 98, 083027 (2018), arXiv:1805.00122 [astro-ph.CO] .
  * Navarro  _et al._ (1996) J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 462, 563 (1996), arXiv:astro-ph/9508025 [astro-ph] .
  * Ricotti and Gould (2009) M. Ricotti and A. Gould, Astrophys. J. 707, 979 (2009), arXiv:0908.0735 [astro-ph.CO] .
  * Bertschinger (1985) E. Bertschinger, ApJS 58, 39 (1985).
  * Bondi (1952) H. Bondi, Mon. Not. Roy. Astron. Soc. 112, 195 (1952).
  * Nobili  _et al._ (1991) L. Nobili, R. Turolla, and L. Zampieri, (1991).
  * Osterbrock and Ferland (2006) D. E. Osterbrock and G. J. Ferland,  _Astrophysics of gaseous nebulae and active galactic nuclei_ (2006).
  * Peebles (1968) P. J. E. Peebles, ApJ 153, 1 (1968).
  * Poulin  _et al._ (2015) V. Poulin, P. D. Serpico, and J. Lesgourgues, JCAP 12, 041 (2015), arXiv:1508.01370 [astro-ph.CO] .
  * Baumann (2022) D. Baumann,  _Cosmology_ (Cambridge University Press, 2022).
  * Stepney and Guilbert (1983) S. Stepney and P. W. Guilbert, MNRAS 204, 1269 (1983).
  * Svensson (1982) R. Svensson, ApJ 258, 335 (1982).
  * Nozawa  _et al._ (2009) S. Nozawa, K. Takahashi, Y. Kohyama, and N. Itoh, A&A 499, 661 (2009).
  * Bondi and Hoyle (1944) H. Bondi and F. Hoyle, MNRAS 104, 273 (1944).
  * Serpico  _et al._ (2020) P. D. Serpico, V. Poulin, D. Inman, and K. Kohri, Phys. Rev. Res. 2, 023204 (2020), arXiv:2002.10771 [astro-ph.CO] .
  * Bosch-Ramon and Bellomo (2020) V. Bosch-Ramon and N. Bellomo, Astron. Astrophys. 638, A132 (2020), arXiv:2004.11224 [astro-ph.CO] .
  * Dvorkin  _et al._ (2014) C. Dvorkin, K. Blum, and M. Kamionkowski, Phys. Rev. D 89, 023519 (2014), arXiv:1311.2937 [astro-ph.CO] .
  * Aloni  _et al._ (2017) D. Aloni, K. Blum, and R. Flauger, JCAP 05, 017 (2017), arXiv:1612.06811 [astro-ph.CO] .
  * Alcock  _et al._ (1998b) C. Alcock, R. A. Allsman, D. Alves, R. Ansari, E. Aubourg, T. S. Axelrod, P. Bareyre, J.-P. Beaulieu, A. C. Becker, D. P. Bennett, and et al., The Astrophysical Journal 499, L9–L12 (1998b).
  * Mroz  _et al._ (2024) P. Mroz _et al._ , (2024), arXiv:2403.02386 [astro-ph.GA] .
  * Niikura  _et al._ (2019c) H. Niikura _et al._ , Nat. Astron. 3, 524 (2019c), arXiv:1701.02151 [astro-ph.CO] .