# Structure functions from renormalization group improved small $x$ evolution

Wanchen Li  Department of Physics, Penn State University, University Park, PA
16802, USA Anna M. Staśto  Department of Physics, Penn State University,
University Park, PA 16802, USA

###### Abstract

We perform the fit to the structure function $F_{2}$ data from HERA including
terms due to the resummation at small $x$. The equation for the unintegrated
gluon density is solved, previously established within the renormalization
group improved small $x$ framework. We find very good description of the
structure function $F_{2}$ and its charm component $F_{2}^{c}$. The resulting
unintegrated gluon density is found to be consistent with the calculations
based on similar approaches available in the literature, with only slightly
higher intercept value.

##  1 Introduction

Deep Inelastic Scattering (DIS), the process of scattering of leptons off
protons is the most precise way to explore the nucleon structure. Since the
groundbreaking experiments at SLAC [1], which have discovered the partonic
structure of the proton, enormous progress has been achieved in understanding
the proton structure thanks to the series of experiments as well as
theoretical developments. HERA collider in Hamburg, the only electron-proton
collider up to date, operated at around $\sqrt{s}=300\;\rm GeV$ of center-of-
mass energy, and provided important insight into the behavior of strong
interactions in the small $x$ Bjorken region. One of the most important
discoveries by HERA experiments was the observation of the strong rise of the
structure function $F_{2}$ in the small $x$ region [2, 3]. This rise is
understood to be driven by the sharp rise of the gluon density towards small
$x$, which is a fundamental property of the strong interactions.

The Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [4, 5, 6, 7] evolution
equation together with the collinear factorization [8] provide a framework for
the perturbative parton density evolution. The DGLAP evolution runs with the
scale $Q^{2}$ thus it resums the powers of $\ln\left(Q^{2}/Q^{2}_{0}\right)$
where $Q_{0}$ is some reference scale. Currently, DGLAP splitting functions
are known up to NNLO (next-to-next-to-leading order) accuracy, and the global
analyses based on the wealth of deep inelastic and hadron-collider data
provide with the important information on parton densities.

