# Efficient Training of Deep Equilibrium Models

Bac Nguyen  Lukas Mauch

###### Abstract

Deep equilibrium models (DEQs) have proven to be very powerful for learning
data representations. The idea is to replace traditional (explicit)
feedforward neural networks with an implicit fixed-point equation, which
allows to decouple the forward and backward passes. In particular, training
DEQ layers becomes very memory-efficient via the implicit function theorem.
However, backpropagation through DEQ layers still requires solving an
expensive Jacobian-based equation. In this paper, we introduce a simple but
effective strategy to avoid this computational burden. Our method relies on
the Jacobian approximation of Broyden’s method after the forward pass to
compute the gradients during the backward pass. Experiments show that simply
re-using this approximation can significantly speed up the training while not
causing any performance degradation.

Deep equilibrium models, implicit layer, efficient training

##  1 Introduction

Recently, deep equilibrium models (DEQs) have been successfully applied to a
variety of tasks, such as semantic segmentation (Bai et al., 2020a), music
source separation (Koyama et al., 2022), reverse imaging (Gilton et al.,
2021), image classification, and language modeling (Bai et al., 2019). Unlike
traditional (explicit) feedforward neural networks, DEQs define their output
as the solution to a fixed-point equation. Alternatively, DEQs can be viewed
as infinite-depth weight-tied feedforward networks. Training and inference
typically rely on fixed-point or root-finding solvers, such as Anderson
acceleration (Walker & Ni, 2011) and Broyden’s method (Bai et al., 2022), for
efficient computation during the forward and backward passes.

From a computational point of view, DEQs have some properties that make them
worth studying for efficient training and inference. 1) DEQs decouple the
network definition and implementation. While the behavior of the DEQ is
defined by the fixed-point equation we can freely choose how to solve it,
meaning that we can adapt to the needs of specific platforms. 2) Fixing the
tolerance with which we want to find the fixed point makes the computational
complexity of the forward pass data-dependent. More specifically, we can often
abort the forward iterations after just a few steps in case of simple inputs.
In practice, this yields a variable depth network (Bai et al., 2019). 3)
Compared to explicit models, training DEQs is memory-efficient. According to
the implicit function differentiation theorem, for a given input, the gradient
of the training loss with respect to the DEQ parameters only depends on the
fixed point of the DEQ but not on intermediate forward iterations. Therefore,
the memory requirement of DEQ training is constant with the network depth and
scales linearly with the layer width (Bai et al., 2019).

Despite these interesting properties, training DEQs is still far from being
efficient. The exact gradient computation in the backward pass involves an
inverse Jacobian matrix. This could be intractable in high-dimensional
features with a cubic time and quadratic memory complexity. To avoid the
explicit computation and storage of the Jacobian matrix, previous methods
proposed solving a Jacobian-based equation using a black-box root finder. In
this way, only vector-Jacobian products are involved (Bai et al., 2022), which
can be efficiently computed via automatic differentiation. Although the
Jacobian matrix is not explicitly constructed, in practice, the backward pass
is still computationally expensive since it depends on a number of iterations.

In this paper, we propose a simple approach to avoid this expensive iterative
procedure during the backward pass. In particular, we re-use the approximation
of the inverse Jacobian matrix given by Broyden’s method (Broyden, 1965) after
the forward pass to perform the gradient estimation. As a result, the backward
pass can be efficiently computed as it only requires a single matrix
multiplication. Empirical experiments show that this simple technique can
achieve very competitive results while being considerably faster than using a
black-box root finder.

##  2 Related works

Despite the success of DEQs in many tasks, reducing their running time is
still challenging. This has motivated a very active research area in designing
efficient forward and backward algorithms to speed up their training. Many
recent works have focused on accelerating the gradient computation during the
backward pass (Fung et al., 2021; Geng et al., 2021). The idea is to replace
the exact gradient computation with an inexpensive gradient approximation.
Besides the computational efficiency, approximated gradients can also
introduce noises, which have proven to be beneficial for training neural
networks (An, 1996; Zhu et al., 2019; Wu et al., 2020).

In particular, Fung et al. (2021) proposed Jacobian-Free backpropagation (JFB)
to avoid solving the Jacobian-based equation. During the backward pass, it
simply replaces the inverse Jacobian-based term with an identity matrix. Under
certain conditions, it was shown that JFB gives ascent directions of the loss
function in the training process. Similarly, Geng et al. (2021) introduced the
Neumann-series-based phantom gradients (NPG) that use a truncated Neumann
series to approximate the costly inverse Jacobian-based term. Interestingly,
JFB is a zeroth-order approximation to the Neumann series, which can be seen
as a special case of NPG. More details will be discussed in Subsection 3.2.

##  3 Proposed method

In this section, we first review DEQs and then discuss how to apply the
implicit function theorem for the gradient computation. Because implicit
differentiation requires an iterative procedure to solve a system of linear
equations, it makes the gradient computation very slow. To reduce this
computational burden, we propose GDEQ, a simple method to compute the backward
pass which only involves a simple matrix computation. As a result, it yields a
faster training scheme and a simpler implementation.

###  3.1 Deep equilibrium models

Let $f_{\bm{\theta}}(\mathbf{z},\mathbf{x})$ parameterized by
$\bm{\theta}\in\mathbb{R}^{d_{\bm{\theta}}}$ be a DEQ layer with input
$\mathbf{x}\in\mathbb{R}^{d_{\mathbf{x}}}$ and state
$\mathbf{z}\in\mathbb{R}^{d_{\mathbf{z}}}$. Furthermore, let $\mathbf{z}^{*}$
be the equilibrium point of $f_{\bm{\theta}}$ that satisfies

$\displaystyle g_{\bm{\theta}}(\mathbf{z^{*}},\mathbf{x})$
$\displaystyle:=f_{\bm{\theta}}(\mathbf{z^{*}},\mathbf{x})-\mathbf{z^{*}}=\mathbf{0}\,.$
(1)

Once the equilibrium is found, it is passed to another post-processing layer
to compute the prediction,
$\mathbf{\hat{y}}=\mathcal{F}(\mathbf{z^{*}})\in\mathbb{R}^{q}$. Let
$\mathcal{L}:\mathbb{R}^{q}\times\mathbb{R}^{q}\to\mathbb{R}$ be a
differentiable function, the network is trained based on the following
expected loss

$\displaystyle\ell=\mathbb{E}_{(\mathbf{x},y)}[\mathcal{L}(\mathbf{\hat{y}},\mathbf{y})]\,,$

where $\mathbf{y}\in\mathbb{R}^{q}$ indicates the ground-truth corresponding
to the training example $\mathbf{x}$.

Forward pass. For a given input $\mathbf{x}$, the goal is to find the
equilibrium point $\mathbf{z^{*}}$ by solving Eq. (1). Broyden’s method
(Broyden, 1965) is used to iteratively update the estimation of the
equilibrium. More specifically, let $\mathbf{z}_{0}$ be an initial random
state of the equilibrium layer. At the $t$-th iteration, we apply the
following update rule

$\displaystyle\mathbf{z}_{t+1}=\mathbf{z}_{t}-\mathbf{B}^{-1}_{t}g_{\bm{\theta}}(\mathbf{z}_{t},\mathbf{x})\,,$

where $\mathbf{B}^{-1}_{t}$ denotes an approximation of the inverse Jacobian
matrix $\nabla^{-1}_{\mathbf{z}}g_{\bm{\theta}}(\mathbf{z}_{t},\mathbf{x})$ at
$\mathbf{z}_{t}$. At the beginning, $\mathbf{B}^{-1}_{0}$ is set to be
$-\mathbf{I}$. To keep this matrix close to the Jacobian, we update it after
each iteration using the Sherman-Morrison formula (Sherman & Morrison, 1950)

$\displaystyle\mathbf{B}^{-1}_{t+1}=\mathbf{B}^{-1}_{t}-\frac{\mathbf{B}^{-1}_{t}g_{\bm{\theta}}(\mathbf{z}_{t+1},\mathbf{x})\Delta\mathbf{z}_{t+1}^{\top}\mathbf{B}^{-1}_{t}}{\Delta\mathbf{z}_{t+1}^{\top}\Big{(}\Delta\mathbf{z}_{t+1}+\mathbf{B}^{-1}_{t}g_{\bm{\theta}}(\mathbf{z}_{t+1})\Big{)}}\,,$

where $\Delta\mathbf{z}_{t+1}=\mathbf{z}_{t+1}-\mathbf{z}_{t}$. The algorithm
stops when a maximum number of iterations is reached or a solution is found.

Backward pass. To train the DEQ, we need to calculate gradients with respect
to $\bm{\theta}$ and $\mathbf{x}$. Using the implicit function theorem (Krantz
& Parks, 2002), we can compute the gradients with respect to the model
parameters $\bm{\theta}$ as

$\displaystyle\frac{\partial\ell}{\partial\bm{\theta}}=\frac{\partial
f}{\partial\bm{\theta}}\mathbf{A}\frac{\partial\ell}{\partial\mathbf{z^{*}}}\,$
(2)

where

$\displaystyle\mathbf{A}=\left(\mathbf{I}-\frac{\partial
f}{\partial\mathbf{z^{*}}}\right)^{-1}\,.$

A similar formulation can be derived to compute the gradients of $\ell$ with
respect to the input representation $\mathbf{x}$. The main computational
burden in Eq. (2) is dominated by the calculation of $\mathbf{A}$. It quickly
becomes intractable when $d_{\mathbf{z}}$ is large. This is due to the cubic
time complexity of the inverse term as well as the quadratic memory complexity
of the Jacobian matrix.

###  3.2 Efficient backward pass

The goal is to compute Eq. (2) in an efficient manner. Several simplifications
have been proposed. More specifically, JFB (Fung et al., 2021) simply ignores
the Jacobian inverse and uses $\mathbf{A}=\mathbf{I}$. This dramatically
reduces the computational complexity. NPG (Geng et al., 2021) uses a truncated
Neumann series to approximate the inverse, yielding

$\displaystyle\bf{A}$ $\displaystyle=\lambda\sum_{i=0}^{k-1}\mathbf{M}^{i}\,,$
$\displaystyle\bf{M}$ $\displaystyle=\lambda\frac{\partial
f}{\partial\mathbf{z^{*}}}+(1-\lambda)\bf{I}\,,$

where $k\in\mathbb{N}$ is the sequence length and $\lambda\in[0,1]$ is a
dampening factor. Instead of calculating $\left(\mathbf{I}-\frac{\partial
f}{\partial\mathbf{z^{*}}}\right)^{-1}$, we only have to compute $k-1$ vetor-
Jacobian products. Typically, we choose $k=5$ and $\lambda=0.5$ for optimal
performance (Geng et al., 2021).

To further reduce the computation, we propose to compute Eq. (2) by re-using
the approximation of the inverse Jacobian term
$\mathbf{B}_{T}^{-1}\approx-\left(\mathbf{I}-\frac{\partial
f}{\partial\mathbf{z^{*}}}\right)^{-1}$ given by Broyden’s method at the last
iteration after the forward pass. This can avoid the intractable computation
of the exact $\mathbf{A}$ in Eq. (2). More specifically, the gradient
computation is formulated as

$\displaystyle\frac{\partial\ell}{\partial\bm{\theta}}\approx-\frac{\partial
f}{\partial\bm{\theta}}\mathbf{B}_{T}^{-1}\frac{\partial\ell}{\partial\mathbf{z^{*}}}\,.$

The main advantage is that the gradient computation in the backward pass only
involves a simple matrix multiplication. This is much simpler than solving a
system of linear equations as proposed by Bai et al. (2019). Besides the
computational simplicity, it can also avoid additional errors that would be
introduced by running an additional root-finding algorithm, as proposed in the
original DEQ paper. Note that the matrix $\mathbf{B}^{-1}_{T}$ is written as a
sum of $\mathbf{B}^{-1}_{0}$ and low-rank updates
$\mathbf{B}^{-1}_{T}=\mathbf{B}^{-1}_{0}+\sum_{i=1}^{T}\mathbf{u}_{i}\mathbf{v}_{i}^{\top}$.
When the maximum memory storage $m$ is reached, we keep the last $m$ low-rank
updates and discard the oldest updates (Bai et al., 2020b). We refer to the
proposed method as GDEQ.

In addition, Table 1 reports the time and memory complexities of different
methods required for the backward pass. Here $K$ and $k$ denote the number of
iterations used in Broyden’s method and the number of unrolling steps used in
NPG, respectively. As shown in the table, JFB and GDEQ are the most efficient
since they both enjoy constant training time and memory complexities.

Table 1: Comparison of time and memory complexities among the competing methods. Method | Time | Memory
---|---|---
Implicit | $\mathcal{O}(K)$ | $\mathcal{O}(1)$
JFB | $\mathcal{O}(1)$ | $\mathcal{O}(1)$
NPG | $\mathcal{O}(k)$ | $\mathcal{O}(1)$
GDEQ | $\mathcal{O}(1)$ | $\mathcal{O}(1)$