The high energy framework provides an alternative approach to the parton
evolution. At sufficiently high center of mass energy $\sqrt{s}$, the
$\ln(s/s_{0})$ terms will become important and one needs a different
resummation approach than DGLAP. Here, $s_{0}$ is the energy scale choice,
which for the case of DIS would correspond to $Q^{2}$, the minus photon
virtuality. The Bjorken $x$ is then defined by $x\simeq Q^{2}/s$, where $s$ is
the photon-proton center of mass energy squared. The resummation in the high
energy limit is then accomplished by the Balitsky-Fadin-Kuraev-Lipatov (BFKL)
equation and thus accounts for the evolution in $\ln s/s_{0}$ or $\ln 1/x$ in
the specific case of DIS. The BFKL evolution equation has been derived at LL
(leading logarithmic in $\ln 1/x$) [9, 10, 11] and NLL [12, 13] accuracy. The
NNLL BFKL kernel has also been derived in the N=4 sYM case [14, 15, 16]. It
has also been derived (up to NLL accuracy in QCD) in the context of the high
energy operator expansion [17, 18, 19], dipole evolution [20, 21, 22, 23] and
in the Color Glass Condensate [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. In
this framework it includes also high density effects, or parton saturation
effects. These are stemming from the non-linear recombination of gluons (in
the context of CGC) or multiple scattering (in the context of the dipole
approach) which lead to the taming of the growth towards small $x$. Together
with $k_{T}$ factorization, the BFKL equation gives prediction for the
observable quantities, like structure functions. The basic object within this
framework is the unintegrated gluon density $f(x,k_{T}^{2})$ which depends on
the transverse momentum $k_{T}$ of the off-shell gluon.

The BFKL evolution at LL accuracy generates the rise of the unintegrated gluon
density with decreasing $x$, with the hard Pomeron intercept
$\omega_{0}=(3\alpha_{s}/\pi)\,4\ln 2$. However, it has been well known since
some time that this is too steep for the experimental data [35].

The NLL BFKL computed in [12, 13] and in [19, 33, 34] (in the context of high
energy processes with saturation) turned out to be very large and negative.
The result even leads to a instability, like negative cross section, and thus
hinted at the neccessity of the calculation of yet higher orders or a
resummation. The main part of the NLL corrections stem from the running of the
coupling, non-singluar terms of the DGLAP splitting function and kinematical
constraint. It can be shown that these contributions exhaust most of the NLL
correction [36, 37].

In order to stabilize the results, various resummation approaches were
constructed in the literature [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47,
48, 49, 50, 51] and their results have been successfully fitted to the
structure function data on $F_{2}$. In particular, the analysis by [51]
demonstrated that the HERA data might actually point to the neccessity of the
resummation, since the fits with the resummation resulted in slightly better
$\chi^{2}$ than the standard fits based on the NNLO DGLAP evolution. While not
a large effect, this may have consequences for processes at yet higher
energies. The resummation of low $x$ terms is especially important in the
context of the next generation high energy colliders as well as for
astroparticle physics. Several next generation Deep Inelastic Scattering
machines are planned, like Electron Ion Collider (EIC) in the US [52], Large
Hadron electron Collider (LHeC) [53], and further in the future and with yet
higher energy Future Circular Collider (FCC), in the electron-hadron option
[54, 55]. The EIC machine is the project based in the US, which will utilize
the existing RHIC complex and add the new electron accelerator. It will have
variable center of mass energy up to about $\sqrt{s}=\,140\,\rm GeV$, high
luminosity ${\cal L}=10^{34}\,\rm cm^{-2}s^{-1}$, and will have various
nuclear targets and polarization of the colliding beams. The LHeC project
plans to utilize the $7\,\rm TeV$ hadron beams from the LHC and collide them
with the electrons accelerated to about $50-60\,\rm GeV$ from the Energy
Recovery Linac (ERL). The center of mass energy would reach about $1\,\rm TeV$
in this collider. The LHeC would also have high luminosity and possibility of
performing electron-ion collisions with beams of lead, up to about $5\,\rm
TeV$ of center of mass energy. The FCC in electron-hadron mode would open up
another territory of the DIS. With the $50\,\rm TeV$ proton beams from the
FCC, colliding with the electron beams from ERL, the center of mass energies
up to $3.5\,\rm TeV$ would be possible, as well as the collisions with the
highly energetic lead beams. The precise knowledge of the evolution towards
small $x$ is thus necessary for these DIS projects as well as for the
processes which are sensitive to the small $x$ dynamics in the high energy
hadronic collisions, also in the context of the astroparticle physics, see
[56, 57].

In this paper, we focus on Ciafaloni-Colferai-Salam-Staśto (CCSS) resummation
[40, 42, 41] and perform the fits to the DIS structure function data based on
this framework. The CCSS resummation features the renormalization group
improved small $x$ equation. In the CCSS scheme the resummation is constructed
based on the collinear splitting function and exact BFKL up to the NLL
accuracy. It contains the kinematical constraint and appropriate subtractions
which avoid double counting. We find very good description of the small $x$
HERA data and the data on charm structure function. The CCSS resummation is
formulated entirely in the momentum space, and performing the fit to the HERA
data one can then extract the resummed unintegrated gluon distribution
function. The latter quantity can thus be used in other processes which make
use of the $k_{T}$ or hybrid factorization. Similar low $x$ resummations,
formulated entirely in the momentum space [58, 59] were successfully applied
to the description of the DIS data on structure functions. In particular, we
shall compare the resulting gluon distribution from CCSS to these earlier
extractions [59]. We find they are consistent, but there are some differences
related to the normalization and also small differences in the evolution
towards small $x$. In particular, CCSS is characterized by slightly faster
growth towards small $x$. This can have implications at yet lower values of
$x$, where saturation corrections (not currently included in CCSS formalism)
may be needed.

The outline of this paper is as follows. In Sec. 2, we introduce the BFKL
formalism and the CCSS resummation, recalling the formulae both in the Mellin
space and in the momentum space. The details of the setup for the calculation
of the structure function $F_{2}$ are presented in Sec. 3. In section Sec. 4
we show the results of the fits and also extracted unintegrated gluon density.
Finally, in Sec. Conclusions and outlook we state brief conclusions and
outlook.

##  2 Resummed Ciafaloni-Colferai-Salam-Stasto (CCSS) scheme

###  2.1 Recap of the BFKL evolution

Let us first recap the basics of the BFKL evolution. The BFKL equation has
been first derived in LL approximation in the high energy limit in series of
seminal papers [11, 10, 9]. It can be written generically in the following
form in the momentum space

$f(x,{\bm{k}})=f^{(0)}(x,{\bm{k}})+\int_{x}^{1}\frac{\mathrm{d}z}{z}\int\mathrm{d^{2}}{\bm{k}}\,\mathcal{K}({\bm{k}},{\bm{k}}^{\prime})\,f(\frac{x}{z},{\bm{k}^{\prime}})\;,$
(1)

where $f(x,{{\bm{k}}})$ is the unintegrated gluon density, which depends on
fraction $x$ of the longitudinal momentum and on the transverse momentum
${{\bm{k}}}$ of the reggeized gluon exchanged in the $t$ channel. The BFKL
kernel $\mathcal{K}(\bm{k},\bm{k}_{0})$ has the following perturbative
expansion

$\mathcal{K}(\bm{k},\bm{k}^{\prime})=\overline{\alpha}_{s}(\mu^{2})\mathcal{K}_{0}(\bm{k},\bm{k}^{\prime})\,+\,\overline{\alpha}_{s}^{2}(\mu^{2})\mathcal{K}_{1}(\bm{k},\bm{k}^{\prime})+\,...\;,$
(2)

where the (rescaled) strong coupling constant is defined as

$\displaystyle\bar{\alpha}_{s}(\mu^{2})=\frac{N_{c}}{\pi}\alpha_{s}(\mu^{2})\;,$
(3)

with $N_{c}$ the number of colors and $\mu^{2}$ being the scale of the strong
coupling. In the above equation (2), $\mathcal{K}_{0}$ is LL and
$\mathcal{K}_{1}$ is NLL BFKL kernel respectively. It is worth noting that, to
the NLL accuracy $\mu^{2}$ in principle can take any form in the argument of
the strong coupling the leading kernel, with appropriate subtractions in the
NLL kernel to match the NLL result. On the other hand, to this level of
accuracy the scale in the coupling in the NLL kernel is not determined. This
would require the knowledge of the kernel at NNLL level.

In the following, we shall also use the notation $|{{\bm{k}}}^{2}|=k^{2}$ and
work under the assumption of the angular averaged density, i.e.
$f(x,{\bm{k}})=f(x,k^{2})$. This is justified for the inclusive case, but
would be important for example in the context of the angular distributions in
particle production. It is usually convenient to study the behavior of the
BFKL kernel in the Mellin space. The Mellin transformation is thus introduced
in the following convention,

$\displaystyle\overline{f}(\omega,k^{2})=\int_{0}^{1}\frac{\mathrm{d}x}{x}\,x^{\omega}\,f(x,k^{2})\,,$
(4)
$\displaystyle{\mathcal{F}}(\omega,\gamma)=\int_{0}^{\infty}\mathrm{d}k^{2}\left(k^{2}\right)^{-\gamma}\,\overline{f}(\omega,k^{2})\,.$
(5)

In Mellin space, the BFKL equation can be rewritten into this general
analytical form

${\mathcal{F}}(\omega,\gamma)={\mathcal{F}}^{(0)}(\omega,\gamma)+\frac{\bar{\alpha}_{s}}{\omega}\,\chi(\gamma)\,{\mathcal{F}}(\omega,\gamma)\,.$
(6)

Similar to the Eq. (2), the kernel $\chi(\gamma)$ can be also computed in
different orders of $\overline{\alpha}_{s}$,

$\chi(\gamma)=\bar{\alpha}_{s}\chi_{0}(\gamma)+\bar{\alpha}_{s}^{2}\chi_{1}(\gamma)+....\,.$
(7)

The expressions for the LL and NLL kernels are well known in QCD and they read
[9, 12, 13]

$\chi_{0}(\gamma)=2\psi(1)-\psi(\gamma)-\psi(1-\gamma)\,,$ (8)
$\displaystyle\chi_{1}(\gamma)$ $\displaystyle=$
$\displaystyle\frac{b}{2}\left[\chi_{0}^{2}(\gamma)+\chi_{0}^{\prime}(\gamma)\right]-\frac{1}{4}\chi_{0}^{\prime\prime}(\gamma)-\frac{1}{4}\left(\frac{\pi}{\sin\pi\gamma}\right)^{2}\frac{\cos\pi\gamma}{3(1-2\gamma)}\left[11+\frac{\gamma(1-\gamma)}{(1+2\gamma)(3-2\gamma)}\right]$
(9) $\displaystyle+$
$\displaystyle\left(\frac{67}{36}-\frac{\pi^{2}}{12}\right)\chi_{0}(\gamma)+\frac{3}{2}\xi(3)+\frac{\pi^{2}}{4\sin\pi\gamma}-\Phi(\gamma)\,.$

where $\psi(\gamma)$ is polygamma function, $b=(33-2N_{f})/(12\pi)$ with
$N_{f}$ the number of quarkflavors, and

$\Phi(\gamma)=\sum_{n=0}^{\infty}(-1)^{n}\left[\frac{\psi(n+1+\gamma)-\psi(1)}{(n+\gamma)^{2}}+\frac{\psi(n+2-\gamma)-\psi(1)}{(n+1+\gamma)^{2}}\right].$
(10)

In N=4 sYM theory the kernel $\chi$ has been calculated up to NNLL order [15,
16, 14]. In the leading logarithmic approximation the functions $(k^{\prime
2})^{\gamma}$ are the eigenfunctions of the equation, and $\chi_{0}(\gamma)$
is the eigenvalue at the leading logarithmic order. This is possible since the
coupling is not running and the equation is scale invariant. At NLL, due to
the running of the coupling the conformal invariance of the lowest order
result is broken. The eigenfunctions become scale dependent and include the
running coupling [60]. We note that, nevertheless, it is possible to analyze
this equation in the Mellin space, even with the running coupling, including
higher order corrections.

The expression in Eq. (9) is given for the symmetric scale choice. In order to
transform to the asymmetric scale choice, for example for $s_{0}=k^{2}$ in the
Mellin space one needs to add to Eq. (9), the term

$-\frac{1}{2}\chi_{0}(\gamma)\frac{\partial\chi_{0}}{\partial\gamma}\;,$ (11)

which is the NLL term due to the scale changing transformation, see [13]. At
NNLL the scale changing transformation is much more complicated as it involves
the NLL kernel as well, as was shown in [61, 62].

It has been known since some time that NLL corrections are numerically very
large [12] and lead to the instabilities in the solution, including
possibility of oscillating and negative cross sections [36]. This has also
been seen in the case of the NLL corrections to the non-linear equation, [63],
meaning in particular that the saturation corrections do not reduce the effect
of the NLL terms. The problem has been identified early on, as originating
from the double and triple collinear poles in the Mellin space $\gamma$. These
are poles are coming from the additional transverse logarithms in the momentum
space [36, 37]. The collinear poles, which correspond to the strong ordering
when $k^{2}\gg k^{\prime 2}$ appear as $1/\gamma$ and the anti-collinear poles
which correspond to $k^{2}\ll k^{\prime 2}$ as $1/(1-\gamma)$. At NLL level in
the high energy limit the double poles arise due to the non-singular (in
$1/z$) part of the DGLAP splitting function as well as due to the running of
the coupling. In addition there are also triple collinear poles which appear
due to the kinematical constraints [40, 41]. The fact that the BFKL equation
receives large corrections from the kinematics was identified early on, before
the full NLL calculation. The kinematical constraint [64] was shown to
originate from the improved kinematics, and more precisely from the
requirement that the exchanged momenta are dominated by the transverse
components [65, 66]. There are different forms of the kinematical constraint
available in the literature, however they all resum the triple collinear poles
at NLL and higher at more subleading orders. For more recent work on different
forms of kinematical constraint see [62]. There it has been shown that
kinematical constraint is consistent with the NNLL calculation in N=4 sYM, and
the resummation of the $1/\gamma^{5}$ poles was demonstrated.

###  2.2 The CCSS resummed equation

The CCSS resummation scheme was formulated in [40], where more details can be
found. Here we only summarize main points and features of the construction,
essential for the subsequent discussion. The resummation was motivated in the
Mellin space, where the analysis of the collinear poles is more transparent.
The final formulation however was constructed in the momentum space, and the
integral equation for the gluon Green’s function solved, see [40] and more
recently [67]. One of the advantages of the formulating the resummation in the
momentum space in the form of the integral equation, is the fact that the
incorporation of the running of the coupling is more straightforward. Also,
one can then solve directly for the unintegrated gluon density, which can then
be used for other phenomenological applications. This is the approach which we
shall adopt in this work too. A similar idea for the resummation was
formulated previously in Ref. [58] by combining the DGLAP and BFKL evolution
with the kinematical constraint. The latter approach was formulated without
full NLL BFKL correction. A fit in the same spirit was performed in [59],
together with the nonlinear term, and we shall perform more detailed
comparison of the results of the CCSS resummation fits with that approach in
Sec. 4.

In the CCSS scheme the starting point is the fixed order LL + NLL calculation.
On top of that, the kinematical constraint as well as the non-singular part of
the DGLAP splitting function are added. Next, appropriate subtractions are
made from the NLL kernel, which remove the double and triple poles and which
are already incorporated in the kinematical constraints and DGLAP terms. This
is to avoid the double counting. Finally, the $\beta_{0}$ function dependent
term is subtracted from the NLL part and included in the running of the
coupling in front of the leading logarithmic kernel. The scale was taken equal
to be the transverse momentum of the emitted gluon. In that way there are no
terms proportional to the beta function in the NLL kernel. The expression for
the NLL part of the resummed kernel in the Mellin space is then

$\chi_{1}^{\rm
subtr}(\gamma)\;=\;\chi_{1}(\gamma)+\frac{1}{2}\chi_{0}(\gamma)\frac{\pi^{2}}{\sin^{2}(\pi\gamma)}-\chi_{0}(\gamma)\frac{A_{1}(0)}{\gamma(1-\gamma)}+\frac{b}{2}(\chi_{0}^{\prime}+\chi_{0}^{2})\;.$
(12)

The first term is the fixed order NLL term, the same as in Eq.(9). The next
three terms are the subtractions: due to the kinematical constraint, DGLAP
terms and the running coupling. In addition, more subtractions were added to
avoid the double counting and ensure the conservation of the momentum sum
rule. Certain flexibility as to how to incorporate such subtractions leads to
the various resummation schemes.

As already mentioned above, the choice of the scale for the coupling constant
in the NLL kernel is arbitrary. This is because any change at this level would
require the knowledge of the NNLL term, which is currently not known in QCD.
In the CCSS scheme [40] the scale at NLL was chosen to be the maximum of the
momenta squared of the exchanged gluons, i.e. $\max(k^{2},k^{\prime 2})$. We
will also adopt this choice for our calculation.

The final result of the CCSS resummation was formulated directly in the
momentum space through integral equation. As mentioned above, this allows for
more control over the implementation of the running coupling corrections.
Together with the DGLAP splitting functions this equation becomes an integral
equation in both the longitudinal $z$ and transverse momentum components
${{{\bm{k}}}}$. The three main contributions

${\mathcal{K}}_{0}^{\rm
kc}(z;{{\bm{k}}},{{\bm{k}}}^{\prime})\stackrel{{\scriptstyle
z,{\bm{q}}}}{{\otimes}}f(\frac{x}{z},k^{\prime})+{\mathcal{K}}_{0}^{\rm
coll}(z;k,k^{\prime})\stackrel{{\scriptstyle
z,k^{\prime}}}{{\otimes}}f(\frac{x}{z},k^{\prime})+{\mathcal{K}}_{1}^{\rm
subtr}(z;k,k^{\prime})\stackrel{{\scriptstyle
z,k^{\prime}}}{{\otimes}}f(\frac{x}{z},k^{\prime})\;,$ (13)

are coming from the leading logarithmic kernel with kinematical constraint,
the so-called collinear (DGLAP) part and the NLL part with subtractions.

The first term in Eq. (13) is

${\mathcal{K}}_{0}^{\rm
kc}(z;{{\bm{k}}},{{\bm{k}}}^{\prime})\stackrel{{\scriptstyle
z,{\bm{q}}}}{{\otimes}}f(\frac{x}{z},k^{\prime})\\\
=\int_{x}^{1}\frac{dz}{z}\int\frac{d^{2}{\bm{q}}}{\pi{\bm{q}}^{2}}\;\bar{\alpha}_{s}({\bm{q}}^{2})\left[f(\frac{x}{z},|{{\bm{k}}}+{\bm{q}}|)\Theta(\frac{k^{2}}{z}-k^{\prime
2})-\Theta(k-q)f(\frac{x}{z},k)\right]\;,$ (14)

where ${\bm{q}}={{\bm{k}}}-{{\bm{k}}}^{\prime}$ corresponds to the transverse
momentum of the emitted gluon. This choice of a scale in the running coupling
is convenient since in this case the $b$ dependent terms in the NLL part of
the kernel are exactly zero. The kinematical (or consistency) constraint is
implemented onto the real emissions only (see discussion in [66] and [67]). It
is here asymmetric, which corresponds to the asymmetric scale choice suitable
for the DIS problem we are considering. It is implemented as

$k^{\prime 2}\leq\frac{k^{2}}{z}\;.$ (15)

To be precise, there are different versions of this constraint which appear in
the literature, see [65, 66, 64]. Detailed analysis (see [62]) showed that all
versions are generating the same leading $1/\gamma^{3}$ poles in the Mellin
space at NLL, $1/\gamma^{5}$ poles in NNLL level (for supersymmetric case) and
they do not generate any double poles, and with the difference starting to
appear in the single pole level.

The second contribution in (13) is

${\mathcal{K}}_{0}^{\rm coll}(z;k,k^{\prime})\stackrel{{\scriptstyle
z,k^{\prime}}}{{\otimes}}f(\frac{x}{z},k^{\prime})=\int_{x}^{1}{dz\over
z}\int_{0}^{k^{2}}\frac{{dk^{\prime}}^{2}}{k^{2}}\;\bar{\alpha}_{s}(k^{2})z\tilde{P}_{gg}(z)f(\frac{x}{z},k^{\prime})\\\
+\int_{x}^{1}{dz\over
z}\int_{k^{2}}^{k^{2}/z}\frac{{dk^{\prime}}^{2}}{k^{\prime}{}^{2}}\;\bar{\alpha}_{s}({k^{\prime}}^{2})z{{k^{\prime
2}}\over k^{2}}\tilde{P}_{gg}(z{k^{\prime 2}\over
k^{2}})f(\frac{x}{z},k^{\prime})\;.$ (16)

It is the sum of the collinear and anticollinear parts with the non-singular
part of the splitting function

$\tilde{P}^{(0)}_{gg}=P^{(0)}_{gg}-{1\over z}\;,$ (17)

where the $P^{(0)}_{gg}$ is the DGLAP gluon-gluon splitting function in LO.

Finally, the last term in Eq.(13) is the NLL part of the BFKL with appropriate
subtractions (corresponding to expression in Eq. (12)) transformed into
momentum space

$\displaystyle\int_{x}^{1}\frac{dz}{z}$ $\displaystyle\int
dk^{\prime}{}^{2}\;\bar{\alpha}^{2}_{s}({k}^{2}_{>})\tilde{K}_{1}(k,k^{\prime})f(\frac{x}{z},k^{\prime})\;=\;{1\over
4}\int_{x}^{1}\frac{dz}{z}\int
d{k^{\prime}}^{2}\;\bar{\alpha}^{2}_{s}({k}^{2}_{>})\bigg{\\{}$
$\displaystyle{\left({67\over 9}-{\pi^{2}\over
3}\right){1\over|{k^{\prime}}^{2}-k^{2}|}\left[f(\frac{x}{z},{k^{\prime}}^{2})-{2k_{<}^{2}\over({k^{\prime}}^{2}+k^{2})}f(\frac{x}{z},k^{2})\right]+}$
$\displaystyle{\bigg{[}-{1\over 32}\left({2\over{k^{\prime}}^{2}}+{2\over
k^{2}}+\left({1\over{k^{\prime}}^{2}}-{1\over
k^{2}}\right)\log\left({k^{2}\over{k^{\prime}}^{2}}\right)\right)+{4{\mathrm{Li}_{2}}(1-k_{<}^{2}/k_{>}^{2})\over|{k^{\prime}}^{2}-k^{2}|}}$
$\displaystyle{-4A_{1}(0){\rm sgn}({k}^{2}-{k^{\prime}}^{2})\left({1\over
k^{2}}\log{|{k^{\prime}}^{2}-k^{2}|\over{k^{\prime}}^{2}}-{1\over{k^{\prime}}^{2}}\log{|{k^{\prime}}^{2}-k^{2}|\over{k}^{2}}\right)}$
$\displaystyle-\left(3+\left({3\over 4}-{({k^{\prime}}^{2}+k^{2})^{2}\over
32{k^{\prime}}^{2}k^{2}}\right)\right)\int_{0}^{\infty}{dy\over
k^{2}+y^{2}{k^{\prime}}^{2}}\log|{1+y\over 1-y}|$
$\displaystyle+{1\over{k^{\prime}}^{2}+k^{2}}\left({\pi^{2}\over
3}+4{\mathrm{Li}_{2}}({k_{<}^{2}\over
k_{>}^{2}})\right)\bigg{]}f(\frac{x}{z},k^{\prime})\bigg{\\}}$
$\displaystyle+{1\over
4}6\zeta(3)\int_{x}^{1}\frac{dz}{z}\;\bar{\alpha}^{2}_{s}(k^{2})f(\frac{x}{z},k)\;.\hskip
113.81102pt$ (18)

The above construction for the resummed kernel needs to be supplemented by
additional subtractions. It turns out, [40] that there are terms which are
giving spurious DGLAP anomalous dimension at NLO. This needs to be canceled by
appropriate subtraction and it was achieved by adding extra terms to the
kernel. Obviously, there is some ambiguity in this procedure, since one is
working with the information up to a fixed order in perturbation theory.
Therefore two different schemes were proposed $A$,$B$ in [40]. In the
following, we shall utilize scheme $B$ from that work.

##  3 Contributions to structure function

The structure function $F_{2}$ can be evaluated by the $k_{T}$ factorization
theorem, which involves an off-shell matrix element and the unintegrated gluon
density. The structure function $F_{2}$ receives however large contributions
from the non-perturbative, or soft, regime. This is parametrized in our
description as the contribution coming from the low momenta of the gluon
$k^{2}$ and with the addition of the soft Pomeron contribution. The setup is
similar to the one presented in [66], without however the matrix formulation
which would involve the evolution of quarks.

###  3.1 Perturbative contribution

The perturbative contribution to the structure function is based on the
$k_{T}$ factorization theorem, together with the unintegrated gluon density
obtained from the CCSS resummed evolution equation discussed in the previous
section. The expression for the structure function $F_{2}$ from the $k_{T}$
factorization is given by

$F_{2}(x,Q^{2})=\sum_{q}\,e_{q}^{2}\,S_{q}(x,Q^{2})\;,$ (19)

where the sum is over the quark flavors and general expression for
$S_{q}(x,Q^{2})$ is

$S_{q}(x,Q^{2})\;=\;\int_{x}^{1}\>\frac{dz}{z}\>\int\>\frac{dk^{2}}{k^{2}}\>S_{\rm
box}^{q}\>(z,m_{q}^{2},k^{2},Q^{2})\>f\left(\frac{x}{z},k^{2}\right)\;.$ (20)

The explicit expression for the convolution of the matrix element with the
unintegrated gluon density is given by [68, 58]

$\displaystyle S_{q}(x,Q^{2})$ $\displaystyle=$
$\displaystyle\frac{Q^{2}}{4\pi^{2}}\int\frac{dk^{2}}{k^{4}}\int_{0}^{1}d\beta\int
d\kappa^{\prime}\alpha_{s}\left\\{\left[\beta^{2}+(1-\beta^{2})\right]\left(\frac{\bm{\kappa}}{D_{1q}}-\frac{\bm{\kappa}-\bm{k}}{D_{2q}}\right)^{2}\right.$
$\displaystyle+$
$\displaystyle\left.\left[m_{q}^{2}+4Q^{2}\beta^{2}(1-\beta)^{2}\right]\left(\frac{1}{D_{1q}}-\frac{1}{D_{2q}}\right)^{2}\right\\}f\left(\frac{x}{z},k^{2}\right)\Theta\left(1-\frac{x}{z}\right).$
(21)

In the above, $\bm{\kappa}$ and $\bm{k}$ are quark and gluon transverse
momenta respectively, and $\beta$ is the variable defined in the Sudakov
decomposition of the quark momentum (longitudinal momentum fraction of the
photon carried by the quark, for details see Ref.[68]). In addition it is
useful to defined the shifted quark transverse momentum is
$\bm{\kappa}^{\prime}=\bm{\kappa}-(1-\beta)\bm{k}$. The energy denominators
are

$\displaystyle D_{1q}$ $\displaystyle=$
$\displaystyle\kappa^{2}+\beta(1-\beta)Q^{2}+m_{q}^{2}\;,$ (22) $\displaystyle
D_{2q}$ $\displaystyle=$
$\displaystyle(\bm{\kappa}-\bm{k})^{2}+\beta(1-\beta)Q^{2}+m_{q}^{2}\;.$ (23)

The argument of the unintegrated gluon density is equal to $x/z$ with

$z=\left[1+\frac{\kappa^{\prime
2}+m_{q}^{2}}{\beta(1-\beta)Q^{2}}+\frac{k^{2}}{Q^{2}}\right]^{-1}\;.$ (24)

This stems from the exact kinematics in the photon-gluon fusion process, see
[68]. As analyzed in detail in [69, 70] the exact kinematics in the impact
factor, goes beyond the leading order approximation in high energy. It has
been demonstrated that it leads to large effect numerically and is important
for phenomenology [71].

The argument of the strong coupling $\alpha_{s}$ is taken to be
$(k^{2}+\kappa^{2}+m_{q}^{2})$ in this analysis. The masses of quarks are
taken to be $m_{u}=m_{d}=m_{s}=0$ and $m_{c}=1.4\ \text{GeV}$. The integration
over the transverse momenta in the $k_{T}$ factorization formula formally
extends down to zero into the non-perturbative region. We assume the validity
of the formula (20) and (3.1) only for the transverse momenta
$k^{2},\kappa^{2}>k_{0}^{2}$ where cutoff $k_{0}^{2}$ parametrizes the
boundary between the perturbative and non-perturbative regions of the
transverse momentum. We took the value of $k_{0}^{2}=1.\,\rm GeV^{2}$ for the
cutoff.

###  3.2 Non-perturbative contribution

The structure function $F_{2}$ receives large soft contribution. For example
in the approach of [72], it has been simply parametrized as the constant
background term in addition to the perturbative small $x$ part. In the
approaches within the dipole model, the non-perturbative contribution is
usually taken automatically into account by integration over the large dipole
sizes with the flat dipole cross section, see for example discussion in [73,
74]. Here we follow the approach of [58] where the non-perturbative
contribution from the low gluon and quark transverse momenta can be
parametrized as follows. In the non-perturbative region, where both quark
momenta and gluon momenta are small $k^{2},\kappa^{\prime 2}<k_{0}^{2}$, we
assume that light quark contribution is phenomenologically evaluated as the
soft Pomeron exchange [75]. That is we assume the soft Pomeron form for
$u,d,s$ contributions

$S^{(a)}\;=\;S_{u}^{I\\!\\!P}\>+\>S_{d}^{I\\!\\!P}\>+\>S_{s}^{I\\!\\!P}\;,$
(25)

and

$S_{u}^{I\\!\\!P}\;=\;S_{d}^{I\\!\\!P}\;=\;2S_{s}^{I\\!\\!P}\;=\;C_{I\\!\\!P}\>x^{-\lambda}\>(1-x)^{8},$
(26)

where coefficient $C_{I\\!\\!P}$ is a free parameter independent of $Q^{2}$.
Here, $0\lesssim\lambda\lesssim 0.1$ is the soft-Pomeron power.

In the region where the quark momenta are higher but the gluon momenta are
very low, i.e. $k^{2}<k_{0}^{2}<\kappa^{\prime 2}$, we assume the strong-
ordered approximation at quark-gluon vertex and used the the collinear
approximation

$S_{\rm box}\;\rightarrow\;S_{\rm box}^{(b)}\>(z,k^{2}=0,Q^{2}),$ (27)

and thus we have

$\displaystyle S^{(b)}$ $\displaystyle=$
$\displaystyle\int_{x}^{1}\>\frac{\mathrm{d}z}{z}\>S_{\rm
box}^{(b)}\>(z,k^{2}=0,Q^{2})\>\int_{0}^{k_{0}^{2}}\>\frac{\mathrm{d}k^{2}}{k^{2}}\>f\left(\frac{x}{z},k^{2}\right)$
(28) $\displaystyle=$
$\displaystyle\int_{x}^{1}\>\frac{\mathrm{d}z}{z}\>S_{\rm
box}^{(b)}\>(z,k^{2}=0,Q^{2})\>\frac{x}{z}\>g\left(\frac{x}{z},k_{0}^{2}\right)\;,$

where $xg\left(x,k_{0}^{2}\right)$ is the non-perturbative input collinear
gluon density at scale $k_{0}^{2}$, and the form of its parametrization will
be specified in Sec. 4.

Thus the complete contribution from the light quarks is simply the sum of

$S^{(a)}_{q}+S^{(b)}_{q}+S^{(c)}_{q}\;,$ (29)

where the last, perturbative contribution is computed from Eq. (3.1) with the
lower bound on the transverse momentum integration given by the cutoff
$k_{0}^{2}$.

###  3.3 Charm quark contribution

In addition to the light quarks, one needs to include the charm quark
contribution, which is also evaluated from the $k_{T}$ factorization. The mass
of the charm quark has been taken to be $m_{c}=1.4\;\rm GeV$.

For the gluon transverse momenta larger than the cutoff $k^{2}>k_{0}^{2}$, we
treat the charm contributions perturbatively, using Eq.(3.1). In the region
$k^{2}<k_{0}^{2}$, we use on-shell approximation, i.e. $k^{2}=0$. Therefore
the contribution is then

$\displaystyle S^{(b)}_{c}$ $\displaystyle=$
$\displaystyle\int_{x}^{a}\frac{\mathrm{d}z}{z}S_{\text{box}}(z,k^{2}=0,Q^{2};m_{c}^{2})\int_{0}^{k_{0}^{2}}\frac{\mathrm{d}k^{2}}{k^{2}}f\left(\frac{x}{z},k_{0}^{2}\right)$
(30) $\displaystyle=$
$\displaystyle\int_{x}^{a}\frac{\mathrm{d}z}{z}S_{\text{box}}(z,k^{2}=0,Q^{2};m_{c}^{2})\frac{x}{z}g\left(\frac{x}{z},k_{0}^{2}\right)\;,$

where $a=(1+4m_{c}^{2}/Q^{2})^{-1}$ and from Ref. [76]

$\displaystyle S_{\text{box}}(z,k^{2}=0,Q^{2};m_{c}^{2})$ $\displaystyle=$
$\displaystyle\frac{z\alpha_{s}}{2\pi}\left\\{\ln\frac{1+\xi}{1-\xi}\left[z^{2}+(1-z)^{2}+z(1-3z)\frac{4m_{c}^{2}}{Q^{2}}-z^{2}\frac{8m_{c}^{4}}{Q^{4}}\right]\right.$
(31)
$\displaystyle\left.+\xi\left[-1+8z(1-z)-z(1-z)\frac{4m_{c}^{2}}{Q^{2}}\right]\right\\}.$

Here, $\xi^{2}=1-(4m_{c}^{2}/Q^{2})z(1-z)^{-1}$. These are two contributions
from the charm quark that are included, which are dynamically generated from
the the photon-gluon fusion. We do not consider any additional ‘intrinsic’
charm contribution.

##  4 Numerical results

In this section we present the results of the fits to the structure functions
and the properties of the extracted unintegrated gluon distribution function.
The fits have been performed to the DIS HERA $F_{2}$[77] and charm data
$F_{2}^{c}$ [78].

In this analysis we focus on high energy, or low-$x$ physics, therefore we
impose the cuts on the data: $Q^{2}>2\ \text{GeV}^{2},\ x<0.01$. Since we are
working with the data on the reduced cross section, in order to exclude the
contribution from the longitudinal structure function $F_{L}$, we also impose
the cut on inelasticity $y<0.6$.

The fitted parameters stem from the input gluon distribution and the soft
Pomeron contribution. The initial condition for the evolution of the
unintegrated gluon density from the CCSS evolution is given by the convolution
of the integrated gluon density with the DGLAP splitting function as in [58]

$\tilde{f}^{(0)}(x,k^{2})\;=\;\frac{\alpha_{S}(k^{2})}{2\pi}\int_{x}^{1}dzP_{gg}(z)\frac{x}{z}g\left(\frac{x}{z},k_{0}^{2}\right)\;,$
(32)

where $xg\left(x,k_{0}^{2}\right)$ is the integrated gluon density at scale
$k_{0}^{2}$. This is done to avoid explicit parametrization of the initial
term for the unintegrated gluon density in the non-perturbative region
$k^{2}<k_{0}^{2}$. That is the input in the resummed renormalization group
improved equation is of the same form as in the conventional DGLAP evolution.
In the case of the unified DGLAP and BFKL evolution it proved to be a very
successful approach and the input could be parametrized in a very economical
form with only few free parameters [58]. For the integrated input gluon
distribution we consider the form

$xg(x,k_{0}^{2})=N(1-x)^{\beta}\left[1+D\left(x+\epsilon\right)^{\alpha}\right]\;,$
(33)

where $\epsilon$ is manually set to be a small positive number, as such we can
prevent a potential negativity of the gluon input function when $D<0$ and
$\alpha<0$.

We have assumed two different fitting scenarios. In the first one, the fit has
been performed to both $F_{2}$ and $F_{2}^{c}$ data simultaneously. In the
second scenario, only $F_{2}$ data have been fitted, thus leaving the charm
structure function $F_{2}^{c}$ as a prediction. The resulting parameters of
the fit, which include the parameters on the initial gluon distribution and
the parameters of the soft Pomeron part, (26), are shown in Table 1. We
observe that the quality of the fits is very good and similar in both cases,
as are also the resulting values of the parameters, which demonstrates the
stability of the approach.

Table 1: Fitting parameters for the gluon distribution and the soft Pomeron part. Data Range | $\chi^{2}$ | $C_{I\\!\\!P}$ | $\lambda$ | $N$ | $\beta$ | $\alpha$ | $D$ | $\epsilon$
---|---|---|---|---|---|---|---|---
Fit $F_{2}$,$F_{2}^{c}$ | 0.9900 | 0.4420 | 0.01562 | 3.894 | 4.951 | -0.4402 | -0.1096 | 0.0003
Fit $F_{2}$ | 1.052 | 0.4427 | 0.01557 | 3.887 | 4.950 | -0.4413 | -0.1084 | 0.0003

In the following we thus only show the results obtained when the charm
structure function is also fitted. In Fig. 1 we show the structure function
$F_{2}(x,Q^{2})$ plotted as a function of $x$ for selected values of
$Q^{2}=2,15,35,90,150,250\;\rm GeV^{2}$. We see that the calculation describes
the experimental data on structure function $F_{2}$ very well. In Fig. 2 the
calculation is shown for the charm structure function $F_{2}^{c}$. The red
points indicate the charm component extracted from experiment assuming
extrapolating factors to the full phase space in the HVQDIS scenario and black
points in the CASCADE scheme, for details see [78]. We see again that the
description of the experimental data on charm structure function is very good.
As mentioned earlier, the quality of the description is comparable to the case
when charm is not fitted but is a prediction.

In Figs. 3,4,5 we show the structure function broken down into separate
components $F_{2}^{(a)},F_{2}^{(b)},F_{2}^{(c)}$ as a function of $x$ for
three values of $Q^{2}=150,15,5\;\rm GeV^{2}$ respectively.

The perturbative component, indicated as $F_{2}^{(c)}$ in figures 4 and 3
dominates at large values of $Q^{2}$ and small $x$. The two non-perturbative
components have very flat dependence on $x$ and contribute to the most of the
cross section at moderate $x$ and lower $Q^{2}$. We see that the non-
perturbative contribution due to the soft Pomeron is substantial. For example,
it is about $25\%$ at small $x\simeq 5\times 10^{-4}$ and $Q^{2}=15\;\rm
GeV^{2}$. Even at high $Q^{2}$ the soft component is still non-negligible. In
the low $Q^{2}$ region, see Fig. 5, the non-perturbative components dominate
the structure function in a very wide range of $x$, down to $x\sim 3\times
10^{-4}$. It is only at smallest $x$ where the perturbative component starts
to dominate at low $Q^{2}=5\,\rm GeV^{2}$.

In Figs. 6 and 7 we show the unintegrated gluon distribution for two fixed
values of $x=0.01,0.001$ as a function of transverse momentum squared $k^{2}$.
The present calculation based on the CCSS resummation is compared with the
extracted gluon distribution from [59] (called KS for short here) in the
linear and non-linear case. The KS calculation is similar in its philosophy to
the one presented here and is widely used in the phenomenology, and thus it is
interesting to compare it with the current extraction. The KS calculation is
based on the LL BFKL with kinematical constraint and includes DGLAP terms, so
the main difference is the NLL part of the BFKL and the type of subtractions
performed. The overall shape in $k^{2}$ is similar for the CCSS and KS
calculation.

In Figs. 8-10 the unintegrated gluon density is shown as a function of $x$ for
fixed values of $k^{2}=10,100,1000\;\rm GeV^{2}$, again compared with the KS-
linear and KS-nonlinear calculation. Unsurprisingly, the small $x$ behavior of
the CCSS calculation is very close to the KS-linear, whereas both calculations
differ substantially from the nonlinear calculation at low $x$ and low $Q^{2}$
where the saturation corrections are the strongest. There are however some
subtle differences between the KS-linear and CCSS calculations which indicate
somewhat stronger small $x$ behavior in the CCSS scenario.

The differences between linear and non-linear scenarios are also visualised by
taking the ratios between the calculations shown in Fig. 11, 12, 13. To be
precise the ratios CCSS/KS-linear and CCSS/KS-nonlinear are shown as a
function of $x$ for three different values of the momentum squared
$k^{2}=10,100,1000\;\rm GeV^{2}$. We see that the CCSS/KS-linear ratio is
close to unity, for most values of $x$, whereas the ratio to the calculation
which includes the nonlinear effects deviates substantially from unity at low
$x$ and low $k^{2}$. Still we see again marked difference, the CCSS
calculation tends to have slightly faster rise towards small $x$. The small
$x$ behavior is also illustrated in Figs. 14, 15 and 16, where the effective
Pomeron intercept is shown by performing the logarithmic derivative of the
unintegrated gluon distribution

$\lambda_{\rm eff}=\frac{\partial\ln f(x,k^{2})}{\partial\ln 1/x}\;.$ (34)

It is seen from all plots that the effective power is very close, especially
at low $x$ to the linear KS calculation, which is expected, and differs from
the non-linear KS, for low $x$ and lowest value of $k^{2}$. We also see that
the power is larger for CCSS calculation than for the linear KS calculation,
which is consistent with previous observations. Asymptotically, for low $x$
and low to moderate $k^{2}$ the value of the power approaches about $0.3$
which is the power observed in HERA data.

## Conclusions and outlook

In this work we have used the renormalization group improved small $x$
evolution equation, which resums important terms in the low $x$ region, and
used it for the first time to describe the structure function $F_{2}$ in DIS.
Very good description of the experimental data in the low $x<0.01$ region was
obtained, with small number of parameters. In addition, the charm structure
function was calculated, in two scenarios, where the charm data were used for
the fit, and in the second scenario where the structure function $F_{2}^{c}$
was not fitted but obtained instead as a prediction. In both cases the
description of the experimental data was very good and the values of the
parameters were very similar. The fit was based on the $k_{T}$ factorization
formula with the decomposition of the integral in transverse momenta into soft
and hard regimes. The unintegrated gluon distribution extracted from this
analysis was compared with the previous calculations, based on the similar
approach and it was found to be consistent, though with some marked
differences. The unintegrated gluon density has somewhat higher power which
governs the low $x$ growth as compared to the resummed schemes, which were
based on the LL BFKL with kinematical constraint and DGLAP but without the
full NLL BFKL terms. This growth should be tamed eventually by including the
non-linear corrections from gluon saturation.

For the outlook, the current analysis should be potentially expanded to
include the DGLAP splitting function up to NNLO, as well as the impact factor
up to NLO [79]. Current approach includes LO impact factor, with exact
kinematics, which go beyond the lowest order calculation. In addition, the
formulation should include the quarks in the evolution, which could be done by
employing the matrix approach to small $x$ evolution as was suggested in [80].
We leave these outstanding problems to the future work.

## Acknowledgments

This work was supported in part by the U.S. Department of Energy grant No. DE-
SC-0002145 and in part by National Science Centre in Poland, grant
2019/33/B/ST2/02588. We thank Dimitri Colferai and Krzysztof Kutak for
discussions and comments.

## References

  * [1] M. Breidenbach, J. I. Friedman, H. W. Kendall, E. D. Bloom, D. H. Coward, H. C. DeStaebler, J. Drees, L. W. Mo, and R. E. Taylor, “Observed behavior of highly inelastic electron-proton scattering,” Phys. Rev. Lett., vol. 23, pp. 935–939, 1969.
  * [2] I. Abt et al., “Measurement of the proton structure function $F_{2}(x,Q^{2})$ in the low x region at HERA,” Nucl. Phys. B, vol. 407, pp. 515–538, 1993.
  * [3] M. Derrick et al., “Measurement of the proton structure function $F_{2}$ in ep scattering at HERA,” Phys. Lett. B, vol. 316, pp. 412–426, 1993\.
  * [4] G. Altarelli and G. Parisi, “Asymptotic Freedom in Parton Language,” Nucl. Phys., vol. B126, pp. 298–318, 1977.
  * [5] Y. L. Dokshitzer, “Calculation of the Structure Functions for Deep Inelastic Scattering and $e^{+}e^{-}$ Annihilation by Perturbation Theory in Quantum Chromodynamics.,” Sov. Phys. JETP, vol. 46, pp. 641–653, 1977.  [Zh. Eksp. Teor. Fiz.73,1216(1977)].
  * [6] V. N. Gribov and L. N. Lipatov, “Deep inelastic ep scattering in perturbation theory,” Sov. J. Nucl. Phys., vol. 15, pp. 438–450, 1972.  [Yad. Fiz.15,781(1972)].
  * [7] V. N. Gribov and L. N. Lipatov, “$e^{+}e^{-}$ pair annihilation and deep inelastic ep scattering in perturbation theory,” Sov. J. Nucl. Phys., vol. 15, pp. 675–684, 1972.  [Yad. Fiz.15,1218(1972)].
  * [8] J. C. Collins, D. E. Soper, and G. F. Sterman, “Factorization of Hard Processes in QCD,” Adv. Ser. Direct. High Energy Phys., vol. 5, pp. 1–91, 1989.
  * [9] L. N. Lipatov, “The Bare Pomeron in Quantum Chromodynamics,” Sov. Phys. JETP, vol. 63, pp. 904–912, 1986.  [Zh. Eksp. Teor. Fiz.90,1536(1986)].
  * [10] E. A. Kuraev, L. N. Lipatov, and R. S. Fadin, “The Pomeranchuk Singularity in Nonabelian Gauge Theories,” Sov. Phys. JETP, vol. 45, pp. 199–204, 1977\.  [Zh. Eksp. Teor. Fiz.72,377(1977)].
  * [11] I. I. Balitsky and L. N. Lipatov, “The Pomeranchuk Singularity in Quantum Chromodynamics,” Sov. J. Nucl. Phys., vol. 28, pp. 822–829, 1978.  [Yad. Fiz.28,1597(1978)].
  * [12] V. S. Fadin and L. N. Lipatov, “BFKL pomeron in the next-to-leading approximation,” Phys. Lett. B, vol. 429, pp. 127–134, 1998.
  * [13] M. Ciafaloni and G. Camici, “Energy scale(s) and next-to-leading BFKL equation,” Phys. Lett. B, vol. 430, pp. 349–354, 1998.
  * [14] S. Caron-Huot and M. Herranen, “High-energy evolution to three loops,” Journal of High Energy Physics, vol. 2018, no. 2, p. 58, 2018.
  * [15] V. N. Velizhanin, “BFKL pomeron in the next-to-next-to-leading approximation in the planar N=4 SYM theory,” 2015.
  * [16] N. Gromov, F. Levkovich-Maslyuk, and G. Sizov, “Pomeron Eigenvalue at Three Loops in $\mathcal{N}=$ 4 Supersymmetric Yang-Mills Theory,” Phys. Rev. Lett., vol. 115, no. 25, p. 251601, 2015.
  * [17] I. Balitsky, “Operator expansion for high-energy scattering,” Nucl. Phys. B, vol. 463, pp. 99–160, 1996.
  * [18] I. Balitsky, “Factorization and high-energy effective action,” Phys. Rev. D, vol. 60, p. 014020, 1999.
  * [19] I. Balitsky and G. A. Chirilli, “Next-to-leading order evolution of color dipoles,” Phys. Rev., vol. D77, p. 014019, 2008.
  * [20] A. H. Mueller, “Small x Behavior and Parton Saturation: A QCD Model,” Nucl. Phys. B, vol. 335, pp. 115–137, 1990.
  * [21] A. H. Mueller, “Soft gluons in the infinite momentum wave function and the BFKL pomeron,” Nucl. Phys. B, vol. 415, pp. 373–385, 1994.
  * [22] Y. V. Kovchegov, “Unitarization of the BFKL pomeron on a nucleus,” Phys. Rev. D, vol. 61, p. 074018, 2000.
  * [23] Y. V. Kovchegov, “Small $x$ $F_{2}$ structure function of a nucleus including multiple pomeron exchanges,” Phys. Rev. D, vol. 60, p. 034008, 1999.
  * [24] L. D. McLerran and R. Venugopalan, “Gluon distribution functions for very large nuclei at small transverse momentum,” Phys. Rev. D, vol. 49, pp. 3352–3355, 1994.
  * [25] L. D. McLerran and R. Venugopalan, “Computing quark and gluon distribution functions for very large nuclei,” Phys. Rev. D, vol. 49, pp. 2233–2241, 1994.
  * [26] L. D. McLerran and R. Venugopalan, “Green’s functions in the color field of a large nucleus,” Phys. Rev. D, vol. 50, pp. 2225–2233, 1994.
  * [27] J. Jalilian-Marian, A. Kovner, L. D. McLerran, and H. Weigert, “The Intrinsic glue distribution at very small x,” Phys. Rev. D, vol. 55, pp. 5414–5428, 1997.
  * [28] J. Jalilian-Marian, A. Kovner, and H. Weigert, “The Wilson renormalization group for low x physics: Gluon evolution at finite parton density,” Phys. Rev. D, vol. 59, p. 014015, 1998.
  * [29] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, “The Wilson renormalization group for low x physics: Towards the high density regime,” Phys. Rev. D, vol. 59, p. 014014, 1998.
  * [30] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, “The BFKL equation from the Wilson renormalization group,” Nucl. Phys. B, vol. 504, pp. 415–431, 1997.
  * [31] E. Iancu, A. Leonidov, and L. D. McLerran, “Nonlinear gluon evolution in the color glass condensate. 1.,” Nucl. Phys. A, vol. 692, pp. 583–645, 2001\.
  * [32] E. Iancu, A. Leonidov, and L. D. McLerran, “The Renormalization group equation for the color glass condensate,” Phys. Lett. B, vol. 510, pp. 133–144, 2001.
  * [33] A. Kovner, M. Lublinsky, and Y. Mulian, “NLO JIMWLK evolution unabridged,” JHEP, vol. 08, p. 114, 2014.
  * [34] M. Lublinsky and Y. Mulian, “High Energy QCD at NLO: from light-cone wave function to JIMWLK evolution,” JHEP, vol. 05, p. 097, 2017.
  * [35] I. Bojak and M. Ernst, “Small x resummations confronted with $F_{2}(x,Q^{2})$ data,” Phys. Lett. B, vol. 397, pp. 296–304, 1997.
  * [36] G. P. Salam, “A Resummation of large subleading corrections at small x,” JHEP, vol. 07, p. 019, 1998.
  * [37] G. P. Salam, “An Introduction to leading and next-to-leading BFKL,” Acta Phys. Polon., vol. B30, pp. 3679–3705, 1999.
  * [38] M. Ciafaloni, D. Colferai, and G. P. Salam, “A collinear model for small x physics,” JHEP, vol. 10, p. 017, 1999.
  * [39] M. Ciafaloni, D. Colferai, and G. P. Salam, “Renormalization group improved small x equation,” Phys. Rev., vol. D60, p. 114036, 1999.
  * [40] M. Ciafaloni, D. Colferai, D. Colferai, G. P. Salam, and A. M. Stasto, “Extending QCD perturbation theory to higher energies,” Phys. Lett. B, vol. 576, pp. 143–151, 2003.
  * [41] M. Ciafaloni, D. Colferai, G. P. Salam, and A. M. Stasto, “Renormalization group improved small x Green’s function,” Phys. Rev. D, vol. 68, p. 114003, 2003.
  * [42] M. Ciafaloni, D. Colferai, G. P. Salam, and A. M. Stasto, “The Gluon splitting function at moderately small x,” Phys. Lett., vol. B587, pp. 87–94, 2004.
  * [43] R. S. Thorne, “The Running coupling BFKL anomalous dimensions and splitting functions,” Phys. Rev., vol. D64, p. 074005, 2001.
  * [44] C. D. White and R. S. Thorne, “A Global Fit to Scattering Data with NLL BFKL Resummations,” Phys. Rev., vol. D75, p. 034005, 2007.
  * [45] G. Altarelli, R. D. Ball, and S. Forte, “Resummation of singlet parton evolution at small x,” Nucl. Phys., vol. B575, pp. 313–329, 2000.
  * [46] G. Altarelli, R. D. Ball, and S. Forte, “Small x resummation and HERA structure function data,” Nucl. Phys., vol. B599, pp. 383–423, 2001.
  * [47] G. Altarelli, R. D. Ball, and S. Forte, “Factorization and resummation of small x scaling violations with running coupling,” Nucl. Phys., vol. B621, pp. 359–387, 2002.
  * [48] G. Altarelli, R. D. Ball, and S. Forte, “An Anomalous dimension for small x evolution,” Nucl. Phys., vol. B674, pp. 459–483, 2003.
  * [49] G. Altarelli, R. D. Ball, and S. Forte, “Small x Resummation with Quarks: Deep-Inelastic Scattering,” Nucl. Phys., vol. B799, pp. 199–240, 2008\.
  * [50] A. Sabio Vera, “An ’All-poles’ approximation to collinear resummations in the Regge limit of perturbative QCD,” Nucl. Phys., vol. B722, pp. 65–80, 2005\.
  * [51] R. D. Ball, V. Bertone, M. Bonvini, S. Marzani, J. Rojo, and L. Rottoli, “Parton distributions with small-x resummation: evidence for BFKL dynamics in HERA data,” Eur. Phys. J., vol. C78, no. 4, p. 321, 2018.
  * [52] R. Abdul Khalek et al., “Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report,” 3 2021.
  * [53] P. Agostini et al., “The Large Hadron-Electron Collider at the HL-LHC,” 7 2020.
  * [54] A. Abada et al., “FCC Physics Opportunities: Future Circular Collider Conceptual Design Report Volume 1,” Eur. Phys. J. C, vol. 79, no. 6, p. 474, 2019.
  * [55] A. Abada et al., “FCC-hh: The Hadron Collider: Future Circular Collider Conceptual Design Report Volume 3,” Eur. Phys. J. ST, vol. 228, no. 4, pp. 755–1107, 2019.
  * [56] A. Bhattacharya, R. Enberg, Y. S. Jeong, C. S. Kim, M. H. Reno, I. Sarcevic, and A. Stasto, “Prompt atmospheric neutrino fluxes: perturbative QCD models and nuclear effects,” JHEP, vol. 11, p. 167, 2016.
  * [57] J. A. Carpio, K. Murase, M. H. Reno, I. Sarcevic, and A. Stasto, “Charm contribution to ultrahigh-energy neutrinos from newborn magnetars,” Phys. Rev. D, vol. 102, no. 10, p. 103001, 2020.
  * [58] J. Kwiecinski, A. D. Martin, and A. M. Stasto, “A Unified BFKL and GLAP description of $F_{2}$ data,” Phys. Rev., vol. D56, pp. 3991–4006, 1997\.
  * [59] K. Kutak and S. Sapeta, “Gluon saturation in dijet production in p-Pb collisions at Large Hadron Collider,” Phys. Rev. D, vol. 86, p. 094043, 2012.
  * [60] G. A. Chirilli and Y. V. Kovchegov, “Solution of the NLO BFKL Equation and a Strategy for Solving the All-Order BFKL Equation,” JHEP, vol. 06, p. 055, 2013.
  * [61] S. Marzani, R. D. Ball, P. Falgari, and S. Forte, “BFKL at next-to-next-to-leading order,” Nucl. Phys., vol. B783, pp. 143–175, 2007\.
  * [62] M. Deak, K. Kutak, W. Li, and A. M. Stasto, “On the different forms of the kinematical constraint in BFKL,” Eur. Phys. J., vol. C79, no. 8, p. 647, 2019.
  * [63] T. Lappi and H. Mäntysaari, “Direct numerical solution of the coordinate space Balitsky-Kovchegov equation at next to leading order,” Phys. Rev. D, vol. 91, no. 7, p. 074016, 2015.
  * [64] M. Ciafaloni, “Coherence Effects in Initial Jets at Small $q^{2}/s$,” Nucl. Phys., vol. B296, pp. 49–74, 1988.
  * [65] B. Andersson, G. Gustafson, and J. Samuelsson, “The Linked dipole chain model for DIS,” Nucl. Phys., vol. B467, pp. 443–478, 1996.
  * [66] J. Kwiecinski, A. D. Martin, and P. J. Sutton, “Constraints on gluon evolution at small x,” Z. Phys., vol. C71, pp. 585–594, 1996.
  * [67] M. Deak, L. Frankfurt, M. Strikman, and A. M. Staśto, “Taming of preasymptotic small $x$ evolution within resummation framework,” Eur. Phys. J. C, vol. 80, no. 4, p. 315, 2020.
  * [68] A. J. Askew, J. Kwiecinski, A. D. Martin, and P. J. Sutton, “QCD predictions for deep inelastic structure functions at HERA,” Phys. Rev. D, vol. 47, pp. 3775–3782, 1993.
  * [69] A. Bialas, H. Navelet, and R. B. Peschanski, “QCD dipole model and $k_{T}$ factorization,” Nucl. Phys. B, vol. 593, pp. 438–450, 2001.
  * [70] A. Bialas, H. Navelet, and R. B. Peschanski, “Virtual photon impact factors with exact gluon kinematics,” Nucl. Phys. B, vol. 603, pp. 218–230, 2001\.
  * [71] K. Golec-Biernat and A. M. Stasto, “$F_{L}$ proton structure function from the unified DGLAP/BFKL approach,” Phys. Rev. D, vol. 80, p. 014006, 2009.
  * [72] A. J. Askew, J. Kwiecinski, A. D. Martin, and P. J. Sutton, “Properties of the BFKL equation and structure function predictions for HERA,” Phys. Rev. D, vol. 49, pp. 4402–4414, 1994.
  * [73] J. Berger and A. M. Stasto, “Small x nonlinear evolution with impact parameter and the structure function data,” Phys. Rev. D, vol. 84, p. 094022, 2011.
  * [74] H. Mäntysaari and B. Schenke, “Confronting impact parameter dependent JIMWLK evolution with HERA data,” Phys. Rev. D, vol. 98, no. 3, p. 034013, 2018.
  * [75] A. Donnachie and P. V. Landshoff, “Total cross-sections,” Phys. Lett. B, vol. 296, pp. 227–232, 1992.
  * [76] M. Gluck, E. Reya, and A. Vogt, “Dynamical parton distributions of the proton and small x physics,” Z. Phys. C, vol. 67, pp. 433–448, 1995.
  * [77] F. D. Aaron et al., “Combined Measurement and QCD Analysis of the Inclusive $e^{\pm}p$ Scattering Cross Sections at HERA,” JHEP, vol. 01, p. 109, 2010.
  * [78] F. D. Aaron et al., “Measurement of $D^{*\pm}$ Meson Production and Determination of $F_{2}^{c\bar{c}}$ at low $Q^{2}$ in Deep-Inelastic Scattering at HERA,” Eur. Phys. J. C, vol. 71, p. 1769, 2011.  [Erratum: Eur.Phys.J.C 72, 2252 (2012)].
  * [79] I. Balitsky and G. A. Chirilli, “Photon impact factor and $k_{T}$-factorization for DIS in the next-to-leading order,” Phys. Rev. D, vol. 87, no. 1, p. 014013, 2013.
  * [80] M. Ciafaloni, D. Colferai, G. P. Salam, and A. M. Stasto, “A Matrix formulation for small-x singlet evolution,” JHEP, vol. 08, p. 046, 2007\.

Figure 1: Structure function $F_{2}(x,Q^{2})$ as a function of $x$ for fixed
values of $Q^{2}=2,15,35,90,150,250\,\rm GeV^{2}$, indicated next to the
curves. Solid red lines correspond to a fit with the CCSS resummed scheme.
Experimental data are from Ref. [77].  Figure 2: Charm structure function
$F_{2}^{c}(x,Q^{2})$ as a function of $x$ for fixed values of
$Q^{2}=6.5,12,20,35,60\,\rm GeV^{2}$, indicated next to the curves. Solid blue
curves indicates a fit using the CCSS resummed scheme. The experimental data
using different phase space extrapolations based on theoretical calculations
CASCADE and HVQDIS are from Ref. [78]. Figure 3: Structure function
$F_{2}(x,Q^{2})$ as a function of $x$ for fixed value of $Q^{2}=150\,\rm
GeV^{2}$ broken down into various contributions. Red dashed: $F_{2}^{(a)}$,
Eq. (25); pink dotted $F_{2}^{(b)}$, Eq. (28); blue dashed-dotted
$F_{2}^{(c)}$, Eq. (20). Finally, black solid indicates the sum of all
contributions. Figure 4: Structure function $F_{2}(x,Q^{2})$ as a function of
$x$ for fixed value of $Q^{2}=15\,\rm GeV^{2}$ broken down into various
contributions. Red dashed: $F_{2}^{(a)}$, Eq. (25); pink dotted $F_{2}^{(b)}$,
Eq. (28); blue dashed-dotted $F_{2}^{(c)}$, Eq. (20). Finally, black solid
indicates the sum of all contributions.  Figure 5: Structure function
$F_{2}(x,Q^{2})$ as a function of $x$ for fixed value of $Q^{2}=5\,\rm
GeV^{2}$ broken down into various contributions. Red dashed: $F_{2}^{(a)}$,
Eq. (25); pink dotted $F_{2}^{(b)}$, Eq. (28); blue dashed-dotted
$F_{2}^{(c)}$, Eq. (20). Finally, black solid indicates the sum of all
contributions.  Figure 6: Unintegrated gluon distribution function extracted
from the fit as a function of the transverse momentum squared $k^{2}$ for
fixed $x=0.01$. Unintegrated gluon from this work based on CCSS resummation
(red solid) is compared with the other two models, from Ref. [59], KS linear
(dashed orange) and KS non-linear (dashed-dotted blue).  Figure 7:
Unintegrated gluon distribution function extracted from the fit as a function
of the transverse momentum squared $k^{2}$ for fixed $x=0.001$. Unintegrated
gluon from this work based on CCSS resummation (red solid) is compared with
the other two models, from Ref. [59], KS linear (dashed orange) and KS non-
linear (dashed-dotted blue). Figure 8: Unintegrated gluon distribution
function extracted from the fit as a function of $x$ for fixed value of the
transverse momentum squared $k^{2}=10\,\rm GeV^{2}$. Unintegrated gluon from
this work based on CCSS resummation (red solid) is compared with the other two
models, from Ref. [59], KS linear (dashed orange) and KS non-linear (dashed-
dotted blue). Figure 9: Unintegrated gluon distribution function extracted
from the fit as a function of $x$ for fixed value of the transverse momentum
squared $k^{2}=100\,\rm GeV^{2}$. Unintegrated gluon from this work based on
CCSS resummation (red solid) is compared with the other two models, from Ref.
[59], KS linear (dashed orange) and KS non-linear (dashed-dotted blue). Figure
10: Unintegrated gluon distribution function extracted from the fit as a
function of $x$ for fixed value of the transverse momentum squared
$k^{2}=1000\,\rm GeV^{2}$. Unintegrated gluon from this work based on CCSS
resummation (red solid) is compared with the other two models, from Ref. [59],
KS linear (dashed orange) and KS non-linear (dashed-dotted blue). Figure 11:
Ratios of unintegrated gluon distributions as a function of $x$ for fixed
value of the transverse momentum $k^{2}=10\;\rm GeV^{2}$. Red solid: ratio of
the gluon extracted in this work to the KS linear; blue dashed, ratio of the
gluon extracted in this work to the KS non-linear.  Figure 12: Ratios of
unintegrated gluon distributions as a function of $x$ for fixed value of the
transverse momentum $k^{2}=100\;\rm GeV^{2}$. Red solid: ratio of the gluon
extracted in this work to the KS linear; blue dashed, ratio of the gluon
extracted in this work to the KS non-linear. Figure 13: Ratios of unintegrated
gluon distributions as a function of $x$ for fixed value of the transverse
momentum $k^{2}=1000\;\rm GeV^{2}$. Red solid: ratio of the gluon extracted in
this work to the KS linear; blue dashed, ratio of the gluon extracted in this
work to the KS non-linear. Figure 14: Effective power $\lambda_{\rm eff}$ from
Eq. (34) of unintegrated gluon distribution as a function of $x$ for fixed
value of the transverse momentum squared $k^{2}=10\;\text{GeV}^{2}$. Solid
red, this work; orange dashed, KS linear; blue dashed-dotted KS non-linear.
Figure 15: Effective power $\lambda_{\rm eff}$ from Eq. (34) of unintegrated
gluon distribution as a function of $x$ for fixed value of the transverse
momentum squared $k^{2}=100\;\text{GeV}^{2}$. Solid red, this work; orange
dashed, KS linear; blue dashed-dotted KS non-linear. Figure 16: Effective
power $\lambda_{\rm eff}$ from Eq. (34) of unintegrated gluon distribution as
a function of $x$ for fixed value of the transverse momentum squared
$k^{2}=1000\;\text{GeV}^{2}$. Solid red, this work; orange dashed, KS linear;
blue dashed-dotted KS non-linear.


# Defrosting and Blast Freezing Dark Matter

Marcos M. Flores  Department of Physics and Astronomy, University of
California, Los Angeles
Los Angeles, California, 90095-1547, USA  Chris Kouvaris  Physics Division,
National Technical University of Athens,
15780 Zografou Campus, Athens, Greece  Alexander Kusenko  Department of
Physics and Astronomy, University of California, Los Angeles
Los Angeles, California, 90095-1547, USA  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), UTIAS
The University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Theoretical Physics
Department, CERN, 1211 Geneva 23, Switzerland