##  4 Experiments

In this section, we conduct experiments on the CIFAR-10 dataset for image
classification to show the effectiveness of our method. Multiscale deep
equilibrium models (MDEQ) (Bai et al., 2020b) is used as the backbone
model111Code available at https://github.com/locuslab/deq. All experiments run
on 4 NVIDIA RTX A6000 GPUs and use the same configuration.

In Table 2, we report the comparison between GDEQ and other state-of-the-art
methods, namely with 1) implicit differentiation that solves the system of
linear equations during the backward pass (Implicit) (Bai et al., 2019), 2)
Jacobian-free backpropagation (JFB) (Fung et al., 2021), and 3) the Neumann-
series-based phantom gradient method (NPG) (Geng et al., 2021). Following Bai
et al. (2020a), we use Broyden’s method with 18 iterations in the forward and
20 iterations in the backward pass.

As seen from Table 2, GDEQ can significantly speed up the training time. It is
2.16x faster than the implicit differentiation. In an early stage of training,
it is common to explicitly unroll the DEQ layers as a weight-tied feedforward
neural network. The pre-trained weights are then used to initialize the
parameters of MDEQ. In our experiment, we report the classification accuracies
under both settings, with and without the pre-training step. The weight
initialization during the pre-training stage helps to stabilize the training
of MDEQ. It consistently improves the classification accuracies for most of
the competing methods. We plan to investigate why this happens in the future.

Table 2: Performances of the competing methods on the CIFAR-10 dataset. We report the training speed and the classification accuracy. Method | Speedup | Accuracy
---|---|---
without | with
Implicit | 1.00x | 91.64 | 92.29
JFB | 1.97x | 91.56 | 91.99
NPG | 1.71x | 93.42 | 92.46
GDEQ | 2.16x | 92.01 | 93.08

Figure 1 compares the convergence speed of the MDEQ with different backward
algorithms. As before, we consider implicit differentiation, JFB, NPG, and our
GDEQ. We show the accuracy over training time. Obviously, GDEQ converges
fastest, performing $220$ training epochs in just $6.7$ hours. This is about
2.16x faster than the vanilla DEQ training that uses implicit differentiation,
which takes $14.6$ hours for $220$ epochs. Note, that this is even faster than
with JFB, although the backward pass is slightly more expensive. We explain
this with the slower convergence caused by a less accurate gradient estimation
and the requirement of more forward iterations associated with it. However,
given more time GDEQ is outperformed by NPG.

Figure 1: Convergence of the MDEQ on the test set for different backward
algorithms, using no pre-training.

Furthermore, we perform an ablation study to see how similar the gradients
approximated by GDEQ are to those computed by Implicit. Figure 2 illustrates
the cosine similarity of gradients computed by the two methods. It shows that
the approximation of the inverse Jacobian given by Broyden’s method is quite
close to that of Implicit. This implies that GDEQ can still provide the ascent
direction of the loss function.

(a) Epoch 1 (b) Epoch 10 (c) Epoch 20 (d) Epoch 50

Figure 2: Cosine similarity between the gradients computed by our method and
implicit differentiation at different training epoch.

##  5 Conclusions

In this paper, we have demonstrated that DEQs can be efficiently trained by
using an approximation of the inverse Jacobian matrix given by Broyden’s
method during the forward pass. Several studies on image classification tasks
have been conducted to verify our claim. In particular, the proposed method
GDEQ achieves 2.16x faster in terms of training time, while being competitive
in terms of classification performances. Clearly, it shows great advantages
for training DEQs on large-scale datasets.

Limitations. Despite the promising results, our work has some limitations.
First, it lacks studies when Broyden’s method does not converge during the
forward pass under a limited number of iterations. Since the proposed method
heavily depends on the approximation given by Broyden’s method, a very biased
approximation may lead to unexpected results. Second, it lacks a theoretical
convergence guarantee of stochastic gradient descent using gradients provided
by our method. Finally, it is also interesting to extend our studies of GDEQ
to other task domains rather than image classification. Our future work will
aim to address these issues.

## References

  * An (1996) An, G.  The effects of adding noise during backpropagation training on a generalization performance.  _Neural computation_ , 8:643–674, 1996.
  * Bai et al. (2019) Bai, S., Kolter, J. Z., and Koltun, V.  Deep equilibrium models.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Bai et al. (2020a) Bai, S., Koltun, V., and Kolter, J. Z.  Multiscale deep equilibrium models.  In _Advances in Neural Information Processing Systems_ , volume 33, pp. 5238–5250, 2020a.
  * Bai et al. (2020b) Bai, S., Koltun, V., and Kolter, J. Z.  Multiscale deep equilibrium models.  In _Advances in Neural Information Processing Systems_ , pp. 5238–5250, 2020b.
  * Bai et al. (2022) Bai, S., Koltun, V., and Kolter, J. Z.  Neural deep equilibrium solvers.  In _Proceedings of the International Conference on Learning Representations_ , 2022.
  * Broyden (1965) Broyden, C. G.  A class of methods for solving nonlinear simultaneous equations.  _Mathematics of computation_ , 19(92):577–593, 1965.
  * Fung et al. (2021) Fung, S. W., Heaton, H., Li, Q., McKenzie, D., Osher, S., and Yin, W.  Jfb: Jacobian-free backpropagation for implicit networks.  _arXiv preprint arXiv:2103.12803_ , 2021.
  * Geng et al. (2021) Geng, Z., Zhang, X.-Y., Bai, S., Wang, Y., and Lin, Z.  On training implicit models.  In _Advances in Neural Information Processing Systems_ , pp. 24247–24260, 2021.
  * Gilton et al. (2021) Gilton, D., Ongie, G., and Willett, R.  Deep equilibrium architectures for inverse problems in imaging.  _IEEE Transactions on Computational Imaging_ , 7:1123–1133, 2021.
  * Koyama et al. (2022) Koyama, Y., Murata, N., Uhlich, S., Fabbro, G., Takahashi, S., and Mitsufuji, Y.  Music source separation with deep equilibrium models.  In _Proceedings of the International Conference on Acoustics, Speech and Signal Processing_ , pp. 296–300, 2022.
  * Krantz & Parks (2002) Krantz, S. G. and Parks, H. R.  _The Implicit Function Theorem: History, Theory, and Applications_.  Springer Science & Business Media, 2002.
  * Sherman & Morrison (1950) Sherman, J. and Morrison, W. J.  Adjustment of an inverse matrix corresponding to a change in one element of a given matrix.  _The Annals of Mathematical Statistics_ , 21:124–127, 1950\.
  * Walker & Ni (2011) Walker, H. F. and Ni, P.  Anderson acceleration for fixed-point iterations.  _SIAM Journal on Numerical Analysis_ , 49:1715–1735, 2011\.
  * Wu et al. (2020) Wu, J., Hu, W., Xiong, H., Huan, J., Braverman, V., and Zhu, Z.  On the noisy gradient descent that generalizes as sgd.  In _Proceedings of the International Conference on Machine Learning_ , pp. 10367–10376, 2020.
  * Zhu et al. (2019) Zhu, Z., Wu, J., Yu, B., Wu, L., and Ma, J.  The anisotropic noise in stochastic gradient descent: Its behavior of escaping from sharp minima and regularization effects.  In _Proceedings of the International Conference on Machine Learning_ , pp. 7654–7663, 2019.


# Constraints on Dark Matter from Dynamical Heating of Stars in Ultrafaint
Dwarfs. Part 2: Substructure and the Primordial Power Spectrum

Peter W. Graham<EMAIL_ADDRESS>Stanford Institute for Theoretical
Physics, Department of Physics, Stanford University, Stanford, CA 94305, USA
Kavli Institute for Particle Astrophysics & Cosmology, Department of Physics,
Stanford University, Stanford, CA 94305, USA  Harikrishnan Ramani
<EMAIL_ADDRESS>Stanford Institute for Theoretical Physics, Department of
Physics, Stanford University, Stanford, CA 94305, USA  Department of Physics
and Astronomy, University of Delaware, Newark, DE 19716, USA

###### Abstract

There is a large and growing interest in observations of small-scale structure
in dark matter. We propose a new way to probe dark matter structures in the
$\sim 10-10^{8}\,M_{\odot}$ range. This allows us to constrain the primordial
power spectrum over shorter distances scales than possible with direct
observations from the CMB. For $k$ in the range $\sim 10-1000\,{\rm Mpc}^{-1}$
our constraints on the power spectrum are orders of magnitude stronger than
previous bounds. We also set some of the strongest constraints on dark matter
isocurvature perturbations. Our method relies on the heating effect such dark
matter substructures would have on the distribution of stars in an ultra-faint
dwarf galaxy. Many models of inflation produce enhanced power at these short
distance scales and can thus be constrained by our observation. Further, many
dark matter models such as axion dark matter, self-interacting dark matter and
dissipative dark matter, produce dense structures which could be constrained
this way.

##  I Introduction

The nature of dark matter (DM) is one of the biggest outstanding mysteries in
physics. We have a wide variety of probes of dark matter including direct and
indirect detection and collider experiments. The study of dark matter
structure in the universe is another important probe. Dark matter over-
densities collapse in the early universe leading to non-trivial structure over
several mass and length scales including clusters, galaxies and dwarf
galaxies. The matter power spectrum obtained from correlations in the observed
structure agrees very well with the $\Lambda$CDM paradigm down to the
$\approx\textrm{Mpc}$ scale. However there has been a large and growing
interest recently in studying shorter scales, as a great deal can be learned
both about the nature of dark matter and about the formation of the universe.
While a scale invariant primordial spectrum of density perturbations should
lead to formation of structures at all length scales, there are many possible
cosmological scenarios which produce even more structure at these shorter
length scales. For example many models of inflation produce enhancements in
the primordial power spectrum at short scales, causing additional, denser
substructures to form (see e.g. [1, 2, 3, 4, 5, 6]). Additionally a period of
early matter dominance can enhance small-scale power (see e.g. [7, 8, 9, 10]).
Further, even with a scale-invariant initial power spectrum, nontrivial
dynamics in the dark sector can cause substructures to form. For example,
axion dark matter, dark photon dark matter, self-interacting dark matter
(SIDM), atomic dark matter, dissipative dark matter, and mirror sectors can
form dense substructures (see e.g. [11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26]). Thus, studying dark matter structures on these short
scales can teach us much about the origin of the universe and the nature of
dark matter.

The study of dark matter structures becomes increasingly difficult at smaller
length/mass scales below dwarf galaxy size, owing to the inability of these
dark matter halos to produce significant numbers of stars. Techniques to study
these structures have instead involved gravitational lensing [27, 28, 29] of
light, photometric lensing [30] and stellar streams [31, 32]. These promising
approaches will hopefully reveal much about the small-scale structure of dark
matter in the future. In this paper we describe a new, complementary approach
to studying dark matter structures in the mass range $\sim
10-10^{8}\,M_{\odot}$, or primordial density perturbations in the range of $k$
from $\sim 10-1000\,{\rm Mpc}^{-1}$.

These smaller structures, which we refer to as clumps in this work, collapse
before galaxies and subsequently fall into the galaxies as they form. Thus,
these clumps make up intra-galactic substructure and would survive until today
barring destruction due to tidal forces present in the host. In this work, we
describe an important probe of such sub-structure down to the $\approx
10M_{\odot}$ mass scale, specifically due to their presence inside some of the
faintest substructure identified yet: ultra-faint dwarf galaxies (UFDs).

Part one of this two-part paper set stringent limits on MACHOs/PBHs making up
some or all of the UFD dark matter [33]. These limits were obtained via
considering the gravitational scattering of MACHOs off stars (for similar
works see e.g. [34, 35, 36, 37, 38, 39, 40, 41]). Such scattering will, on
average, transfer heat from the MACHO’s kinetic energy to the stellar kinetic
energy which would result in the expansion of the stellar scale radius. This
can be compared to the observed stellar scale radius to set constraints. That
paper was essentially in the ‘point-mass’ limit, taking the MACHOs/PBHs to be
arbitrarily small. In this sequel we extend the analysis to diffuse dark
matter clumps of arbitrary density. Importantly, we still expect a similar
heat transfer rate from the dark matter clumps to the stars even if the dark
matter clumps are significantly larger than the MACHOs, and in fact even if
the dark matter clump size approaches the size of the stellar scale radius
($\sim 30\,{\rm pc}$). This is because the minimum impact parameter for
gravitational clump-star scattering enters only via a log. Thus even though
moving from point-like MACHOs to diffuse clumps will cause an increase by many
orders of magnitude in the minimal impact parameter, this only enters the log
in the scattering and so has very little effect on the heat transfer rate.
Thus this effect should be a powerful probe of even relatively diffuse DM
substructures. However for low enough densities, such clumps will not survive
for long in the environment of the UFD, thus cutting off the heating effect at
some minimal density for the clumps. In this work we analyze the survival
prospects of these diffuse objects in the UFD environment. For the clumps that
survive, we estimate the effects of the finite size of the clumps on the heat
transfer rate thus deriving limits on the existence of clumps in the UFD.

The rest of this paper is organized as follows. Section II describes our
calculation of the gravitational heat transfer rate from clumps to stars.
Section III analyzes the survival probability of clumps against various
disruption processes. We present results in Section IV and discuss
consequences in Section V.

##  II Heat Transfer for Clumps

In this section we describe our assumptions about the clumps that are inside
the UFD and compute the heat transfer rate due to gravitational scattering
with stars as well as the smooth DM.

###  II.1 Parametrization of clumps

First, we consider dark matter to reside in subhalos all of the same mass
$M_{\rm clump}$. While this assumptions simplifies the limit setting
procedure, the limit calculation is straightforward for extended halo mass
functions which we leave for future work.

We consider only NFW clumps and characterize them as

$\displaystyle\rho_{\rm
NFW}(r)=4\rho_{s}\left(\frac{r}{r_{s}}\right)^{-1}\left(\frac{r}{r_{s}}+1\right)^{-2}$
(1)

Here, $r_{s}$ is the scale radius, $\rho_{s}$ is the density at the scale
radius. Owing to the fact that the NFW cores are more robust to tidal
stripping, we follow [30] and truncate the NFW profile at the scale radius
$r_{s}$ with masses

$M_{\rm clump}=8\pi r_{s}^{3}\rho_{s}(-1+\log 4)$ (2)

These cores are assumed to make up a fraction $f_{\rm clump}$ of the total
dark matter mass of the UFD. The rest of the dark matter is assumed to be
smooth for simplicity. Thus, the three parameters $\\{M_{\rm
clump},\rho_{s},f_{\rm clump}\\}$ fully describe the clumps we consider. For
later use we define the mass enclosed inside a radius $b$ as

$\displaystyle M_{\rm enc}(b)=\int_{0}^{b}4\pi r^{2}\rho_{\rm NFW}(r)$ (3)

since the clumps are truncated at $r_{s}$, we assume that $M_{\rm
enc}(b>r_{s})=M_{\rm clump}$.

Finally, we assume the clumps themselves are distributed following a Dehnen
profile in the UFD. So the averaged mass density of all the clumps in the UFD
is:

$\displaystyle\rho_{\rm clumps}(r,t)=f_{\rm clump}\frac{3M_{\rm UFD}}{4\pi
R_{\rm clumps}^{3}(t)}\left(1+\frac{r}{R_{\rm clumps}(t)}\right)^{-4}$ (4)

with the scale radius $R_{\rm clumps}$ being time dependent due to heat
transfer between clumps and the smooth DM component causing the clumps to
migrate inwards as derived in the next subsection. The smooth DM itself also
follows a Dehnnen profile but with a constant scale radius $R_{\rm
UFD}=446.3~{}\textrm{pc}$.

$\displaystyle\rho_{\rm DM,smooth}(r,t)=\left(1-f_{\rm
clump}\right)\frac{3M_{\rm UFD}}{4\pi R_{\rm UFD}^{3}}\left(1+\frac{r}{R_{\rm
UFD}}\right)^{-4}$ (5)

We assume the smooth DM component is the dominant component of DM so
technically our results would not apply if $f_{\rm clump}$ is too close to 1.
We calculate the time dependent $R_{\rm clumps}(t)$ next.

###  II.2 Migration of clumps

Although our main effect is the heating of stars, in this subsection we
consider the heat transfer from clumps to the dominant, diffuse component of
the DM. This is the dominant heat loss mechanism for the clumps. It causes the
clumps to migrate inwards towards the center of the UFD, decreasing the scale
radius of the clump distribution over time. This increases the density of the
clumps at the location of the stars, and is thus an important effect to
consider for the heating of the stars.

In Paper 1 [33] we considered the scenario that dark matter had two
components: MACHOs (with mass $M_{\rm MACHO}$) and diffuse DM (with mass
$m_{\rm DM}$). The heat transfer rate from MACHOs to diffuse DM per unit MACHO
mass, in the $m_{\rm MACHO}\gg m_{\rm DM}$ limit was given as

$\displaystyle H_{\rm MACHO}=2\sqrt{2\pi}G^{2}\left(1-f_{\rm
MACHO}\right)\rho_{\rm DM}M_{\rm MACHO}$ $\displaystyle\times\frac{\sigma_{\rm
MACHO}^{2}}{\left(\sigma_{\rm DM}^{2}+\sigma_{\rm
MACHO}^{2}\right)^{\frac{3}{2}}}\log\left(\frac{b_{90}^{2}+b_{\rm
max}^{2}}{b_{90}^{2}+b_{\rm min}^{2}}\right)$ (6)

Here we define

$\displaystyle\rho_{\rm DM}\equiv\rho_{\rm DM,smooth}\left(r\ll R_{\rm
UFD}\right)\approx\frac{3M_{\rm UFD}}{4\pi R_{\rm UFD}^{3}}$ (7)

Eqn. 6 assumes the MACHOs are effectively point masses (as would be the case
for primordial black holes for example).

In this paper we wish to consider what happens when we replace the MACHOs with
finite-size dark matter clumps (bound together by gravity). So we now make
modifications to the above expressions in order to capture the finite size
effects of the clump. First $b_{\rm max}$ is set to the current scale radius
of the distribution of clumps in the UFD, $R_{\rm clumps}(t)$, because this is
roughly the orbit size of a typical clump. For point particles, $b_{\rm min}$
was set to zero because there are no sampling limitations for smooth DM. For
the present case case, we can conservatively restrict the impact parameter to
be larger than the scale radius of each clump, such that any DM particle sees
the whole mass of the clump $M_{\rm clump}$ i.e. we set $b_{\rm min}=r_{s}$.
Hence, we have,

$\displaystyle H_{\rm clump}=2\sqrt{2\pi}G^{2}\left(1-f_{\rm
clump}\right)\rho_{\rm DM}M_{\rm clump}$ $\displaystyle\times\frac{\sigma_{\rm
clump}^{2}}{\left(\sigma_{\rm DM}^{2}+\sigma_{\rm
clump}^{2}\right)^{\frac{3}{2}}}\log\left(\frac{b_{90}^{2}+R_{\rm
clumps}(t)^{2}}{b_{90}^{2}+r_{s}^{2}}\right)$ (8)

We then integrate this heat transfer rate over the distribution of clumps to
find the total rate of energy loss of the clumps:

$\displaystyle\frac{d\mathcal{E}_{\rm clump}(R_{\rm clumps})}{dt}=\int
dr\,4\pi r^{2}\rho_{\rm clumps}(r)\,$ $\displaystyle\times H_{\rm
clump}\left(\sigma_{\rm clump}(r),\sigma_{\rm DM}(r),\rho_{\rm DM}(r)\right)$
(9)

And then we plug in the expression for the energy of the clumps
$\mathcal{E}_{\rm clump}$ which is identical to $\mathcal{E}_{\rm MACHO}$ in
Ref. [33] with $M_{\rm MACHO}$ set by $M_{\rm clump}$. This is because the
potential energy and kinetic energy of an individual clump due to its presence
in the UFD should be the same as the corresponding MACHO since the clump
radius $r_{s}$ is much smaller than the size of the UFD, $R_{\rm UFD}$. The
initial condition is given by $R_{\rm clumps}(t=0)=R_{\rm UFD}$. And we can
then solve for the evolution of the scale radius of the clump distribution in
time, $R_{\rm clumps}(t)$.

What we really care about is the enhancement in the number density of clumps
seen by the stars (i.e. at the stellar scale radius $R_{0,\star}$) averaged
over the age of the UFD, $T_{\rm UFD}$ compared to the average DM density
without the effect of migration. This is given by

$\displaystyle\eta_{\rm NT}=\frac{1}{T_{\rm UFD}}\int_{0}^{T_{\rm
UFD}}dt\frac{\rho_{\rm clumps}(R_{0,\star}(t),t)}{\rho_{\rm clumps}(R,t=0)}$
(10)

Here, the subscript NT emphasizes that no tidal effects have been included.
Note that $R_{\rm clumps}(t)$ starts at $R_{\rm UFD}$ at $t=0$ and stops
contracting when it reaches [33] $R_{\rm clumps}^{\rm cutoff}\approx f_{\rm
clump}^{\frac{1}{3}}R_{\rm UFD}$ i.e. when the density of clumps matches the
ambient smooth DM density. Note that almost the entire range of
$10^{-4}<f_{\rm clump}<1$ that we consider in this paper corresponds to
$R_{\rm clumps}^{\rm cutoff}>24.73~{}\textrm{pc}$. Furthermore, since
$R_{0,\star}<24.73~{}\textrm{pc}$ during its evolution, it is true that
$R_{0,\star}<R_{\rm clumps}$, i.e. the stars are always inside the region
containing the maximal density of clumps. Since the core density inside
$R_{\rm clumps}$ is roughly a constant owing to the assumption of a Dehnen
profile with $\gamma=0$, we can make the approximation

$\displaystyle\rho_{\rm clumps}(R_{0,\star}(t),t)\approx f_{\rm
clump}\frac{3M_{\rm UFD}}{4\pi R_{\rm clumps}^{3}(t)}$ (11)

In other words, we define the enhancement $\eta_{\rm NT}$ at the stellar scale
radius $R_{0,\star}$ but it would be the same when evaluated at any radius
inside the scale radius of the clumps.

Figure 1: The enhancement factor $\eta(R=24.7~{}\textrm{pc})$. This is the
enhancement of the clump density at the location of the stars (averaged over
the age of the UFD). The enhancement arises from the heat transfer from clumps
to diffuse DM which causes the clumps to migrate inwards towards the center of
the UFD over time. The enhancement ignoring clump destruction $\eta_{\rm NT}$
as given in Eqn. (10) is plotted in dashed lines. The effective enhancement
$\eta_{\rm HHH}$ as given in Eqn. (30) is plotted in solid lines.

The result of solving the differential equation for $R_{\rm clumps}(t)$ and
evaluating $\eta_{\rm NT}$ is plotted in Fig. 1 as a function of the clump
mass $M_{\rm clump}$ for different choices of the density of the clump
$\rho_{s}$ for $f_{\rm clump}=10^{-3}$ in dashed lines. The solid lines are
the result after clump survival in tidal heating is included and will be
discussed in Sec. III.1. Note that the $\rho_{s}$ dependence comes through
$r_{s}$ in Eqn. 8. Since $r_{s}$ only appears in the $\log$ in Eqn. 8, the
enhancement is only weakly dependent on the clump density $\rho_{s}$. A five
order of magnitude drop in the clump density from
$10^{6}\,M_{\odot}\,\textrm{pc}^{-3}$ to $10\,M_{\odot}\,\textrm{pc}^{-3}$ can
be roughly compensated by a factor of 2 increase in the clump mass to achieve
the same enhancement.

Thus we now know what to take for the local density of clumps when we evaluate
the heating of the stars by the clumps.

###  II.3 Direct heating rate of stars by clumps

In this subsection we calculate the direct heat transfer rate from clumps to
stars. By direct heat transfer we mean clumps that pass directly through the
stellar scale radius (as opposed to tidal heating which we consider in the
next subsection).

In Paper 1 [33] the heat transfer between stars (typical mass $m_{\star}$) and
MACHOs (taken to be point objects) per unit stellar mass, in the $m_{\rm
MACHO}\gg m_{\star}$ limit, was given as

$\displaystyle H_{\star}=2\sqrt{2\pi}\rho_{\rm MACHO}(t)\frac{G^{2}m_{\rm
MACHO}\sigma_{\rm MACHO}^{2}}{\left(\sigma_{\star}^{2}+\sigma_{\rm
MACHO}^{2}\right)^{\frac{3}{2}}}$
$\displaystyle\times\log\left(\frac{b_{90}^{2}+b_{\rm
max}^{2}}{b_{90}^{2}+b_{\rm min}^{2}}\right)$ (12)

where $\eta$ is the enhancement in local number density due to migration,
$b_{\rm max}=3R_{0,\star}$ and $b_{\rm min}=b_{\rm samp}$ and here
$R_{0,\star}$ is the stellar scale radius and $b_{\rm samp}$ is the sampling
radius, i.e. the smallest radius around a star through which at least 3 MACHOs
will pass in the age of the UFD. We now modify this formula to change from
MACHOs (point particles) to extended clumps with scale radius $r_{s}$.