###### Abstract

We show that the present-day dark matter abundance can be produced through a
novel mechanism that involves a very rapid thermal freeze-out caused by
inhomogeneous heating and successive fast cooling of small fireballs in the
early Universe. The fireballs can be produced from energy deposited in small
scale structure growth induced by Yukawa interactions in certain particle
species. Yukawa interactions are known to cause growth of halos even during a
radiation dominated era, and the same interactions facilitate cooling and
collapse of the halos by the emission of scalars. Energy deposited in the
Standard Model plasma at the locations of the halo collapse can heat the
plasma, re-establishing thermal equilibrium. The subsequent expansion and
cooling of plasma fireballs leads to freeze-out of dark matter on timescales
much shorter than the Hubble time. This mechanism can produce the right
abundance of dark matter for masses and annihilation cross sections previously
thought to be ruled out.

††preprint: IPMU23-0015

The weakly interacting massive particle (WIMP) is a well-motivated dark matter
(DM) candidate. The most commonly assumed production scenario is based on
freeze-out: the DM abundance is frozen at the temperature at which the WIMP
annihilation rate becomes slower than the expansion rate of the Universe.
Thus, it is the Hubble rate that determines the WIMP abundance.

However, there may be another relevant timescale affecting freeze-out. A
recently discovered phenomenon of halo formation in some particle species
during the radiation dominated era Amendola  _et al._ (2018); Savastano  _et
al._ (2019); Flores and Kusenko (2021); Domènech and Sasaki (2021); Flores
_et al._ (2022a); Domènech  _et al._ (2023) can create inhomogeneous heating
of plasma, with subsequent cooling of the produced fireballs, which introduces
a new timescale, much shorter than the Hubble timescale. The WIMP freeze-out
is then determined by this shorter timescale, rather than the Hubble rate,
leading to a different dependence of the DM abundance on the annihilation
cross sections. In this paper we will explore the implications of defrosting
and blast-freezing plasma for WIMP abundance. We will show that this
possibility opens a new range of WIMP parameters, which has important
implications for direct and indirect DM searches.