In the limit $b_{\rm max}>b_{\rm samp}>r_{s}$, there is no need for any
modification since for scattering purposes, the clump is a point particle due
to Newton’s shell theorem. We just make the obvious replacements $m_{\rm
MACHO}\rightarrow M_{\rm clump}$, $f_{\rm MACHO}\rightarrow f_{\rm clump}$ and
$\sigma_{\rm MACHO}\rightarrow\sigma_{\rm clump}$.

When $b_{\rm max}>r_{s}>b_{\rm samp}$, we can conservatively take $b_{\rm
min}=r_{s}$, requiring the star to pass outside the dark matter clump.
Although there will be some extra heating effect from scatterings where the
star passes through the dark matter clump, we ignore this extra heating effect
which gives us a conservative bound.

Finally, when $r_{s}>b_{\rm max}$, we can conservatively choose to use the
heating effect of only a smaller part of each clump (the part within a certain
radius $b_{\rm opt}$ of its center). We choose this radius to be small enough
that this part of the clump does fit within the stellar distribution so
$b_{\rm opt}<b_{\rm max}$. This is conservative because we are only using the
heating effect of a fraction of each clump, the rest of the clump will add
some extra heating but we neglect that to find a conservative bound. Thus for
impact parameters between $b_{\rm opt}$ and $b_{\rm max}$ we can still use our
old point mass formula Eqn. (12). We thus take the minimum impact parameter to
be $b_{\rm min}=b_{\rm opt}$. Then we also have to replace $M_{\rm
clump}\rightarrow M_{\rm enc}(b_{\rm opt})$ which is the mass enclosed inside
$b_{\rm opt}$, i.e. the mass of the smaller part of the clump we are now
considering. Since we are only considering a concentric subsphere of the
original clump and this part has a smaller mass, we would also need to modify
the total mass density of clumps in Eqn. (12) by replacing $\rho_{\rm
clumps}\rightarrow\rho_{\rm clumps}\frac{M_{\rm enc}(b)}{M_{\rm clump}}$. In
principle we could choose any $b_{\rm opt}$ so long as it was less than
$b_{\rm max}$ and it would give us some conservative lower bound on the
heating rate, however we choose the optimum $b_{\rm opt}$ which maximizes the
heating rate of the stars.

To summarize, the direct heat transfer rate from clumps to stars per unit
stellar mass is

$\displaystyle H_{\star,{\rm direct}}=2\sqrt{2\pi}\frac{G^{2}\sigma_{\rm
clump}^{2}}{\left(\sigma_{\star}^{2}+\sigma_{\rm
clump}^{2}\right)^{\frac{3}{2}}}\rho_{\rm clumps}(t)\frac{M_{\rm enc}(b_{\rm
min})^{2}}{M_{\rm clump}}$
$\displaystyle\times\log\left(\frac{b_{90}^{2}+b_{\rm
max}^{2}}{b_{90}^{2}+b_{\rm min}^{2}}\right)$ (13)

Here, the radial dependence of $\rho_{\rm clumps}$ is suppressed owing to the
choice of a Dehnen profile with $\gamma=0$ such that the DM density is near
constant inside the DM scale radius. Furthermore, as discussed in Ref. [33],
the only time dependence arises from $\rho_{\rm clumps}(t)$ and to a good
approximation, it can be replaced by its time average i.e. $\rho_{\rm
clumps}(t)\rightarrow\rho_{\rm clumps}(t=0)\times\eta_{\rm NT}$ which we do to
simplify the computation.

The quantity $b_{\rm min}$ is given by,

$\displaystyle b_{\rm min}=\begin{cases}b_{\rm samp}&b_{\rm samp}>r_{s}\\\
r_{s}&b_{\rm max}>r_{s}>b_{\rm samp}\\\ b_{\rm opt}&r_{s}>b_{\rm
max}\end{cases}$ (14)

To get a feeling for where the optimization of $b_{\rm opt}$ pushes us in the
case of $r_{s}>b_{\rm max}$, let us estimate $b_{\rm opt}$. The heating rate
is maximized when $M_{\rm enc}^{2}(b)\times\log\left(\frac{b_{90}^{2}+b_{\rm
max}^{2}}{b_{90}^{2}+b^{2}}\right)$ is maximized. For the density profile in
Eqn. 2, we find $b_{\rm opt}\approx\frac{b_{\rm max}}{\sqrt{e}}\approx
0.6b_{\rm max}$ when $b_{\rm max}\gg b_{90}$ and $b_{\rm opt}\approx 0.7b_{\rm
max}$ when $b_{\rm max}\ll b_{90}$. So in other words, the optimum pushes us
to take new, smaller clumps which are roughly the same size as the stellar
scale radius. Any smaller than this and we would gain logarithmically in the
range of impact parameters allowed for the scattering but lose like a power
law in the mass of the smaller clumps being considered. So the optimum is
pushed up very close to the maximum impact parameter $b_{\rm max}$.

###  II.4 Tidal heating rate of stars by clumps

In this subsection we calculate the tidal heat transfer rate from the clumps
to the stars. By tidal heating we mean scatterings where the clump passes
outside of the stellar scale radius and so their dominant effect on all the
stars is the same. Thus they primarily cause a change to the center-of-mass of
the whole stellar population, but not internal heating of the gas of stars. In
this subsection we calculate this suppressed heating rate.

We follow [33] and make a straight line approximation for clumps passing by
and causing tidal heating of the star system. If a clump with mass $M_{\rm
clump}$ passes by with impact parameter $b$ (distance from the center of the
stellar distribution), the energy injected into the star cluster by one pass
of a clump per unit stellar mass is [42],

$\displaystyle\frac{\langle\Delta E_{\rm
dt}(b)\rangle}{M_{s}}=\frac{4G^{2}M_{\rm clump}^{2}}{3v_{\rm rel}^{2}}\langle
r^{2}\rangle_{\star}A_{\rm corr}(b)\frac{\chi_{\rm st}(b)}{b^{4}},$ (15)

Here $\langle r^{2}\rangle_{\star}$ is evaluated over the star configuration
(defined in Ref. [33]), and $A_{\rm corr}$ is the adiabatic correction factor
since the stars are in orbit. $v_{\rm rel}$ is the relative velocity between
stars and clumps. The function $\chi_{\rm st}$ depends on the density profile
of the clump and is described in Ref. [42]. As a result the heating rate can
be determined as,

$\displaystyle H_{\star,\rm dt}=\int_{b_{\rm min}}^{b_{\rm max}}\Delta E_{\rm
dt}(b)\,\,2\pi b\,db\,\frac{\rho_{\rm clumps}(t)}{M_{\rm clump}}v_{\rm clump}$
(16)

Similar to the direct heating case, we drop out the radial dependence as well
as the time dependence of $\rho_{\rm clumps}$ via the replacement $\rho_{\rm
clumps}(t)=\rho_{\rm clumps}(t=0)\times\eta_{\rm NT}$, which is a good
approximation and simplifies the computation. The lower limit of the integral
is given by,

$\displaystyle b_{\rm min}=\textrm{Max}\left(b_{\rm samp},3R_{0,\star}\right)$
(17)

and $b_{\rm max}=R_{\rm UFD}$.

###  II.5 Limits without incorporating clump destruction effects

In this subsection we use the above stellar heating rates to set limits on
what clumps can exist inside a UFD. We use the UFD Segue-I with stellar scale
radius of 24.7 pc as an example. Roughly, we will call a point ruled out if
the clumps cause so much heating that in order to reproduce the currently
observed stellar scale radius, the stars initially needed to start in a tight
cluster of size less than 2 pc.

We solve the same differential equation derived in Ref. [33] for the evolution
of the stellar scale radius $R_{0,\star}$.

$\displaystyle\frac{dR_{0,\star}}{dt}=H_{\star}(R_{0,\star})\times\left\\{\frac{3\pi
GM_{\star}}{64R_{0,\star}^{2}}+\frac{GM_{\rm UFD}R_{0,\star}}{4R_{\rm
UFD}^{3}}\right.$ $\displaystyle\left.\left(6\log\frac{R_{\rm
UFD}}{R_{0,\star}}+(3\log(4)-17)\right)\right\\}^{-1}$ (18)

where $H_{\star}=H_{\star,{\rm direct}}+H_{\star,\rm dt}$. Setting
$R_{0,\star}(\rm today)=24.73~{}\textrm{pc}$ and $\eta=\eta_{\rm NT}$, we set
limits on the heating rate, and hence limits on the clumps by requiring that
$R_{0,\star}(t=0)>2~{}\textrm{pc}$. A parameter point that predicts
$R_{0,\star}(t=0)<2~{}\textrm{pc}$ is ruled out. We plot the inferred limits
in the $\rho_{s}$ vs $M_{\rm clump}$ plane in Fig. 2. It is important to
stress that these are not limits on this parameter space, since survival under
clump destruction effects, covered in Sec. III are not incorporated yet. We
display this figure solely to explain the dependence of the limits on the
density of the clumps.

We see that as the clump density $\rho_{s}$ increases towards
$10^{6}M_{\odot}~{}\textrm{pc}^{-3}$ we see that the boundaries become almost
vertical. On the right side the boundaries are exactly vertical and are set by
sampling effects as explained in [33]. On the left side the boundaries are
changing only logarithmically with the density. As the density increases at a
fixed mass the clumps become smaller. This increases the range of impact
parameters allowed in the log in the heating rate, thus slightly increasing
the heating rate and allowing us to rule out slightly lower masses. But for
the higher densities we consider, the low mass clumps are already quite small
and this effect is small. For lower densities we can see the lines on the left
side bending due to this logarithmic effect. If we kept plotting up to very
high densities this plot would of course reproduce the results for (point-
like) MACHOs/PBHs from [33]. Note that fractions of DM as low as $f_{\rm
clump}\sim 10^{-4}$ can be constrained.

Figure 2: This figure shows what the limits from stellar heating in the clump
density vs clump mass plane would have been if there were no effects that
caused destruction of clumps in UFDs. So note that these are NOT the real
limits, the real limits including destruction of clumps from tidal
heating/stripping effects are shown in Figure 5. We show curves for several
different fractions $f_{\rm clump}$.

##  III Clump survival in UFDs

In hierarchical structure formation, the clumps formed via collapse act as
particles when larger halos begin to collapse later. These clumps that are now
gravitationally bound to the UFD can undergo mass loss due to various
processes. We use the formalism developed in [42] to quantify clump survival
due to the different processes.

In all processes of interest, we will consider a clump “A” transiting another
object “B”. Here “B” could be individual stars, other individual clumps or the
host halo.

Let the particles of clump “A” have individual velocities $v_{i}$ and $\Delta
v_{i}$ be the velocity kick caused by the transit of “B”. Then, the energy
increase in clump “A” per unit mass is given by,

$\displaystyle\Delta E_{i}=v_{i}\cdot\Delta v_{i}+\frac{1}{2}(\Delta
v_{i})^{2}$ (19)

The first term averages to zero, but a diffusion term survives. After $N$
encounters,

$\displaystyle\Delta E_{i}(N)=\sqrt{N}v_{i}\cdot\langle\Delta
v_{i}\rangle+N\frac{1}{2}\langle\Delta v_{i}\rangle^{2}$ (20)

We will next show that it is justified to consider only the second term in
order to calculate subhalo harrasment. It was shown in [42] that $\Delta
E_{i}\geq E_{i}$ is needed for the clump to get appreciable tidal heating. To
be conservative we require that the total energy transferred $\Delta E_{\rm
tot}\equiv\sum_{i}\Delta E_{i}$ does not exceed the binding energy
$E_{b}\equiv|\sum_{i}E_{i}|$, i.e.

$\displaystyle E_{b}\geq\Delta E_{\rm tot}$ (21)

in order to avoid disruption. Right on the edge of our disruption condition,
when $\Delta E_{i}\sim E_{i}$, the two terms in equation 20 are comparable.
For small energy transfer $\Delta E_{i}$ and small $N$, the first term
dominates. The second term starts becoming important when $\Delta E_{i}$ is
comparable to the initial energy per mass, $E_{i}\sim v_{i}^{2}\approx N\Delta
v_{i}^{2}\sim\Delta E_{i}$. After this, the second term dominates. Since Ref.
[42] found that $E_{b}\gg\Delta E_{\rm tot}$ for significant mass loss in
clumps. In this limit the second term s the relevant one and hence we consider
only the second term going forward.

Next let us estimate both the LHS and RHS of Eqn. 21. The LHS is given by,

$\displaystyle E_{b}=4\pi G\int drr\rho(r)M_{\rm enc}(r)$ (22)

The RHS, for a given impact parameter $b$, is given by [42],

$\displaystyle\langle\Delta E(b)\rangle=\frac{4G^{2}M_{B}^{2}M_{A}}{3v_{\rm
rel}^{2}}\langle r^{2}\rangle A_{\rm corr}\frac{\chi_{\rm st}(b)}{b^{4}}$ (23)