The traditional DM formation scenario involves a heavy particle $X$ which is
weakly coupled to the Standard Model (SM) early in the evolution of the
Universe. At high temperatures, the $X$ population is initially in thermal
equilibrium with the SM. As the Universe expands, the DM abundance is diluted
until $XX\leftrightarrow{\rm SM}$ interactions occur slowly compared to the
Hubble rate. Once interactions become rare, the comoving number density of $X$
particles remains fixed to the present day. This “freeze-out” process is
described by the Boltzmann equation,

$\dot{n}_{X}+3Hn_{X}=-\langle\sigma_{\rm
ann}v\rangle\left(n_{X}^{2}-(n_{X}^{\rm eq})^{2}\right),$ (1)

where $\langle\sigma_{\rm ann}v\rangle$ is the thermally averaged cross
section times the relative particle velocity. The temperature at which the
final DM abundance is frozen out, $T_{\rm FO}^{X}$, can be approximated by
solving

$\Gamma_{\rm ann}(T_{\rm FO}^{X})\equiv\langle\sigma_{\rm ann}v\rangle
n_{X}^{\rm eq}(T_{\rm FO}^{X})=H(T_{\rm FO}^{X}).$ (2)

The present-day $X$ abundance is given by Baumann (2022)

$\Omega_{X}\simeq 5.2\times 10^{-2}\frac{x_{\rm
FO}^{X}}{\sqrt{g_{*}(M_{X})}}\left(\frac{10^{-8}\ {\rm
GeV}^{-2}}{\langle\sigma_{\rm ann}v\rangle}\right),$ (3)

where $x_{\rm FO}^{X}=m_{X}/T_{\rm FO}^{X}$. The above result is insensitive
to the mass of the heavy particle, and primarily determined by the cross
section. The fact that the observed DM particle density is reproduced if

$\sqrt{{\langle\sigma_{\rm ann}v\rangle}}\sim 0.1\sqrt{G_{F}},$ (4)

where $G_{F}$ is the Fermi four-interaction strength, is referred to as the
WIMP miracle.

Simultaneously, Eq. (3) illustrates the fact that larger cross sections can
lead to an underabundance of DM. The principle goal of this paper is to
introduce a scenario which can produce the proper DM abundance, particularly
in the high-annihilation cross section parameter space where the traditional
freeze-out formalism produces too little DM, thus failing to generate the
abundance observed today.

This scenario relies on the recently explored principle of early structure
formation. The generation of overdensities by long-range forces, like Yukawa
interactions, has been examined in many contexts Amendola  _et al._ (2018);
Savastano  _et al._ (2019); Flores and Kusenko (2021); Domènech and Sasaki
(2021); Flores  _et al._ (2022a); Domènech  _et al._ (2023). Yukawa forces are
generally stronger than gravity, thus allowing for the formation of structure
during both the matter and radiation dominated eras. The growth of
overdensities through Yukawa forces can lead to bound states Wise and Zhang
(2014); Gresham  _et al._ (2018); Gresham and Zurek (2019), or in the presence
of radiative cooling via the same Yukawa interaction, collapse and formation
of primordial black holes Flores and Kusenko (2021, 2023). Alternatively, an
overdensity may evaporate due to annihilation of its constituent particles.
However, if these particles couple to the SM, the formation, collapse, and
annihilation of an overdense region can locally heat the SM plasma. This
inhomogeneous, local heating has proven useful when applied to either the
matter antimatter asymmetry of the Universe Asaka  _et al._ (2004); Flores
_et al._ (2022b) or the generation of primordial magnetic fields Durrer and
Kusenko (2022).