Here $\langle r^{2}\rangle$ is the density averaged $r^{2}$ over the clump
“A”. The function $\chi_{\rm st}$ is also defined in [42] and it captures the
density profile of object “B”. It is smaller than unity for an extended mass
distribution and one for a point mass. Here $A_{\rm corr}$ is the adiabatic
correction factor to account for the dilution in tidal heating when the
orbital frequency of particles in “A” about “A”’s center is larger than the
time for “A” to transit “B”. We describe the particular form of $A_{\rm corr}$
in each individual subsection.

In the rest of this section, we adapt this formalism for different choices of
“B”.

###  III.1 Host Halo Heating

We discuss the heating of clumps due to the gravitational potential of the
host halo in this subsection. We identify “A” with the clumps and “B” with the
host halo.

For a host halo that satisfies a Dehnen profile with $\gamma=0$, we can
simplify the function $\chi_{\rm st}$ that captures the profile information of
the host in the $b\ll R_{\rm UFD}$ limit. Taking the expression $\chi_{\rm
st}$ from [42], we find in the $b\ll R_{\rm UFD}$ limit,

$\displaystyle\chi_{\rm st}(b)=\frac{b^{4}}{4R_{\rm UFD}^{4}}.$ (24)

Furthermore, for NFW cores, $\langle r^{2}\rangle=0.41r_{s}^{2}$. Thus,

$\displaystyle\langle\Delta E_{\rm HHH}\rangle=0.41\frac{4G^{2}M_{\rm
UFD}^{2}M_{\rm clump}}{3v_{\rm rel}^{2}}r_{s}^{2}\frac{A_{\rm
corr,HHH}}{4R_{\rm UFD}^{4}}$ (25)

Here the subscript HHH stands for host halo heating.

To our knowledge, thus far, the adiabatic correction for host halo heating has
only been derived in the $b\approx R_{\rm host}$ limit, i.e. when the typical
pericenter of the clump orbit is comparable to the scale radius of the host
halo. In our case, we are interested in the survival of clumps that pass near
the stars in the UFD, rather than the typical clump, since it is these clumps
which are responsible for the heating of the stars. Thus we are interested in
the case $b\approx R_{0,\star}\ll R_{\rm host}=R_{\rm UFD}$ where
$R_{0,\star}$ is the stellar scale radius of the UFD111Our discussion holds
even in the case where pericenter changes from orbit to orbit. This is because
as we show later the pericenter drops out of our calculations entirely. .

The adiabatic correction kicks in when the timescale of scattering with the
host halo $\tau$ is longer than the orbital period of particles about the
clump $\frac{1}{\omega}$. $\tau$ is typically taken as $\tau=\frac{b}{v_{\rm
rel}}$ where $b$ is the impact parameter. However, since the clump spends all
its orbit inside the UFD, it is possible that $\tau$ is much larger,
$\mathcal{O}\left(\frac{R_{\rm clumps}}{v_{\rm rel}}\right)$. Here $R_{\rm
clumps}(t)$ is the scale radius of the clump distribution which evolves in
time since the clumps feel dynamical friction. To estimate this, we construct
a toy model of a mass on a spring that undergoes simple harmonic motion with
frequency $\omega$. This approximately captures the orbit of individual
particles about the clump center. We then integrate the tidal force as the
center of this harmonic oscillator traverses the UFD with pericenter $b$. For
the $\gamma=0$ Dehnen profile, we confirm that the adiabatic correction is
indeed very weakly dependent on $b$ but instead depends on $R_{\rm clumps}$ in
the $b\ll R_{\rm clumps}$ limit. Hence we take the adiabatic factor to be

$\displaystyle A_{\rm corr,HHH}=\left(1+\frac{v_{\rm
orb}^{2}}{r_{s}^{2}}\frac{R_{\rm clumps}^{2}}{v_{\rm
rel}^{2}}\right)^{-\frac{3}{2}}$ (26)

Here $v_{\rm orb}=\sqrt{\frac{0.45GM_{\rm clump}}{r_{s}}}$ is the typical
orbital velocity of a particle about the clump Thus, the heating due to one
passage given in Eqn. 25 is independent of the pericenter $b$.

Hence, the instantaneous heating rate for an individual clump is given by

$\displaystyle H_{\rm HHH}(t)=\langle\Delta E_{\rm HHH}(R_{\rm
clumps}(t))\rangle\frac{v_{\rm clumps}}{R_{\rm clumps}(t)}$ (27)

$R_{\rm clumps}(t)$ is obtained by solving the migration differential Eqn. 9.

The total amount of energy imparted after time $t$ is given by,

$\displaystyle\Delta E_{\rm HHH}^{\rm tot}(t)=\int^{t}_{0}H_{\rm
HHH}(t^{\prime})dt^{\prime}$ (28)

Using this quantity we can estimate the time taken for a clump to be destroyed
by host halo heating, $t_{\rm HHH}$. This is the time taken for the total
energy imparted to equal the binding energy, i.e.

$\displaystyle\Delta E_{\rm HHH}^{\rm tot}(t=t_{\rm HHH})=E_{b}$ (29)

We next define an effective density enhancement after incorporating host halo
heating $\eta_{\rm HHH}$,

$\displaystyle\eta_{\rm HHH}=\frac{1}{t_{\rm UFD}}\int_{0}^{T_{\rm
HHH}}dt\frac{\rho_{\rm clumps}(R_{0,\star},t)}{\rho_{\rm clumps}(R,t=0)}$ (30)

the difference between this equation and Eqn. 10 is that the upper limit of
the integral is now set to $t_{\rm HHH}$. We generally want to know the
density enhancement (or reduction) in the location of the stars, so we
evaluate this at the scale radius of the stars. However once again, we can
actually evaluate this anywhere inside the clumps scale radius so we replace
$\rho_{\rm clumps}(R_{0,\star},t)$ with its approximate value derived in Eqn.
(11).

$\eta_{\rm HHH}$ is plotted in Fig. 1 as a function of the clump mass $M_{\rm
clump}$ for different choices of the density of the clump $\rho_{s}$ for
$f_{\rm clump}=10^{-3}$ in solid lines. One can see that for core densities
down to $\rho_{s}\approx 10^{3}M_{\odot}\textrm{pc}^{-3}$, $\eta_{\rm HHH}$ is
identical to $\eta_{\rm NT}$ (dashed lines), i.e. the effect of host halo
heating is negligible. However there is a large deviation in $\eta_{\rm HHH}$
from $\eta_{\rm NT}$ for $\rho_{s}=100\,M_{\odot}\,\textrm{pc}^{-3}$ for
masses $3\times 10^{4}\,M_{\odot}$ and above. This can be understood as the
effect of migration, which reduces $R_{\rm clumps}$ which in turn exacerbates
host halo heating. Nonetheless $\eta_{\rm HHH}$ is still larger than unity.
For $\rho_{s}=10\,M_{\odot}\,\textrm{pc}^{-3}$, the density is so low that
$t_{\rm HHH}\ll t_{\rm UFD}$ even for small clump masses where migration is
negligible. As a result $\eta_{\rm HHH}$ is smaller than unity for these
densities, meaning that the host halo heating effect destroys the clumps fast
enough that this effect is larger than the migration effect and thus there is
less total heating.

###  III.2 Clump-Clump Harassment

The issue of subhalo-subhalo harassment was treated in detail in [42]. Hence
our discussion mirrors the one found in the above reference. Specifically we
identify both “A” and “B” with clumps in this subsection. Therefore the energy
imparted per clump-clump collision is given by,

$\displaystyle\langle\Delta E_{\rm CCH}(b)\rangle=0.41\frac{4G^{2}M_{\rm
clump}^{3}}{3v_{\rm rel}^{2}}r_{s}^{2}A_{\rm corr,CCH}\frac{\chi_{\rm
st}(b)}{b^{4}}$ (31)

Here, the subscript CCH stands for clump-clump harassment. Here,

$\displaystyle A_{\rm corr,CCH}=\left(1+\frac{v_{\rm
orb}^{2}}{r_{s}^{2}}\frac{b^{2}}{v_{\rm rel}^{2}}\right)^{-\frac{3}{2}}$ (32)

The total energy imparted over the age of the UFD is then,

$\displaystyle\Delta E_{\rm tot,CCH}=\int_{b_{\rm min}}^{b_{\rm max}}$
$\displaystyle\langle\Delta E_{\rm CCH}(b)\rangle 2\pi bdb$
$\displaystyle\times\eta_{\rm HHH}f_{\rm clump}\frac{\rho_{\rm DM}}{M_{\rm
clump}}T_{\rm UFD}v_{\rm rel}$ (33)

Note that we use the enhancement $\eta_{\rm HHH}$ that was evaluated at the
scale radius of the stars which was shown to be a good approximation for any
radius smaller than $R_{\rm clumps}$. Since a majority of the clumps sit
inside $R_{\rm clumps}$ as well, $\eta_{\rm HHH}$ is an excellent
approximation for the enhancement in clumps that any individual clump sees.
Here $b_{\rm min}$ is given by solving $\pi b_{\rm min}^{2}\eta_{\rm
HHH}f_{\rm clump}\frac{\rho_{\rm DM}}{M_{\rm clump}}v_{\rm rel}=1$ and $b_{\rm
max}=R_{\rm UFD}$. Unlike the other processes discussed thus far, subhalo
harassment depends on the fraction of DM in clumps i.e. $f_{\rm clump}$.
Afterall, $\Delta E_{\rm tot,CCH}=0$ for $\eta_{\rm HHH}f_{\rm clump}=0$ and
monotonically increases with $\eta_{\rm HHH}f_{\rm clump}$. Hence a $\eta_{\rm
HHH}f_{\rm clump}$ exists for which $\Delta E_{\rm tot,CCH}=E_{b}$. We define
this quantity as $f_{\rm clump}^{\rm safe}$, i.e. the maximum effective DM
fraction $\eta_{\rm HHH}f_{\rm clump}$ that is safe from clump-clump
harassment. We plot contours of $f_{\rm clump}^{\rm safe}$ in the $\rho_{s}$
vs $M_{\rm clump}$ plane in Fig. 3. Scale densities $\rho_{s}$ above a
particular $f_{\rm clump}^{\rm safe}$ contour are safe from harassment for
$f_{\rm clump}^{\rm safe}=\eta_{\rm HHH}f_{\rm clump}$. In the small $M_{\rm
clump}$ limit, the heating rate is proportional to $M_{\rm clump}$ and hence
larger $\rho_{s}$ are required to keep the clumps harassment free for larger
masses. However, at some large mass that is $f_{\rm clump}$ dependent, the low
sampling rate leads to large impact parameters which in turn leads to a large
adiabatic correction. This dilutes the clumps harassment rate which leads to
the turn around behavior. We see that for $f_{\rm clump}^{\rm safe}=0.01$
scale densities $\rho_{s}=0.01~{}M_{\odot}\textrm{pc}^{-3}$ are safe from
harassment for the entire range of clump masses we consider.

Clump-clump harassment is incorporated by taking the minimum of $\eta_{\rm
HHH}f_{\rm clump}$ and $f_{\rm clump}^{\rm safe}$ for any parameter point
$\\{\rho_{s},M_{\rm clump}\\}$, i.e the effective clump fraction we take never
exceeds the safe fraction with regard to clump-clump harassment. The final
result of this section is that we have derived a replacement rule

$\displaystyle\rho_{\rm clumps}(t)\rightarrow\rho_{\rm
DM}\times\textrm{Min}\left(\eta_{\rm HHH}f_{\rm clump},f_{\rm clump}^{\rm
safe}\right)$ (34)

which should be substituted in Eqn. (13) for direct heating and Eqn.(16) for
tidal heating222Note that the enhancement factor $\eta_{\rm HHH}$ technically
applies only inside the clumps’ scale radius while the tidal heating comes
from clumps at all distances. However the contribution to the tidal heating
coming from clumps outside the tidal radius is negligible here, partly because
everywhere that tidal heating is even relevant (see Fig. 5 Right of [33]) the
clumps’ tidal radius remains always much bigger than 24.7 pc. respectively.

Figure 3: Contours of $f_{\rm clump}^{\rm safe}$ are plotted in the $\rho_{s}$
vs $M_{\rm clump}$ plane. Scale densities $\rho_{s}$ above a particular
$f_{\rm clump}^{\rm safe}$ contour are safe from clump-clump harassment for
$f_{\rm clump}^{\rm safe}=\eta_{\rm HHH}f_{\rm clump}$.

We also calculated harassment from individual stars but it is far subdominant
to all other clump destruction processes so we do not discuss it further.

###  III.3 Effects of clump destruction on constraints

In this subsection we show the effects that these clump destruction processes
have on our constraints. Without these processes we would have set the
constraints in Figure 2. However these destruction processes cut off the
curves below some density.