Following Ref. Flores and Kusenko (2021) we consider a dark (sub)sector with a
heavy fermion $\psi$ and a light scalar $\chi$ interacting via Yukawa
coupling:

$\mathcal{L}\supset\frac{1}{2}m_{\chi}^{2}\chi^{2}+y\chi\bar{\psi}\psi+\mathcal{L_{\rm
Y-SM}}.$ (5)

This sector, which we will refer to as the Yukawa sector throughout, is
introduced in addition to the SM and the WIMP $X$ (which may be accompanied by
some additional new physics). The interactions in the Yukawa sector of Eq. (5)
are designed to create the inhomogeneous heating. The Yukawa sector is weakly
coupled to the SM via the cross terms in $\mathcal{L_{\rm Y-SM}}$ and is not
directly linked to the WIMP sector. A schematic of the individual sectors can
be found in Fig. 1 . We will parametrize the strength of the SM-to-Yukawa
sector coupling below, when we discuss the energy transfer from the Yukawa
fireballs to the SM plasma.

Figure 1: Schematic diagram of the sectors utilized in our scenario. Here,
$\xi_{s}$ is the efficiency of energy transfer from the Yukawa sector to the
SM [see discussion below Eq. (17)], and $\alpha_{X}$ is the WIMP to SM
coupling.

We require that the fermions $\psi$ are either stable or have a total decay
width $\Gamma\ll m_{\psi}^{2}/M_{\rm Pl}$ where $M_{\rm Pl}$ is the reduced
Planck mass with numerical value, $M_{\rm Pl}\approx 2.4\times 10^{18}$ GeV.
This ensures there is a cosmological epoch where the $\psi$ particles can
become nonrelativistic, decoupled from equilibrium, and interact via the long-
range force mediated by the $\chi$ field.

The strength of the Yukawa interaction is generally much larger than gravity.
This is demonstrated by comparing the strength of each force, i.e., through
$\beta\equiv yM_{\rm Pl}/m_{\psi}\gg 1$. It should also be briefly noted that
another key difference between Yukawa interactions and gravity is the fact
that Yukawa interactions couple to the number density of $\psi$ rather than
its energy density.