We demonstrate the effect of these clump destruction processes for a single
value of $f_{\rm clump}$ in Fig. 4.

First, the constraint that would be derived if all the clumps were deemed to
survive without incorporating tidal heating processes is shown in dashed red.
This is the same red dashed line from Fig. 2. Next, we show the contour after
incorporating host halo heating alone in dotted line (labeled HHH only). We
see that this roughly cuts of the constraint at $\rho_{s}\approx
10M_{\odot}\textrm{pc}^{-3}$. Finally, incorporating both host halo heating as
well as clump-clump harassment (CCH) produces the solid red curve. We see that
CCH is not a significant process as compared to HHH for $M_{\rm clump}\gtrsim
10^{4.5}M_{\odot}$. This is the mass range where host halo heating
significantly affects $\eta_{\rm HHH}$ as seen in Fig. 1 owing to significant
migration. Coincidentally, it is also the mass range where CCH peaks as seen
in Fig. 3. For smaller masses where CCH dominates, we see a further dent in
constraints for “HHH+CCH” as compared to the “HHH only” curve. We also replot
the $f_{\rm safe}^{\rm clump}=0.1$ contour from Fig. 3 to illustrate the fact
that CCH is an issue only below the black contour. This is the reason why the
“no HHH or CCH” (dashed) curve and the “HHH + CCH” (solid) curve share a
common intersection point with the $f_{\rm safe}^{\rm clump}=0.1$ curve
(black). The solid lines after incorporating both host halo heating as well as
clump clump harassment are the main result of this paper. We discuss only
these contours moving forward. Below we discuss contours for many values of
$f_{\rm clump}$.

Figure 4: Limits on clump parameter space after incorporating various clump
destruction processes, taking $f_{\rm clump}=0.1$ as an example. The “no HHH
or CCH” dashed line contains no clump destruction and is identical to Fig. 2.
The “HHH only” dotted line incorporates only host halo harassment. The
“HHH+CCH” solid line additionally incorporates clump-clump harassment. The
$f_{\rm clump}^{\rm safe}=0.1$ black line is identical to the red line in Fig.
3 and reproduced here for reference. The solid red line heres is reproduced as
the solid red line in Fig. 5 along with contours for other $f_{\rm clump}$.

##  IV Results

In this section we present our results for constraints on the DM
subhalos/clumps in a UFD arising from their heating effect on the stellar
population. We use the heating rate calculation given in Section II to find
the parts of clump parameter space (clump mass, density, and dark matter
fraction) which produce excessive heating, assuming such clumps exist in the
UFD. This produces the region shown in Figure 2 where, if such clumps existed
in the UFD, they would have caused too much heating of the stars. We then
combine this with the calculations of Section III which determine the minimum
density necessary for a clump to survive destruction in the environment of the
UFD (as a function of clump mass and fraction of dark matter). This puts a
lower bound on the regions of Figure 2 where clumps could exist and thus cuts
off our constraints at the lowest densities. Combining these gives our main
result shown in Figure 5.

Figure 5: Our constraints on clump parameter space in clump mass $M_{\rm
clump}$, clump scale density $\rho_{s}$, and fraction of dark matter made up
of clumps $f_{\rm clump}$. In the left plot, the different lines represent
contours of different fraction $f_{\rm clump}$ in the $\rho_{s}$ \- $M_{\rm
clump}$ plane. In the right plot the different lines represent contours of
different $\rho_{s}$ in the $f_{\rm clump}$ \- $M_{\rm clump}$ plane.

Figure 6:  Our limits (blue) on the power spectrum of adiabatic primordial
density perturbations $P(k)$ vs the wave-number $k$. These are significantly
stronger than existing bounds in this range from FIRAS (red) via
$\mu$-distortions and y-distortions [43]. At longer length scales the power is
measured by CMB (green) [44] and Lyman-$\alpha$ (orange) [45]. For comparison
we also show the nearly scale invariant power measured by Planck extrapolated
to higher $k$ (black-dashed). The blue shaded region corresponds to the top-
hat window function with the Press-Schechter [46] multiplicity function and
$\Delta=200$. In the Left panel, we illustrate the uncertainty associated with
a choice of filter, the Gaussian (dot-dashed blue), as well as an alternate
choice of $\Delta=600$. In the Right panel, we show alternate choices of the
multiplicity function: Sheth-Tormen [47] (Dashed) and Watson et al. [48] (Dot-
dashed).

###  IV.1 Limits on clumps in UFDs

Figure 5 shows our constraints on the clump parameter space in the three
variables: clump mass $M_{\rm clump}$, clump density $\rho_{s}$, and fraction
of dark matter made up of these clumps $f_{\rm clump}$.

The left plot of Figure 5 appears similar to Figure 2 except the limits are
cut off below a certain density because clumps of such low densities are
destroyed by tidal forces in the UFD before they can heat the stars
excessively. For masses above $\sim 10^{5}\,M_{\odot}$ the lower limit on the
density that we can constrain comes from host halo heating. For masses below
$\sim 10^{5}\,M_{\odot}$ the lower limit arises from clump-clump harassment.
The constraints cut off sharply at high masses (the vertical edges on the
right side of the curves) because of sampling; past those edges there simply
are not any clumps of that size in the UFD. The left edge of each curve at low
masses comes from the heating rate calculation, at masses below those edges
such clumps would not heat the stars significantly.

The right plot of Figure 5 is similar to the main results in Figure 6 of Paper
1 [33]. Thus we see that for clumps of high enough density (above $\sim
10^{6}\,M_{\odot}\,\text{pc}^{-3}$) there is essentially no difference from
the point mass (MACHO/PBH) case. Such dense clumps do not get destroyed by the
environment of the UFD. And the fact that they are spatially extended (much
more so than a PBH for example) makes almost no difference to the heating rate
since gravitational scattering is inherently long-range and the only loss is a
small reduction in the log (see Eqns. (13) and (14) for example). As
$\rho_{s}$ is dialed down, we see there is an abrupt drop in limit at a small
enough $M_{\rm clump}$ . This is due to clump-clump harassment and occurs when
the contour hits $f_{\rm clump}=f_{\rm safe}$. For much smaller densities
$\rho_{s}\approx 10~{}M_{\odot}~{}\textrm{pc}^{-3}$, we finally see the effect
of host halo heating which further suppresses the constraints.

###  IV.2 Primordial Power Spectrum

We can now translate our limits on the existence of clumps in UFDs into limits
on the primordial power spectrum of density perturbations from inflation. The
primordial power spectrum determines the distribution of masses and densities
of dark matter structure that forms in the early universe. There are of course
significant uncertainties in predicting the spectrum of clumps in a UFD that
will arise from a given primordial power spectrum of density perturbations.
For simplicity we will take an analytic procedure for estimating this, and try
to keep our estimated limits on the primordial power spectrum as conservative
as possible.

Current observations are consistent with a nearly scale invariant primordial
power spectrum upto comoving wave-number $k\approx 3\,\textrm{Mpc}^{-1}$

$\displaystyle P_{\rm CMB}(k)=P(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}$
(35)

with $\log\left(10^{10}P(k_{*})\right)=3.047\pm 0.014$ and $n_{s}=0.966\pm
0.004$ with $k_{*}=0.05~{}\textrm{Mpc}^{-1}$ obtained from Ref. [49]. Limits
on low mass halos from strong lensing was used to constrain power down to the
$10~{}\textrm{Mpc}^{-1}$ scale in Ref. [50]. There are only very weak limits
$P_{k}\approx 10^{-4}$ at larger scales arising from $y$ and $\mu$ distortions
of the CMB as measured by FIRAS [43].

Here we translate the derived limits on clumps onto limits on the primordial
power spectrum at higher scales. We choose to place constraints on a power
spectrum which is the nearly flat spectrum given above plus a spike in the
spectrum at some $k$ which is higher than the range of observations.

So we model the excess power on small scales as a delta function at scale
$k_{\rm high}$

$\displaystyle P(k)=P_{\rm CMB}(k)+P_{k}k_{\rm high}\delta(k-k_{\rm high})$
(36)

We will then place limits on the amount of power at short scales in the
$P_{k}$ vs $k_{\rm high}$ parameter space.

Given a primordial power spectrum $P(k)$, we can obtain the linear matter
spectrum $P_{\rm lin}(k)$ via,

$\displaystyle P_{\rm
lin}(k)=\frac{1}{4}2\pi^{2}P(k)T(k)^{2}\left(\frac{k}{k_{\rm eq}}\right)^{4}$
(37)

Here, $k_{\rm eq}$ is the wave-number of the perturbation mode that enters the
horizon at matter-radiation equality and $T(k)$ is the transfer function,
which we take from Eisenstein and Hu [51]. We can next calculate the variance,

$\displaystyle\sigma^{2}(R,z)=\frac{D^{2}(z)}{2\pi^{2}}\int k^{2}P_{\rm
lin}(k)\tilde{W}^{2}(kR)dk$ (38)

where $D(z)$, is the growth function, which we obtain from Ref. [52] and
$\tilde{W}(kR)$ is the filter function and $R$ is the comoving length scale.

The halo mass function at a redshift $z$ is given by,

$\displaystyle\frac{dn_{\rm clump}(z)}{d\ln M_{\rm
clump}}=f(\sigma)\frac{\rho_{M}(z)}{M_{\rm clump}}\frac{d\ln\sigma^{-1}}{d\ln
M_{\rm clump}}$ (39)

Here, $n_{\rm clump}$ is the number density of clumps, with clump mass $M_{\rm
clump}$ and $\rho_{M}(z)$ is the non-relativistic matter density at redshift
$z$. Since $R$ is a comoving scale, we set $M_{\rm clump}=\frac{4}{3}\pi
R^{3}\rho_{\rm crit}\Omega_{m}$, where $\rho_{\rm crit}$ is the energy density
of the universe today and $\Omega_{m}$ is the matter fraction today. This can
be used to write $\sigma$ as a function of $M_{\rm clump}$. The multiplicity
function $f(\sigma)$ is chosen to be that given by Press-Schechter [46] for
limit setting purposes. We also provide limits with other choices such as
Sheth-Tormen [47] and Watson et al [48] in order to show how much our limits
depend on this choice.

Now that we have the number density given a primordial power spectrum, we can
compare it to the limits derived in the previous section, to finally obtain
the limits on primordial power. Our procedure for setting this limit is as
follows. Given a primordial power spectrum $P(k)$, we obtain the fraction of
DM that has collapsed into clumps with masses within a decade of $M_{\rm
clump}$ mass by redshift $z$ as

$\displaystyle f_{\textrm{clump},P(k)}(M_{\rm clump},z)\approx M_{\rm
clump}\frac{dn_{\rm clump}(z)}{d\ln M_{\rm clump}}$ (40)

For a clump that collapsed at redshift $z$, its average density is given by

$\displaystyle\langle\rho_{\rm clump}\rangle(z)\approx\Delta\rho_{\rm
crit}\Omega_{m}(1+z)^{3}$ (41)

Here $\Delta\approx 200$ is the resultant over density above the ambient
matter density due to gravitational collapse. If we assume that this results
in NFW cores, then the density at the scale radius of the clump is

$\displaystyle\rho_{s}(z)\approx\frac{\Delta}{2.317}\rho_{\rm
crit}\Omega_{m}(1+z)^{3}$ (42)

The limits we set in Fig. 5 give a maximal allowed fraction $f$ as a function
of mass and density, namely $f_{\rm clump,lim}\left(M_{\rm
clump},\rho_{s}\right)$. Thus, we consider a $P(k)$ distribution ruled out if
there exists any $z$ and any $M_{\rm clump}$ for which,

$\displaystyle f_{\textrm{clump},P(k)}\left(M_{\rm clump},z\right)>f_{\rm
clump,lim}\left(M_{\rm clump},\rho_{s}(z)\right)$ (43)

So to decide whether a particular $P(k)$ is ruled out, we scan over all
$z$333Actually we only include redshifts that are before the time of formation
of a UFD which we assume to be around $z\sim 6-10$. However this constraint
does not affect our result at all because we will only use much higher
redshifts in order to get to the clump densities that our results can
constrain. and all $M_{\rm clump}$.

If a particular $P(k)$ from Eqn. (36) is ruled out, then we will say that
value of $P_{k}$ is ruled out at the scale $k=k_{\rm high}$.

We plot the limits on $P_{k}$ vs $k$ obtained this way as the blue curve in
Fig. 6.

Our limits on the power spectrum of primordial density perturbations $P_{k}$
vs the wave-number $k$ are more than two orders of magnitude stronger than
existing bounds in this range from FIRAS (red) via $\mu$-distortions and
y-distortions [43]. They extend beyond wavenumbers
$k\approx\textrm{kpc}^{-1}$. At longer length scales, the power is measured by
CMB (green) [44] and Ly$\alpha$ (orange) [45]. For comparison we also show
$P_{\rm CMB}$ from Eqn. 35 i.e. the nearly scale-invariant Planck constrained
power extrapolated to much higher $k$ in black-dashed. The blue shaded region
corresponds to the top-hat window function with the Press-Schechter
multiplicity function and $\Delta=200$.