Figure 2: A visualization of $\psi$ halo collapse which leads to inhomogeneous
heating.

The presence of additional long-range scalar forces generically leads to the
rapid development of structure as the overdensities
$\Delta\equiv(n_{\psi}-\bar{n}_{\psi})/\bar{n}_{\psi}$ grow considerably fast.
In particular, it has be demonstrated that Amendola  _et al._ (2018);
Savastano  _et al._ (2019); Flores and Kusenko (2021); Domènech and Sasaki
(2021); Flores  _et al._ (2022a)

$\Delta\propto a^{\beta},\quad\beta\gg 1,$ (6)

where $a$ is the scale factor, even during radiation domination. For
reference, matter perturbations under the influence of gravity grow as
$\delta\propto\ln a$ during radiation domination and $\delta\propto a$ during
matter domination. The rapid growth of structures is generally faster than the
Hubble rate, implying that the overdensities become nonlinear within a Hubble
time. This is followed by the formation of virializied $\psi$ halos comprised
of $\psi$ particles.

The details of this early structure formation have been explored both
analytically and numerically in Ref. Domènech  _et al._ (2023). In this study,
the quartic terms in the $\chi$ potential are included. Both analytical
results and $N$-body simulations point to the possibility of a rapid structure
formation in the presence of the background dynamics of the scalar field
$\chi$.

For the formation of overdensities to occur, we require that
$\bar{\psi}\psi\leftrightarrow\chi\chi$ interactions freeze out so that a
fixed population of $\psi$ particles can be captured into $\psi$ halos. To do
so, we will use a similar framework as described above and define the $\psi$
freeze-out temperature as the solution to the following equation:

$\frac{\Gamma(T_{\rm FO}^{\psi})}{H(T_{\rm FO}^{\psi})}=1\ {\rm for}\
\Gamma(T)\simeq\frac{y^{4}}{4\pi(T^{2}+m_{\psi}^{2})}n_{\psi}^{\rm eq}(T).$
(7)

Once the temperature of the $\psi$ fluid reaches $T_{\rm FO}^{\psi}$, scalar
forces can begin to coalesce material into $\psi$ halos, as previously
described.

Before calculating this temperature we must note that unlike with
gravitational forces, the binding energy of the Yukawa interactions can
contribute to the total energy budget in a nontrivial fashion. To accommodate
for this possibility, we include an additional energy component into the
equation describing the evolution of the Hubble parameter,

$3M_{\rm Pl}^{2}H^{2}=\rho_{\rm rad}+\rho_{\psi}+\rho_{y},$ (8)

where $\rho_{y}$ accounts for the energy density of Yukawa potential energy.
The length scale of the scalar force $m_{\chi}^{-1}$ requires that we consider
two regimes, namely when $H^{-1}<m_{\chi}^{-1}$ or $H^{-1}>m_{\chi}^{-1}$.
When the horizon is smaller than the Compton wavelength of the mediator, the
entire Hubble volume is subject to influence from the scalar interaction.
Alternatively, when the horizon grows beyond $m_{\chi}^{-1}$ only subhorizon
regions can communicate via the Yukawa force. In this case the number of
regions subject to scalar interactions within the horizon is
$N_{h}=(m_{\chi}/H)^{3}$. The relationship between the horizon size and the
mediator mass gives two expressions for the Yukawa energy density,

$\rho_{y}(T)=\frac{3y^{2}}{4\pi m_{\psi}^{2}H^{-3}}\begin{cases}M_{\rm
hor}^{2}/H^{-1}&H^{-1}<m_{\chi}^{-1}\\\ N_{h}M_{\rm
hal}^{2}/m_{\chi}^{-1}&H^{-1}>m_{\chi}^{-1}\end{cases},$ (9)

where

$\begin{pmatrix}M_{\rm hor}\\\ M_{\rm
hal}\end{pmatrix}=\frac{4\pi}{3}m_{\psi}n_{\psi}^{\rm
eq}(T)\begin{pmatrix}H(T)^{-3}\\\ m_{\chi}^{-3}\end{pmatrix}.$ (10)

Depending on the selection of $\\{m_{\psi},y,m_{\chi}\\}$ we have three
relevant temperatures. The Universe originally begins from a radiation
dominated era which eventually transitions to Yukawa domination at $T_{\rm
eq}^{\rm RD\to YD}$. After this, the horizon size grows to exceed
$m_{\chi}^{-1}$ at $T_{m_{\chi}=H}$. Later on as the Universe keeps expanding,
the number density of the $\psi$ fluid rapidly decreases as the temperature
falls below $m_{\psi}$. This allows for the reestablishment of radiation
domination at $T_{\rm eq}^{\rm YD\to RD}$. This leads to an evolution of the
Hubble parameter given by

$H(T)^{2}=\begin{cases}\frac{\pi^{2}}{90}g_{*}\frac{T^{4}}{M_{\rm
Pl}^{2}}&T\lesssim T_{\rm eq}^{\rm YD\to RD}\ \&\ T\gtrsim T_{\rm eq}^{\rm
RD\to YD}\\\\[7.11317pt] \frac{2\pi^{1/2}}{3M_{\rm Pl}}yn_{\psi}^{\rm
eq}(T)&T_{m_{\chi}=H}\lesssim T\lesssim T_{\rm eq}^{\rm RD\to
YD}\\\\[7.11317pt] \frac{4\pi}{9M_{\rm pl}^{2}}\frac{y^{2}n_{\psi}^{\rm
eq}(T)^{2}}{m_{\chi}^{2}}&T_{\rm eq}^{\rm YD\to RD}\lesssim T\lesssim
T_{m_{\chi}=H}.\end{cases}$ (11)

Our choice of parameters might lead to the situation where $T_{\rm eq}^{\rm
YD\to RD}>T_{m_{\chi}=H}$. In this case the evolution of the Hubble parameter
instead follows

$H(T)^{2}=\begin{cases}\frac{\pi^{2}}{90}g_{*}\frac{T^{4}}{M_{\rm
Pl}^{2}}&T\lesssim T_{\rm eq}^{\rm YD\to RD}\ \&\ T\gtrsim T_{\rm eq}^{\rm
RD\to YD}\\\\[7.11317pt] \frac{2\pi^{1/2}}{3M_{\rm Pl}}yn_{\psi}^{\rm
eq}(T)&T_{\rm eq}^{\rm YD\to RD}\lesssim T\lesssim T_{\rm eq}^{\rm RD\to
YD}.\end{cases}$ (12)

Having determined the Hubble rate in this general framework, we can now
determine when the $\bar{\psi}\psi\leftrightarrow\chi\chi$ freeze out using
Eq. (7).

Before discussing scalar radiation, we will briefly reiterate the thermal
history up to this point. The Yukawa sector is lightly coupled to the SM, for
example, through high-dimensional operators which result from a higher energy,
UV complete theory. At some temperature $T\gtrsim T_{\rm FO}^{\psi}$, the
Yukawa sector and the SM decouple. Once $T\lesssim T_{\rm FO}^{\psi}$, early
structure formation and collapse can proceed.

The statement that the SM and the Yukawa sector decouple for $T\gtrsim T_{\rm
FO}^{\psi}$ can be understood as follows. Firstly, the smallness of the SM to
Yukawa sector coupling may lead to early freeze-out of
$\bar{\psi}\psi\leftrightarrow$ SM interactions simply by weakness of this
interaction strength. However, we only require that
$\bar{\psi}\psi\leftrightarrow$ SM interactions are frozen out once $T\sim
T_{\rm FO}^{\psi}$. As is usual for freeze-out calculations, $T_{\rm
FO}^{\psi}\sim\mathcal{O}(0.1-0.01)\ m_{\psi}$, implying that the equilibrium
number density $n_{\psi}^{\rm eq}(T_{\rm FO}^{\psi})$ is within the regime of
Boltzmann suppression. This fact, and the smallness of the SM-Yukawa sector
coupling imply that these two sectors should decouple at a temperature for
$T\gtrsim T_{\rm FO}^{\psi}$. The exact details of this decoupling process
will require an exact specification of the SM-Yukawa coupling, which we will
leave to future work.

Without energy dissipation the formation of virialized $\psi$ halos through
Yukawa interactions would be the end of the story, with the newly formed halos
either remaining stable or evaporating once the constituent particles decay.
However, the same Yukawa interaction allows for the emission of scalar
radiation much in the same way as electromagnetic interactions between charged
particles allow for energy dissipation. Initially the energy is carried away
primarily through $\psi$ pair interactions, i.e., through free-free or
bremsstrahlung emission. As the halo continues collapsing, it becomes
optically thick to the scalar mediator $\chi$, and radiation becomes trapped.
This restricts radiation to be emitted only from the surface. Given that
surface cooling is the least efficient channel of energy removal, its
associated timescale determines whether radiative cooling is rapid enough to
facilitate collapse within a Hubble time. The characteristic timescale
associated with the energy loss during the surface radiation stage is

$\tau_{\rm cool}\equiv\frac{E}{|dE/dt|}\sim R_{h}\ll H^{-1},$ (13)

thus implying that radiative collapse is swift.

For our exploration of DM generation, we assume that the formation and
collapse of the $\psi$ halos occur quickly after
$\bar{\psi}\psi\leftrightarrow\chi\chi$ freeze out. Specifically, we assume
that the formation and collapse occur rapidly so the change in average
background temperature outside of the collapsing halos is negligible.
Therefore, we will set the background temperature $T_{\rm bg}$ equal to the
$\psi$ freeze-out temperature $T_{\rm FO}^{\psi}$. At formation, we will
assume that the halos initially have radius $R_{h}=m_{\chi}^{-1}$ and masses

$M_{h}=\frac{4\pi}{3}m_{\psi}n_{\psi}(T_{\rm bg})R_{h}^{3}.$ (14)

In principle, the masses and radii of the halos should be derived from an
underlying distribution such as that described in the Press-Schechter
formalism. Existing $N$-body simulations have yet to determine the precise
nature of the mass distribution of the $\psi$ halos. Our assumption that the
halos have a similar composition is motivated by both simplicity and the fact
that the strength of Yukawa interactions will facilitate a rich merger
history, as demonstrated in Ref. Domènech  _et al._ (2023), which form halos
of a maximal radius given by the Compton wavelength of the mediator $\chi$.

Without an asymmetry in the $\psi$ population, the halos will annihilate after
the initial stage of collapse. Annihilations will begin when the average
distance between particles within the halo is less than the Compton
wavelength, i.e.,

$R_{\rm
ann}\equiv\frac{1}{m_{\psi}}\left(\frac{3}{4\pi}\frac{M_{h}}{m_{\psi}}\right)^{1/3}.$
(15)

The energy released through scalar quanta during the initial collapse is given
by

$\Delta E_{\rm emis}=\frac{y^{2}M_{h}^{2}}{m_{\psi}^{2}R_{\rm
ann}}\left(1-\frac{R_{\rm ann}}{R_{h}}\right).$ (16)

Here, $\Delta E_{\rm emis}$ designates the energy released through scalar
radiation alone. Annihilation of the halo will also release energy into the
ambient plasma, $\Delta E_{\rm ann}=\epsilon_{\rm ann}M_{h}$ where
$\epsilon_{\rm ann}\leq 1$ parametrizes the efficiency of annihilation. In
particular, $\epsilon_{\rm ann}$ encapsulates deviations from perfect
annihilation, which may be caused by many sources like, for example, the
nonsphericity of the $\psi$ halo. The total energy emitted through scalar
particles is the sum, $\Delta E\equiv\Delta E_{\rm emis}+\Delta E_{\rm ann}$.
We will assume that the Yukawa sector, which contains $\psi$ and $\chi$, is
weakly coupled to the SM (see Fig. 1). The sudden release of a large amount of
energy from collapse and halo annihilation locally heats the SM plasma (see
Fig. 2). We will assume that the collapsing halos become relativistic so that
the initial temperature of the heated region is

$T_{i}^{4}=\frac{90\xi_{s}\Delta E}{4\pi^{3}g_{*}(T_{i})R_{i}^{3}},$ (17)

where $R_{i}$ is the initial radius of the heated region, $g_{*}(T_{i})$ are
the relativistic degrees of freedom at $T_{i}$ and, crucially, $\xi_{s}$ is
the efficiency of energy transfer from the Yukawa sector $\chi$ particles to
the SM plasma. Determining the quantity $\xi_{s}$ from first principles is
nontrivial. This parameter could be viewed to encapsulate a number of effects,
and not just the coupling from the Yukawa sector to the SM, and further
investigations will take a model-dependent approach in order to determine a
better understanding of the energy transfer between the SM and the Yukawa
sector.