In the left panel, we illustrate the uncertainty associated with the choice of
filter. A Gaussian filter (dot-dashed blue) shifts the limits somewhat to
lower $k$. While we have chosen $\Delta=200$ to be conservative, far higher
values of $\Delta$ are found in literature [30, 16]. To illustrate the
constraining power for a larger $\Delta$, we show the limit for a top-hat
window and $\Delta=1390$ (which corresponds to $C_{\rho}=600$ in Ref. [30]) in
dashed lines. As expected, this puts deeper constraints on the primordial
power.

In the right panel, we show alternate choices of the multiplicity function:
Sheth-Tormen (Dashed) and Watson13 (Dot-dashed). As seen, the constraints are
fairly independent of the choice of multiplicity function.

###  IV.3 Isocurvature Perturbations

Figure 7: Our limits (blue) on the power spectrum of dark matter isocurvature
primordial density perturbations $P_{k}$ vs the wave-number $k$. These are the
only bounds we know in this range of scales. Also shown are the PIXIE
projections via $\mu$-distortions obtained by appropriately rescaling the
projections for the adiabatic power spectrum found in [43] (Refer to text for
more details). At longer length scales the isocurvature power cannot exceed
the power measured in the CMB by Planck (green) [44] and Lyman-$\alpha$
(orange) [45]. At even longer scales, the limits on isocurvature set by Planck
Ref. [44] are shown (dark green).

Our limits apply not just to adiabatic density perturbations as discussed
above, but also to isocurvature perturbations in the dark matter density
alone. Anything that causes dark matter to form structures at the relevant
scales $\sim 10-10^{8}\,M_{\odot}$ will be constrained by our approach. We
show our limit on the power spectrum of dark matter isocurvature perturbations
$P_{\rm DM,iso}$ in Figure 7. Again we have taken a delta function spike at
$k$ for the DM isocurvature power spectrum to make this plot. The solid blue
line here for limits on the isocurvature power spectrum is identical to the
solid blue line in Fig. 6 for the adiabatic power spectrum.

We also show existing constraints on the dark matter isocurvature power
spectrum arising from Planck [44] (dark green), which require dark matter
isocurvature to be no more than $4\%$ of the power in the adiabatic density
perturbations at those scales. Further, as a rough approximation we assume the
dark matter isocurvature perturbations cannot be larger than the measurements
of the total power spectrum in adiabatic perturbations from CMB (light green)
[44] and Ly$\alpha$ (orange) [45].

Note that the FIRAS constraints on the adiabatic perturbations do not appear
in Figure 7. The FIRAS constraints are actually constraints on the photon
power spectrum $P_{\gamma}$. Dark matter isocurvature perturbations will alter
the photon power spectrum at horizon reentry on the corresponding scales, but
this is suppressed by $P_{\gamma}(k\approx aH)\sim\mathcal{O}(10)\frac{k_{\rm
eq}^{2}}{k^{2}}P_{\rm DM,iso}$ [53] where $k_{\rm eq}\sim 0.01\,{\rm
Mpc}^{-1}$ is the comoving scale of matter-radiation equality. Thus the FIRAS
constraint is suppressed significantly and does not appear on the plot. A
projection for the sensitivity of PIXIE to dark matter isocurvature
perturbations [54] is shown in Figure 7 which we derived using the same
suppression factor as for FIRAS.

We are not aware of any other constraints besides ours on dark matter
isocurvature perturbations in the same range of scales as our limits.

Figure 8: Our limit (blue) on a power law isocurvature dark matter power
spectrum $P_{\rm DM,iso}=A_{\rm iso}\left(\frac{k}{k_{0}}\right)^{n_{\rm
iso}-1}$ is shown in the $A_{\rm iso}$ vs $n_{\rm iso}$ plane. We calculate
existing limits from Planck’s isocurvature measurement (green) and
Lyman-$\alpha$ (orange), see text for more details. Limits from FIRAS (red)
[54] as well as projections from PIXIE [54] are also shown.

We can also ask how our results constrain a power law spectrum for the DM
isocurvature perturbations rather than a delta function spike at a given $k$.
We take the parametrization of [54]

$\displaystyle P_{\rm DM,iso}=A_{\rm iso}\left(\frac{k}{k_{0}}\right)^{n_{\rm
iso}-1}$ (44)

where $k_{0}=0.002\,{\rm Mpc}^{-1}$. We show limits on $A_{\rm iso}$ as a
function of $n_{\rm iso}$ in Figure 8. Our limits are orders of magnitude
stronger than existing constraints, which arise dominantly from Ly$\alpha$
constraints (orange) at lower $k$. This previous constraint is obtained by
requiring that $P_{\rm DM,iso}(k)$ is smaller than the Ly$\alpha$ limit at
$k\approx 3~{}\textrm{Mpc}^{-1}$. A similar procedure but invoking the
isocurvature limits from Planck [44] at $k=0.1~{}\textrm{Mpc}^{-1}$ produces
the green line. Also shown are limits from FIRAS and projections from PIXIE as
derived in Ref. [54].

Note that when $n_{\rm iso}=4$ we have a power spectrum $\propto k^{3}$ and we
set a limit $A_{\rm iso}=8\times 10^{-24}$. This is a universal, white noise
form that arises in many production mechanisms for ultralight dark matter (see
e.g. [15, 53, 55]). Our results give the strongest constraints on this generic
form of the DM isocurvature spectrum which applies to many models.

##  V Discussion

We have placed limits on the existence of dark matter subhalos by the heating
effect they would have on the stars in an ultrafaint dwarf galaxy, see Figure
5. This bounds subhalos of mass between about $10\,M_{\odot}$ and
$10^{8}\,M_{\odot}$. There are many ramifications of this bound. For example
we demonstrated that this places limits on the primordial power spectrum over
a range of higher $k$ (shorter distances) than CMB or Ly$\alpha$ measurements
can reach, see Figure 6. These limits are orders of magnitude stronger than
existing limits in this range. And at $k\sim 100\,{\rm Mpc}^{-1}$ the limit we
place currently with UFDs is almost as strong as the projected limit that the
PIXIE mission would be able to place (see e.g. [56, 43]). Further, we also
constrain dark matter isocurvature perturbations, giving the strongest limit
in this range of $k$, see Figure 7. Many models produce a rising power-law
spectrum of dark matter isocurvature perturbations. For such models our
constraint is stronger than the Planck or Ly$\alpha$ constraints on
isocurvature, see Figure 8.

Although these limits are not currently able to constrain a straight-line
extrapolation from lower-k CMB measurements, it may be possible to improve the
reach of this technique in the future. The lowest power that we can constrain
(Figure 6) is being limited by the lowest density clumps that we can constrain
(in Figure 5 Left). This lower edge of the density constraints is primarily
limited by what density clumps actually survive in the environment of the UFD.
In order to estimate this we have made multiple conservative assumptions. Thus
it is possible that a simulation of subhalos inside a UFD could improve these
bounds. Further, considering host objects other than UFDs might improve these
bounds since the high densities of the UFD can lead to greater destruction of
clumps.

Additionally, as discussed in [33], it might even be possible to look for
positive signals of this heating effect from clumps on stars, thus possibly
providing a way to actually detect the presence of dark matter subhalos at
masses below $\sim 10^{8}\,M_{\odot}$. This would require a large sample of
UFDs and more accurate predictions for the effect of such heating on the
stellar distributions as a function of the various parameters of the UFD
(density profiles, size, age, etc.).

Even without reaching the straight-line extrapolation, this significant
improvement in the limits on the primordial power spectrum (or really on the
existance of subhalos) can have important consequences.

For example, it was recently pointed out that our current knowledge of the
power spectrum puts a lower limit on the possible mass of the dark matter
particle, under certain assumptions about the DM formation mechanism [55]. Our
bounds are at significantly higher $k$ and thus likely significantly
strengthen this bound on the mass of dark matter.

And of course, our new bounds limit certain models of inflation which would
generate additional power at these higher $k$ (there are many papers on
inflationary models which generate higher power at high $k$ but for some
examples see [1, 2, 3, 4, 5, 6]).

Further, even if the initial power spectrum from inflation is flat, non-
trivial dynamics of dark matter could generate more substructure than expected
for standard CDM (see e.g. [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
23, 24, 25, 26]). For example axion dark matter can form large, dense clumps
[14], so some axion DM parameter space can be explored/constrained using our
results. Or for example self-interacting dark matter (SIDM) is known to change
dark matter structures on short scales. SIDM can greatly increase the density
of the cores of dark matter subhalos (see e.g. [17, 19, 20, 18, 21]). These
cores would then be much more likely to survive the various tidal destruction
forces at work in the UFD environment. Thus, our limits could put strong
constraints on SIDM. We leave consideration of this for future work.

In [33] we applied this heating effect to set limits on effectively point-like
objects (MACHOs/PBHs). In this paper we applied the effect to find limits on
gravitationally-bound clumps of any kind of dark matter. However, our results
will also apply to extended dark objects bound together by forces stronger
than gravity such as axion halos (see e.g. [14, 57, 41]). And in fact our
results will be stronger when applied to such objects, assuming the binding
force is stronger than gravity, since the main effect limiting our results was
the tidal destruction of the dark matter clumps. The extended size of the
object does not decrease the heating effect significantly since it only
affects the log as discussed above. For objects that do not get tidally
disrupted, our limits would become as strong as the curves in Figure 2 rather
than the curves in Figure 5 (left). Our bounds are complementary to existing
bounds (see e.g. [57]). Previously there was a gap in the limits on extended
dark objects making up all of dark matter for masses above about
$10\,M_{\odot}$ which our bounds close [57].

Studying structures smaller than currently observable dwarf galaxies is an
exciting way to learn more about cosmology and dark matter. We believe that
studying the heating effect such structures have on stars is one very
promising technique for learning more about such small-scale structure in our
universe.

###### Acknowledgements.

We thank M. Amin, Z. Bogorad, A. Ireland, D.E. Kaplan, M. Kaplinghat, P.
Mansfield, O. Slone, J. Thompson, R. Wechsler for useful discussions. The
authors acknowledge support by NSF Grants PHY-2310429 and PHY-2014215, Simons
Investigator Award No. 824870, DOE HEP QuantISED award #100495, the Gordon and
Betty Moore Foundation Grant GBMF7946, and the U.S. Department of Energy
(DOE), Office of Science, National Quantum Information Science Research
Centers, Superconducting Quantum Materials and Systems Center (SQMS) under
contract No. DEAC02-07CH11359.

## References

  * Leach  _et al._ [2001] S. M. Leach, M. Sasaki, D. Wands, and A. R. Liddle, Enhancement of superhorizon scale inflationary curvature perturbations, Phys. Rev. D 64, 023512 (2001), arXiv:astro-ph/0101406 .
  * Pajer and Peloso [2013] E. Pajer and M. Peloso, A review of Axion Inflation in the era of Planck, Class. Quant. Grav. 30, 214002 (2013), arXiv:1305.3557 [hep-th] .
  * Hertzberg and Yamada [2018] M. P. Hertzberg and M. Yamada, Primordial Black Holes from Polynomial Potentials in Single Field Inflation, Phys. Rev. D 97, 083509 (2018), arXiv:1712.09750 [astro-ph.CO] .
  * Mishra and Sahni [2020] S. S. Mishra and V. Sahni, Primordial Black Holes from a tiny bump/dip in the Inflaton potential, JCAP 04, 007, arXiv:1911.00057 [gr-qc] .
  * Inomata  _et al._ [2022] K. Inomata, E. McDonough, and W. Hu, Amplification of primordial perturbations from the rise or fall of the inflaton, JCAP 02 (02), 031, arXiv:2110.14641 [astro-ph.CO] .
  * Hooper  _et al._ [2023] D. Hooper, A. Ireland, G. Krnjaic, and A. Stebbins, Supermassive Primordial Black Holes From Inflation, (2023), arXiv:2308.00756 [astro-ph.CO] .
  * Erickcek and Sigurdson [2011] A. L. Erickcek and K. Sigurdson, Reheating Effects in the Matter Power Spectrum and Implications for Substructure, Phys. Rev. D 84, 083503 (2011), arXiv:1106.0536 [astro-ph.CO] .
  * Fan  _et al._ [2014] J. Fan, O. Özsoy, and S. Watson, Nonthermal histories and implications for structure formation, Phys. Rev. D 90, 043536 (2014), arXiv:1405.7373 [hep-ph] .
  * Erickcek [2015] A. L. Erickcek, The Dark Matter Annihilation Boost from Low-Temperature Reheating, Phys. Rev. D 92, 103505 (2015), arXiv:1504.03335 [astro-ph.CO] .
  * Dror  _et al._ [2018] J. A. Dror, E. Kuflik, B. Melcher, and S. Watson, Concentrated dark matter: Enhanced small-scale structure from codecaying dark matter, Phys. Rev. D 97, 063524 (2018), arXiv:1711.04773 [hep-ph] .
  * Kolb and Tkachev [1993] E. W. Kolb and I. I. Tkachev, Axion miniclusters and Bose stars, Phys. Rev. Lett. 71, 3051 (1993), arXiv:hep-ph/9303313 .
  * Visinelli  _et al._ [2018] L. Visinelli, S. Baum, J. Redondo, K. Freese, and F. Wilczek, Dilute and dense axion stars, Phys. Lett. B 777, 64 (2018), arXiv:1710.08910 [astro-ph.CO] .
  * Eggemeier and Niemeyer [2019] B. Eggemeier and J. C. Niemeyer, Formation and mass growth of axion stars in axion miniclusters, Phys. Rev. D 100, 063528 (2019), arXiv:1906.01348 [astro-ph.CO] .
  * Arvanitaki  _et al._ [2020] A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, Large-misalignment mechanism for the formation of compact axion structures: Signatures from the QCD axion to fuzzy dark matter, Phys. Rev. D 101, 083014 (2020), arXiv:1909.11665 [astro-ph.CO] .
  * Graham  _et al._ [2016] P. W. Graham, J. Mardon, and S. Rajendran, Vector Dark Matter from Inflationary Fluctuations, Phys. Rev. D 93, 103520 (2016), arXiv:1504.02102 [hep-ph] .
  * Gorghetto  _et al._ [2022] M. Gorghetto, E. Hardy, J. March-Russell, N. Song, and S. M. West, Dark photon stars: formation and role as dark matter substructure, JCAP 08 (08), 018, arXiv:2203.10100 [hep-ph] .
  * Gilman  _et al._ [2021] D. Gilman, J. Bovy, T. Treu, A. Nierenberg, S. Birrer, A. Benson, and O. Sameie, Strong lensing signatures of self-interacting dark matter in low-mass haloes, Mon. Not. Roy. Astron. Soc. 507, 2432 (2021), arXiv:2105.05259 [astro-ph.CO] .
  * Slone  _et al._ [2023] O. Slone, F. Jiang, M. Lisanti, and M. Kaplinghat, Orbital evolution of satellite galaxies in self-interacting dark matter models, Phys. Rev. D 107, 043014 (2023), arXiv:2108.03243 [astro-ph.CO] .
  * Zeng  _et al._ [2023] Z. C. Zeng, A. H. G. Peter, X. Du, S. Yang, A. Benson, F.-Y. Cyr-Racine, F. Jiang, C. Mace, and R. B. Metcalf, Till the core collapses: the evolution and properties of self-interacting dark matter subhalos, (2023), arXiv:2310.09910 [astro-ph.GA] .
  * Gad-Nasr  _et al._ [2023] S. Gad-Nasr, K. K. Boddy, M. Kaplinghat, N. J. Outmezguine, and L. Sagunski, On the Late-Time Evolution of Velocity-Dependent Self-Interacting Dark Matter Halos, (2023), arXiv:2312.09296 [astro-ph.GA] .
  * Mace  _et al._ [2024] C. Mace, Z. C. Zeng, A. H. G. Peter, X. Du, S. Yang, A. Benson, and M. Vogelsberger, Convergence Tests of Self-Interacting Dark Matter Simulations, (2024), arXiv:2402.01604 [astro-ph.GA] .
  * Gemmell  _et al._ [2023] C. Gemmell, S. Roy, X. Shen, D. Curtin, M. Lisanti, N. Murray, and P. F. Hopkins, Dissipative Dark Substructure: The Consequences of Atomic Dark Matter on Milky Way Analog Subhalos, (2023), arXiv:2311.02148 [astro-ph.GA] .
  * Roy  _et al._ [2023] S. Roy, X. Shen, M. Lisanti, D. Curtin, N. Murray, and P. F. Hopkins, Simulating Atomic Dark Matter in Milky Way Analogs, Astrophys. J. Lett. 954, L40 (2023), arXiv:2304.09878 [astro-ph.GA] .
  * Chang  _et al._ [2019] J. H. Chang, D. Egana-Ugrinovic, R. Essig, and C. Kouvaris, Structure Formation and Exotic Compact Objects in a Dissipative Dark Sector, JCAP 03, 036, arXiv:1812.07000 [hep-ph] .
  * Curtin and Setford [2020] D. Curtin and J. Setford, Signatures of Mirror Stars, JHEP 03, 041, arXiv:1909.04072 [hep-ph] .
  * Dvali  _et al._ [2020] G. Dvali, E. Koutsangelas, and F. Kuhnel, Compact Dark Matter Objects via $N$ Dark Sectors, Phys. Rev. D 101, 083533 (2020), arXiv:1911.13281 [astro-ph.CO] .
  * Zumalacarregui and Seljak [2018] M. Zumalacarregui and U. Seljak, Limits on stellar-mass compact objects as dark matter from gravitational lensing of type Ia supernovae, Phys. Rev. Lett. 121, 141101 (2018), arXiv:1712.02240 [astro-ph.CO] .
  * Dhawan and Mörstell [2023] S. Dhawan and E. Mörstell, Type ia supernova constraints on compact object dark matter, Monthly Notices of the Royal Astronomical Society 524, 5762 (2023).
  * DeRocco  _et al._ [2023] W. DeRocco, N. Smyth, and V. Takhistov, New Light on Dark Extended Lenses with the Roman Space Telescope, (2023), arXiv:2312.14782 [astro-ph.CO] .
  * Van Tilburg  _et al._ [2018] K. Van Tilburg, A.-M. Taki, and N. Weiner, Halometry from Astrometry, JCAP 07, 041, arXiv:1804.01991 [astro-ph.CO] .
  * Banik  _et al._ [2021] N. Banik, J. Bovy, G. Bertone, D. Erkal, and T. J. L. de Boer, Evidence of a population of dark subhaloes from $Gaia$ and Pan-STARRS observations of the GD-1 stream, Mon. Not. Roy. Astron. Soc. 502, 2364 (2021), arXiv:1911.02662 [astro-ph.GA] .
  * Ando  _et al._ [2022] S. Ando, N. Hiroshima, and K. Ishiwata, Constraining the primordial curvature perturbation using dark matter substructure, Phys. Rev. D 106, 103014 (2022), arXiv:2207.05747 [astro-ph.CO] .
  * Graham and Ramani [2023] P. W. Graham and H. Ramani, Constraints on Dark Matter from Dynamical Heating of Stars in Ultrafaint Dwarfs. Part 1: MACHOs and Primordial Black Holes, (2023), arXiv:2311.07654 [hep-ph] .
  * Koulen  _et al._ [2024] J. M. Koulen, S. Profumo, and N. Smyth, Constraints on Primordial Black Holes from $N$-body simulations of the Eridanus II Stellar Cluster, (2024), arXiv:2403.19015 [astro-ph.CO] .
  * Brandt [2016] T. D. Brandt, Constraints on MACHO Dark Matter from Compact Stellar Systems in Ultra-Faint Dwarf Galaxies, Astrophys. J. Lett. 824, L31 (2016), arXiv:1605.03665 [astro-ph.GA] .
  * Lu  _et al._ [2021] P. Lu, V. Takhistov, G. B. Gelmini, K. Hayashi, Y. Inoue, and A. Kusenko, Constraining Primordial Black Holes with Dwarf Galaxy Heating, Astrophys. J. Lett. 908, L23 (2021), arXiv:2007.02213 [astro-ph.CO] .
  * Koushiappas and Loeb [2017] S. M. Koushiappas and A. Loeb, Dynamics of Dwarf Galaxies Disfavor Stellar-Mass Black Holes as Dark Matter, Phys. Rev. Lett. 119, 041102 (2017), arXiv:1704.01668 [astro-ph.GA] .
  * Zhu  _et al._ [2018] Q. Zhu, E. Vasiliev, Y. Li, and Y. Jing, Primordial black holes as dark matter: constraints from compact ultra-faint dwarfs, Monthly Notices of the Royal Astronomical Society 476, 2 (2018).
  * Stegmann  _et al._ [2020] J. Stegmann, P. R. Capelo, E. Bortolas, and L. Mayer, Improved constraints from ultra-faint dwarf galaxies on primordial black holes as dark matter, Mon. Not. Roy. Astron. Soc. 492, 5247 (2020), arXiv:1910.04793 [astro-ph.GA] .
  * Zoutendijk  _et al._ [2020] S. L. Zoutendijk, J. Brinchmann, L. A. Boogaard, M. L. Gunawardhana, T.-O. Husser, S. Kamann, A. F. R. Padilla, M. M. Roth, R. Bacon, M. Den Brok, _et al._ , The muse-faint survey-i. spectroscopic evidence for a star cluster in eridanus 2 and constraints on machos as a constituent of dark matter, Astronomy & Astrophysics 635, A107 (2020).
  * Wadekar and Wang [2023] D. Wadekar and Z. Wang, Constraining axion and compact dark matter with interstellar medium heating, Phys. Rev. D 107, 083011 (2023), arXiv:2211.07668 [hep-ph] .
  * van den Bosch  _et al._ [2018] F. C. van den Bosch, G. Ogiya, O. Hahn, and A. Burkert, Disruption of Dark Matter Substructure: Fact or Fiction?, Mon. Not. Roy. Astron. Soc. 474, 3043 (2018), arXiv:1711.05276 [astro-ph.GA] .
  * Chluba  _et al._ [2012] J. Chluba, A. L. Erickcek, and I. Ben-Dayan, Probing the inflaton: Small-scale power spectrum constraints from measurements of the CMB energy spectrum, Astrophys. J. 758, 76 (2012), arXiv:1203.2681 [astro-ph.CO] .
  * Akrami  _et al._ [2020] Y. Akrami, F. Arroja, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. Barreiro, N. Bartolo, S. Basak, _et al._ , Planck 2018 results-x. constraints on inflation, Astronomy & Astrophysics 641, A10 (2020).
  * Bird  _et al._ [2011] S. Bird, H. V. Peiris, M. Viel, and L. Verde, Minimally parametric power spectrum reconstruction from the lyman $\alpha$ forest, Monthly Notices of the Royal Astronomical Society 413, 1717 (2011).
  * Press and Schechter [1974] W. H. Press and P. Schechter, Formation of galaxies and clusters of galaxies by self-similar gravitational condensation, Astrophysical Journal, Vol. 187, pp. 425-438 (1974) 187, 425 (1974).
  * Sheth and Tormen [1999] R. K. Sheth and G. Tormen, Large-scale bias and the peak background split, Monthly Notices of the Royal Astronomical Society 308, 119 (1999).
  * Watson  _et al._ [2013] W. A. Watson, I. T. Iliev, A. D’Aloisio, A. Knebe, P. R. Shapiro, and G. Yepes, The halo mass function through the cosmic ages, Monthly Notices of the Royal Astronomical Society 433, 1230 (2013).
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Gilman  _et al._ [2022] D. Gilman, A. Benson, J. Bovy, S. Birrer, T. Treu, and A. Nierenberg, The primordial matter power spectrum on sub-galactic scales, Monthly Notices of the Royal Astronomical Society 512, 3163 (2022).
  * Eisenstein and Hu [1999] D. J. Eisenstein and W. Hu, Power spectra for cold dark matter and its variants, The Astrophysical Journal 511, 5 (1999).
  * Diemer [2018] B. Diemer, Colossus: A python toolkit for cosmology, large-scale structure, and dark matter halos, The Astrophysical Journal Supplement Series 239, 35 (2018).
  * Redi and Tesi [2023] M. Redi and A. Tesi, Meso-inflationary QCD axion, Phys. Rev. D 107, 095032 (2023), arXiv:2211.06421 [hep-ph] .
  * Chluba and Grin [2013] J. Chluba and D. Grin, CMB spectral distortions from small-scale isocurvature fluctuations, Mon. Not. Roy. Astron. Soc. 434, 1619 (2013), arXiv:1304.4596 [astro-ph.CO] .
  * Amin and Mirbabayi [2022] M. A. Amin and M. Mirbabayi, A lower bound on dark matter mass, (2022), arXiv:2211.09775 [hep-ph] .
  * Byrnes  _et al._ [2019] C. T. Byrnes, P. S. Cole, and S. P. Patil, Steepest growth of the power spectrum and primordial black holes, JCAP 06, 028, arXiv:1811.11158 [astro-ph.CO] .
  * Croon and Sevillano Muñoz [2024] D. Croon and S. Sevillano Muñoz, Cosmic microwave background constraints on extended dark matter objects, (2024), arXiv:2403.13072 [astro-ph.CO] .