Once heated above the background temperature, the excessive energy spreads out
via both a shockwave and diffusion. In the first case the expanding shockwave
travels through the SM plasma at the speed of sound. Using energy
conservation, we determine the characteristic timescale associated with the
explosion to be

$\tau_{\rm
exp}\equiv\frac{T}{|dT/dt|}=\frac{4R_{i}}{\sqrt{3}}\left(1+\frac{t-t_{i}}{\sqrt{3}R_{i}}\right).$
(18)

Figure 3: Predicted DM density for a variety of DM masses and couplings. The
red line indicates values within the fireball scenario where the abundance of
$X$ is exactly the observed abundance of DM today. Black contours are the
predicted DM density in the conventional freeze-out scenario, i.e., Eq. (3).
The green region is excluded by the depletion limit Eq. (24). Here $y=0.75$,
$m_{\psi}=10^{0.3}$ GeV, $m_{\chi}=10^{-16}$ GeV, and $\xi_{s}=10^{-3}$.

Most of the energy released during the collapse occurs right before
annihilation. Therefore, we approximate $R_{i}\sim R_{\rm ann}$. The region
may also cool through diffusion. The timescale associated with this process is
approximately Asaka  _et al._ (2004)

$\tau_{\rm diff}\sim\frac{R_{i}^{2}}{4D}\left(\frac{T_{i}}{T}\right)^{8/3},$
(19)

where $D$ is a diffusion constant. As in Ref. Asaka  _et al._ (2004) we take
$D\sim 1/\gamma_{g}$ where $\gamma_{g}\sim 0.3g_{s}^{2}T$ and $g_{s}$ is the
strong coupling Braaten and Pisarski (1990). The dissipation timescale is
defined as $\tau_{\rm diss}\equiv\min\\{\tau_{\rm exp},\tau_{\rm diff}\\}$.
Generally, $\tau_{\rm exp}\ll\tau_{\rm diff}$. This means that the expanding
fireball is the most rapid method of energy transport.

Before the early structure formation and collapse, we require that the WIMP
sector is decoupled from the SM. Therefore we establish the requirement that
the background temperature $T_{\rm bg}\lesssim m_{X}$. For the parameters we
consider the $X$ density produced from traditional freeze-out must also be
unable to explain the DM abundance today.

With these conditions in mind, the collapse of the $\psi$ halos reheats local
regions above $T=m_{X}$. Within these heated regions, thermal equilibrium is
reestablished allowing DM production followed by a rapid imhomogeneous
refreeze-out. The evolution of the $X$ number density in these heated regions
is described by

$T\frac{d}{dT}\left(\frac{n_{X}(T)}{T^{3}}\right)=-\frac{\Gamma_{X}(T)}{\tau_{\rm
diss}^{-1}(T)}\left[\frac{n_{X}^{2}}{T^{3}}-\frac{(n_{X}^{\rm
eq})^{2}}{T^{3}}\right],$ (20)

where

$\Gamma_{X}(T)\simeq\frac{\alpha_{X}^{2}}{(m_{X}^{2}+T^{2})}n_{X}^{\rm
eq}(T).$ (21)

Unlike the standard paradigm, here DM freeze-out reoccurs at the temperature
where $\Gamma_{X}(T_{f})=\tau_{\rm diss}^{-1}(T_{f})$ where $T_{f}$ is the
freeze-out temperature within an expanding heated region found by solving the
above equation. The fact that the fireballs expand at a much faster rate than
the Hubble parameter causes a rapid DM (re-)freeze-out leaving a significant
DM abundance even for annihilation cross sections much larger than the weak
one required in the standard freeze-out paradigm. The resultant DM energy
density is thus given by

$\rho_{X}(T_{\rm bg})\simeq f\cdot m_{X}n_{X}^{\rm eq}(T_{f}),$ (22)

where $f$ is a volume filling factor defined as

$f\equiv N_{h}H^{3}(T_{\rm bg})R_{i}^{3}\left(\frac{T_{i}}{T_{f}}\right)^{4}.$
(23)

We will restrict our parameter space such that $f<1$, $T_{i}<M_{\rm Pl}$, and
$T_{i}>m_{\psi}$, $T_{\rm bg}$. Furthermore, to ensure that the DM produced is
cold, or nonrelativistic we require that $T_{f}<m_{X}$. Lastly, we need to
impose a condition for the produced DM. Once the fireball stops expanding,
i.e., once its temperature drops to the background value, we have to guarantee
that the generated DM density does not annihilate with a rate higher than the
Hubble expansion. In such a case we will once again have a depletion of the DM
population. This condition reads as

$\frac{3\sqrt{10}}{\pi}g_{*}^{-1/2}\frac{T_{\rm bg}M_{\rm Pl}}{\tau_{\rm
diss}(T_{f})}<T_{f}^{3}.$ (24)

In Fig. 3 we show the present-day DM abundance as a function of the DM mass,
$m_{X}$, and annihilation coupling, $\alpha_{X}$. The black contours
illustrate the abundance predicted by the traditional freeze-out scenario,
where we have taken the thermal averaged cross section to be

$\langle\sigma_{\rm ann}v\rangle\sim\alpha_{X}^{2}/m_{X}^{2}.$ (25)

The red line denotes the region of parameter space within our setup where
$\Omega_{X}=\Omega_{\rm DM}$. The green region indicates parameter space which
fails the depletion limit, Eq. (24). Figure 3 illustrates that our scenario
can generate the full abundance of DM in regions of parameter space which are
ruled out in the standard freeze-out paradigm. In the presented example, we
have used $y=0.75$, $m_{\psi}=10^{0.3}$ GeV, $m_{\chi}=10^{-16}$ GeV, and
$\xi_{s}=10^{-3}$. This slice of parameter space satisfies the conditions
discussed previously, but is far from the only acceptable set of parameters.

As touched upon in Flores and Kusenko (2021); Flores  _et al._ (2022b), the
scalar $\chi$ may act as an additional relativistic degree of freedom thus
leading to a modification to $\Delta N_{\rm eff}$. This is permissible, and
may even help resolve the Hubble tension Bernal  _et al._ (2016); Gelmini  _et
al._ (2021); Anchordoqui and Perez Bergliaffa (2019); Vattis  _et al._ (2019);
Escudero and Witte (2020); Gelmini  _et al._ (2020); Vagnozzi (2020); Wong
_et al._ (2020); Perivolaropoulos and Skara (2022). Crucially, the collapse of
$\psi$ halos is generically asymmetric, thus leading to a nonvanishing
quadrupole moment. This naturally implies the generation of gravitational
waves, and studies have shown that such a signal may be detectable by future
gravitational wave observatories Flores  _et al._ (2022a).

We note that any inhomogeneities in the production of DM might occur on such
small scales that they have no implications for conventional structure
formation or isocurvature fluctuations. The mechanism described above occurs
far before matter-radiation equality. In fact, for the parameter space we have
considered, the generation of DM occurs before neutrino decoupling. The scales
upon which the DM homogeneities are produced are smaller than the horizon size
at that time, and are significantly smaller than those which could ultimately
impact the formation of large scale structures.

Here we demonstrate that early structure formation from Yukawa interactions
can lead to the generation of DM. This occurs due to the collapse of halos of
heavy fermions due to scalar radiation, which exchanges energy with SM, thus
allowing for local heating of the background SM plasma. These inhomogenous hot
spots re-establish thermal equilibrium, allowing for the local formation of
DM. Freeze-out, now dependent on the timescale pertaining to the rapid
expansion of the heated region, produces a non-negligible abundance of matter
which may present itself as the DM observed today. This is particularly
applicable to regions of parameter space which are ruled out in the
traditional freeze-out paradigm.

###### Acknowledgements.

The work of A.K. and M.M.F was supported by the U.S. Department of Energy
(DOE) Grant No. DE-SC0009937. M.M.F was also supported by the University of
California, Office of the President Dissertation Year Fellowship and donors to
the UCLA Department of Physics & Astronomy. A.K. was also supported by World
Premier International Research Center Initiative (WPI), MEXT, Japan, and by
Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No.
JP20H05853. This work used computational and storage services associated with
the Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research
and Education’s Research Technology Group.

## References

  * Amendola  _et al._ (2018) L. Amendola, J. Rubio, and C. Wetterich, Phys. Rev. D 97, 081302 (2018), arXiv:1711.09915 [astro-ph.CO] .
  * Savastano  _et al._ (2019) S. Savastano, L. Amendola, J. Rubio, and C. Wetterich, Phys. Rev. D 100, 083518 (2019), arXiv:1906.05300 [astro-ph.CO] .
  * Flores and Kusenko (2021) M. M. Flores and A. Kusenko, Phys. Rev. Lett. 126, 041101 (2021), arXiv:2008.12456 [astro-ph.CO] .
  * Domènech and Sasaki (2021) G. Domènech and M. Sasaki, JCAP 06, 030 (2021), arXiv:2104.05271 [hep-th] .
  * Flores  _et al._ (2022a) M. M. Flores, A. Kusenko, and M. Sasaki, (2022a), arXiv:2209.04970 [astro-ph.CO] .
  * Domènech  _et al._ (2023) G. Domènech, D. Inman, A. Kusenko, and M. Sasaki, (2023), arXiv:2304.13053 [astro-ph.CO] .
  * Baumann (2022) D. Baumann,  _Cosmology_ (Cambridge University Press, 2022).
  * Wise and Zhang (2014) M. B. Wise and Y. Zhang, Phys. Rev. D 90, 055030 (2014), [Erratum: Phys.Rev.D 91, 039907 (2015)], arXiv:1407.4121 [hep-ph] .
  * Gresham  _et al._ (2018) M. I. Gresham, H. K. Lou, and K. M. Zurek, Phys. Rev. D 97, 036003 (2018), arXiv:1707.02316 [hep-ph] .
  * Gresham and Zurek (2019) M. I. Gresham and K. M. Zurek, Phys. Rev. D 99, 083008 (2019), arXiv:1809.08254 [astro-ph.CO] .
  * Flores and Kusenko (2023) M. M. Flores and A. Kusenko, JCAP 05, 013 (2023), arXiv:2108.08416 [hep-ph] .
  * Asaka  _et al._ (2004) T. Asaka, D. Grigoriev, V. Kuzmin, and M. Shaposhnikov, Phys. Rev. Lett. 92, 101303 (2004), arXiv:hep-ph/0310100 .
  * Flores  _et al._ (2022b) M. M. Flores, A. Kusenko, L. Pearce, and G. White, (2022b), arXiv:2208.09789 [hep-ph] .
  * Durrer and Kusenko (2022) R. Durrer and A. Kusenko, (2022), arXiv:2209.13313 [astro-ph.CO] .
  * Braaten and Pisarski (1990) E. Braaten and R. D. Pisarski, Phys. Rev. D 42, 2156 (1990).
  * Bernal  _et al._ (2016) J. L. Bernal, L. Verde, and A. G. Riess, JCAP 10, 019 (2016), arXiv:1607.05617 [astro-ph.CO] .
  * Gelmini  _et al._ (2021) G. B. Gelmini, A. Kusenko, and V. Takhistov, JCAP 06, 002 (2021), arXiv:1906.10136 [astro-ph.CO] .
  * Anchordoqui and Perez Bergliaffa (2019) L. A. Anchordoqui and S. E. Perez Bergliaffa, Phys. Rev. D 100, 123525 (2019), arXiv:1910.05860 [astro-ph.CO] .
  * Vattis  _et al._ (2019) K. Vattis, S. M. Koushiappas, and A. Loeb, Phys. Rev. D 99, 121302 (2019), arXiv:1903.06220 [astro-ph.CO] .
  * Escudero and Witte (2020) M. Escudero and S. J. Witte, Eur. Phys. J. C 80, 294 (2020), arXiv:1909.04044 [astro-ph.CO] .
  * Gelmini  _et al._ (2020) G. B. Gelmini, M. Kawasaki, A. Kusenko, K. Murai, and V. Takhistov, JCAP 09, 051 (2020), arXiv:2005.06721 [hep-ph] .
  * Vagnozzi (2020) S. Vagnozzi, Phys. Rev. D 102, 023518 (2020), arXiv:1907.07569 [astro-ph.CO] .
  * Wong  _et al._ (2020) K. C. Wong _et al._ , Mon. Not. Roy. Astron. Soc. 498, 1420 (2020), arXiv:1907.04869 [astro-ph.CO] .
  * Perivolaropoulos and Skara (2022) L. Perivolaropoulos and F. Skara, New Astron. Rev. 95 (2022), 10.1016/j.newar.2022.101659, arXiv:2105.05208 [astro-ph.CO] .