# Nonrelativistic formation of scalar clumps as a candidate for dark matter

Philippe Brax, Patrick Valageas, Jose A.R. Cembranos

► To cite this version:

Philippe Brax, Patrick Valageas, Jose A.R. Cembranos. Nonrelativistic formation of scalar clumps as a candidate for dark matter. Physical Review D, 2020, 102 (8), pp.083012. 10.1103/Phys-RevD.102.083012. hal-02908368

HAL Id: hal-02908368

https://hal.science/hal-02908368

Submitted on 23 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Nonrelativistic formation of scalar clumps as a candidate for dark matter

Philippe Brax and Patrick Valageas

Institut de Physique Théorique, Université Paris-Saclay,
CEA, CNRS, F-91191 Gif-sur-Yvette Cedex, France

Jose A. R. Cembranos

Departamento de Física Teórica and IPARCOS, Universidad Complutense de Madrid,
E-28040 Madrid, Spain

<figure><img src="image_2.png" /></figure>

(Received 29 July 2020; accepted 23 September 2020; published 13 October 2020)

We propose a new mechanism for the formation of dark matter clumps in the radiation era. We assume that a light scalar field is decoupled from matter and oscillates harmonically around its vacuum expectation value. We include self-interactions and consider the nonrelativistic regime. The scalar dynamics are described by a fluid approach where the fluid pressure depends on both quantum and self-interaction effects. When the squared speed of sound of the scalar fluid becomes negative, an instability arises and the fluctuations of the scalar energy-density field start growing. They eventually become nonlinear and clumps form. Subsequently, the clumps aggregate and reach a universal regime. Afterwards, they play the role of cold dark matter. We apply this mechanism first to a model with a negative quartic term stabilized by a positive self-interaction of order six, and then to axion monodromy, where a subdominant cosine potential corrects a mass term. In the first case, the squared speed of sound becomes negative when the quartic term dominates, leading to a tachyonic instability. For axion monodromy, the instability starts very slowly after the squared speed of sound first becomes negative and then oscillates around zero. Initially the density perturbations perform acoustic oscillations due to the quantum pressure. Eventually, they start growing exponentially due to a parametric resonance. The shape and the scaling laws of the clumps depend on their formation mechanism. When the tachyonic phase takes place, the core density of the clumps is uniquely determined by the energy density at the beginning of the instability. On the other hand, for axion monodromy, the core density scales with the soliton mass and radius. This difference comes from the crucial role that the quantum pressure plays in both the parametric resonance in the linear regime and in the nonlinear formation regime of static scalar solitons. In both scenarios, the scalar-field clumps span a wide range of scales and masses, running from the size of atoms to that of galactic molecular clouds, and from $10^{-3}$ gram to thousands of solar masses. Because of finite-size effects, both from the source and the lens, these dark matter clumps are far beyond the reach of microlensing observations. We find that the formation redshift of the scalar clumps can span a large range in the radiation era; the associated background temperature can vary from $10$ eV to $10^5$ GeV, and the scalar-field mass from $10^{-26}$ GeV to $10$ GeV.

DOI: 10.1103/PhysRevD.102.083012

## I. INTRODUCTION

Dark matter is an essential ingredient of both astrophysics and cosmology. On very large scales, its description as a cold and pressureless fluid fits all our observations and is part of the standard model of cosmology [1]. From a more fundamental point of view, particle physicists have tried in the last few decades to find natural models of dark matter involving weakly interacting massive particles (WIMPs). Despite very promising theoretical candidates, such as neutralinos in supersymmetric models [2], no experimental evidence of WIMPs has emerged from data so far. This has triggered a recent revival of alternative models, where the origin of dark matter does not follow from the freezing out of particle interactions below a certain temperature. Indeed,

ever since the analysis of the strong CP problem and the design of the axion mechanism, the possibility that dark matter could result from the coherent oscillations of a scalar field around the minimum of its potential has been considered [3–5]. Since then, axions [6–9] and axionlike particles (ALPs) [10] have also been under intense scrutiny, both from the theoretical and the experimental sides. However, they are just one particular example of more general scalar dark-matter models, which can be described from an effective-field-theory point of view as parametrized by their masses and the set of their self-interactions [11]. In the simplest case, called fuzzy dark matter [12,13], a massive scalar field oscillating around its vacuum expectation value (vev), and with a sufficiently low mass
$m \lesssim 10^{-21}$ eV, could play the role of dark matter. The resulting properties of these scalar dark-matter models are similar to the standard cold dark matter (CDM) for the formation of large-scale structures [14–21], but not for small scales, where distinctive features such as a non-vanishing speed of sound can leave different observational signatures [19,20,22–35].

In this paper, we focus on one scalar field of mass $m$ whose self-interactions are subdominant compared to the quadratic mass term. We also consider that the oscillating scalar field is mostly time-dependent, with small space-dependent effects. In momentum space, this amounts to considering that momenta are small compared to the mass and that therefore the nonrelativistic approximation of the dynamics is valid. In this case, and after reducing the field to its equivalent quantum-mechanical picture, the Madelung transformation allows one to describe the evolution of the scalar field in terms of a fluid with nontrivial pressure terms. The first one, which is present even in the absence of self-interaction, is the so-called quantum pressure. It has a repulsive effect and allows for the formation of large solitons where the quantum pressure can balance the gravitational attraction [13]. These quantum-pressure effects are at the heart of the fuzzy dark matter models, and require a mass term of low value, $m \lesssim 10^{-21}$ eV. When self-interactions are present and overcome the quantum pressure, which can happen for masses larger than $10^{-21}$ eV, it is known that a repulsive pressure can also result from positive interaction potentials like $\phi^4$. In this case, large solitons can also form by gravitational instability and eventually stabilize when the gravitational attraction becomes balanced by the scalar self-repulsion [33–36]. In all these cases, dark matter is essentially a smooth fluid with large overdensities in the core of solitonic objects, which have galactic sizes and could play the role of galactic dark matter halos with a smooth inner region. This may alleviate some of the galactic-scale tensions with the data encountered by the standard CDM scenario.

On the other hand, and inspired by the possible representation of dark matter in the form of primordial black holes [37] or massive compact halo objects (MACHOs) [38], it can be envisaged that scalar-field clumps of much smaller sizes could exist and form all (or a large portion) of the dark matter content. In this case, the sign of the quartic self-interactions is crucial [39,40]. For axions, this sign is negative leading to an attractive interaction between particles in the nonrelativistic regime. In the relativistic regime, the self-interaction can become dominant and lead to the formation of “axitons” as the squared mass of the axion can become negative for large excursions of the field, leading to a potential instability and the formation of clumps [41,42]. Another type of phenomenon, which leads to the creation of “oscillons” in some scalar field theories [43–46], has been attributed to an instability where the

modes can undergo a parametric resonance [45] and the creation of nonlinear overdensities. In this paper, we will deal with similar mechanisms in the nonrelativistic regime. In this case, a tachyonic instability can be triggered when the speed of sound squared coming from the self-interactions becomes negative. This happens for simple models with a negative quartic interaction compensated by a positive one of degree six. This instability is counterbalanced by the quantum pressure on small scales and by the repulsive effects of the order-six term on larger scales, leading to the creation of clumps mostly determined by the higher-order operator contribution. In another scenario, the model being of the axion monodromy type, with a scalar potential where a large mass term is modulated by small oscillations, the instability felt by the perturbations of the scalar energy density is due to a parametric resonance triggered when the speed of sound squared becomes negative too. In this case, the ensuing parametric-resonance growth of the instability is delayed by a long period of acoustic oscillations sustained by the quantum pressure. In most of these cases, gravity does not play a role and the clumps are formed in the radiation era. Their subsequent evolution first involves their 2-body collisional aggregation and relaxation toward stable halos, which are next diluted by the expansion of the Universe. At lower redshifts, these small scalar-field solitons play the role of dark matter particles, in the same manner as primordial black holes or MACHOs, and we recover the standard CDM cosmology. Although these scalar clumps are usually much smaller than galactic cores that can form in the fuzzy dark matter models, e.g., they can be as small as one angstrom, they can also reach sub-galactic sizes of the parsec scale, like galactic molecular clouds. Hence these scenarios lead to a wide range of possible dark matter scales.

The two types of formation mechanisms that we consider lead to very different properties for the clumps. In the tachyonic case, with a polynomial potential, the density in the core of the clumps is determined by the features of the potential, i.e., the energy density where the self-interactions change from being attractive to repulsive (which also sets the background energy density at the beginning of the instability). For axion monodromy this is not the case, as the clumps can accommodate a continuous distribution of energy densities in their core. This sharp difference follows from the nature of the energy functional of the clumps as a function of the energy density. In the tachyonic case, the potential energy of the clumps admits a minimum which characterizes the density of the clumps, giving a mass-radius relation $M \sim R^3$. In the axion monodromy setup, the potential energy is a decreasing function which does not select a unique equilibrium density, resulting in a $M \sim R^5$ mass-radius relation when the self-interaction dominates, and $M \sim R^4$ when gravity becomes the relevant interaction after the nonlinear collapse of the structures triggered by the parametric resonance instability.
As already stated, the dynamics comprise two steps. The first one, which we have just described, results from the type of instability of the fundamental model describing the physics of the scalar field, e.g., a polynomial interaction potential vs axion monodromy. The second stage happens postformation and follows a short aggregation phase, which can influence the final mass and radius of the clumps. We describe in detail how this aggregation process depends on the mass-radius relationship of the clumps and therefore on the initial formation mechanism. Whereas in the polynomial interaction case the aggregation process leads to a significant growth of the size and mass of the clumps, in the axion monodromy case the aggregation is not very efficient and the mass and radii are unaffected.

The results that we present in this paper use two main ingredients. The first one is the leading-order fast harmonic motion of the field, with a frequency set by its mass $m$, and the second one is the existence of an instability in the growth of the energy density contrast, which is triggered by the negative sign of the speed of sound squared.

The leading-order harmonic motion is guaranteed by the smallness of the perturbations to the scalar potential compared to the leading quadratic term. This corresponds to models with typically two scales associated with two contributions of different origins to the scalar-field potential. The first contribution, with a large amplitude, is given by a quadratic term and gives rise to the leading-order fast harmonic motion. The second contribution, with a small amplitude, is such that its nonlinear orders cannot be neglected. We will consider two cases, (a) when the small-amplitude self-interaction corrections to the quadratic term are slow varying functions such as a low-order polynomial, and (b) when they show fast oscillations, such as a cosine term. The leading-order harmonic oscillations due to the quadratic term in the scalar potential ensure that the scalar field behaves like dark matter (with a mean density decaying as $1/a^3$ with the expansion of the Universe). The subleading self-interactions however play a critical role, as they can lead to instabilities and the fragmentation of the homogeneous dark-matter distribution.

The two types of instabilities that we exemplify, i.e., the tachyonic and parametric resonance, have been considered in the literature in several contexts. In a recent paper, the case of the “large-misalignment mechanism” [40] was presented. In this scenario, and taking the cosine axion potential as an example, if the field starts initially close enough to the top of the potential, the instability due to the negative quartic term of the cosine function near the origin is delayed and a parametric resonance instability sets in. This leads to the formation of clumps which can be described as “solitons” when the gravitational attraction is balanced by the kinetic pressure and “oscillons” when gravity is irrelevant. Their (meta)-stability is entirely due to the scalar self-interactions. This scenario applies to the QCD axion and certain axion monodromy potentials which

are flatter than quadratic for large field values. In our analysis of the axion monodromy models, with large quadratic potentials perturbed by a small cosine interaction, we preserve the harmonic motion at the leading order throughout our description of the parametric resonance instability. In this dominant-quadratic-term scenario, the nonrelativistic approximation applies throughout. In this case, the parametric resonance instability appears well before the argument of the cosine potential becomes small. Moreover, the speed of sound squared becomes negative well before the parametric resonance starts too. Contrary to the “large-misalignment mechanism,” where the delay in the growth of perturbations is due to the flatness of the interaction potential initially, in our case the delay is due to the effects of the quantum pressure, which drives initial acoustic oscillations before becoming low enough and allowing the onset of the parametric resonance. This delay can also be understood as the time it takes for these acoustic oscillations to become tuned to the frequency set by the cosine self-interaction potential (thanks to their time dependence, due to the expansion of the Universe and the decrease of the background density), so that a resonance can develop. Nonrelativistic clumps and their formation have been analyzed numerically in a recent paper [47], where a potential with a negative $\phi^4$ interaction term close to the origin was completed by higher order terms, eventually leading to a bounded potential for large field values. In this setting, the tachyonic instability plays a prominent role in the formation of the nonrelativistic clumps. Eventually, nonlinear effects take over and individual clumps form with little scalar interactions between each other. Later in the evolution, this gas of clumps is affected by the gravitational attraction and they start moving toward each other. In this paper, we also present a similar mechanism for the formation of clumps through a tachyonic instability and their stabilization by higher order terms in the scalar potential. Then, we analyze the early aggregation process before the dilution by the expansion of the Universe. We also describe the same process in a thermodynamic way. We pay particular attention to the parameter space combining theoretical self-consistency conditions and standard requirements (the formation of the dark matter clumps should occur before matter-radiation equality and their size should not exceed the parsec scale). We also check that their gravitational potential well is too weak to form black holes. Finally, finite-size effects imply that they cannot be detected by microlensing observations. We find that the scalar field can have a mass $m$ ranging from $10^{-17}$ eV to 10 GeV, giving rise to dark matter clumps that range from the size of atoms to that of galactic molecular clouds.

The paper is arranged as follows. In Sec. II, we review the classical field model associated with such a scalar field, and its nonrelativistic regime. In Sec. III, we describe our first scenario, associated with the tachyonic instability
where the speed of sound squared becomes negative at low background densities. We first use a perturbative approach in Sec. III B, to follow the growth of the scalar-field density perturbations. In Sec. III C, we study the stable isolated scalar-field configurations that arise in such a model, i.e., the “solitons” that correspond to the final dark matter clumps. We estimate in Sec. III D the efficiency of the collisional aggregation of these scalar clouds, shortly after their formation and before they are diluted by the expansion of the Universe, and we check in Sec. III E that they do not collapse to black holes. Then, in Sec. III F, we take into account theoretical constraints to compute the parameter space of this scenario. In Sec. III G we compute the scales spanned by the scalar dark-matter clumps and in Sec. III H we check that they are far beyond the reach of microlensing observations.

Next, in Sec. IV, we present a different mechanism for clump formation, associated with a parametric resonance. We take as an example a Lagrangian inspired from axion monodromy, where a dominant mass term is corrected by a subleading cosine term. The parametric resonance then arises from the interplay between this oscillating self-interaction term, the quantum pressure, and the kinetic terms of the scalar field. We again describe the perturbative growth of the scalar-field density fluctuations and the stable solitons that can arise. We also compute the parameter space of this second scenario and the size of the scalar clumps. Again, we check that they do not collapse into black holes and are much below the observational threshold of microlensing observations.

We present our main conclusions in Sec. V. We finally complete our discussion with different Appendices on thermodynamical phase transitions, parametric resonance, and soliton profiles.

## II. CLASSICAL FIELDS AND THEIR NONRELATIVISTIC LIMIT

### A. Classicality

In the following, we shall be interested in models of scalar dark matter where the dark-matter field can be described classically. This is a reasonable approximation for the quantum field $\phi$, whose nonrelativistic behavior will give rise to dark matter, if the occupation number $N$ of the associated quantum state is very large. Denoting by $\rho$ the energy density of the field and by $n = \rho/m$ the number density, where $m$ is the mass of the scalar, the occupation number can be estimated as [48]

$$N \approx \frac{\rho}{m} \lambda_{\text{dB}}^3, \qquad \lambda_{\text{dB}} = \frac{2\pi}{mv}, \qquad (1)$$

where $\lambda_{\text{dB}}$ is the de Broglie wavelength of the scalar particles associated to $\phi$. Here $v$ is their typical velocity. This gives the condition for classicality

$$N \sim \frac{\rho}{m^4 v^3} \gg 1. \qquad (2)$$

We can envisage two types of situations. In the first one, the energy density of the scalar field is nearly homogeneously distributed in the Universe and behaves like $\rho \approx \rho_0/a^3$, where $\rho_0$ is the present dark-matter density in the Universe. Inside large-scale inhomogeneities such as galaxy halos, the typical velocity of dark-matter particles $v_0$ is small and the classical regime is attained when

$$m^4 v_0^3 \ll \rho_0 \sim 10^{-48} \text{ GeV}^4, \qquad (3)$$

where we consider low redshifts in the matter era. As we expect $v_0 \approx 10^{-3}$, this is the case when

$$\text{cosmological inhomogeneities only: } m \ll 0.1 \text{ eV.} \qquad (4)$$

In this mass range the field can be treated classically. This also applies at higher redshifts, as $\rho \propto a^{-3}$ and typically $v \sim a^{-1}$ because of the expansion of the Universe.

Another scenario is the one that we consider in this paper: dark matter is made of scalar-field clumps created in the radiation era and forming a bound state of dark-matter fluid. Then, in a fashion similar to primordial black holes, these clumps play the role of dark matter particles and behave at late times as in standard CDM cosmologies. In this case, the density $\rho$ is large inside the clumps, reflecting the large energy densities at the time of their formation, and the velocity is negligible as these clumps are equilibrium configurations. Hence, for such clumps $N$ will be very large and we can treat $\phi$ as a classical field. In fact, the classicality condition (2) will provide a self-consistency constraint on the parameter space of the scenarios we study in this paper.

### B. Equations of motion

We focus on scalar-field models characterized by canonical kinetic terms and an interaction potential $V_1(\phi)$. Thus, they are governed by the action

$$S[\phi] = \int d^4x \sqrt{-g} \left[ -\frac{1}{2} g^{\mu\nu} \partial_{\mu} \phi \partial_{\nu} \phi - V(\phi) \right], \qquad (5)$$

with

$$V(\phi) = \frac{1}{2} m^2 \phi^2 + V_1(\phi). \qquad (6)$$

In this paper, we restrict our study to the nonrelativistic regime, when the self-interactions are small as compared with the quadratic part,

$$V_1 \ll \frac{1}{2} m^2 \phi^2. \qquad (7)$$
At linear order in the gravitational potential $\Phi$ and for $m \gg H$, where $H$ is the Hubble expansion rate, the equation of motion of the real scalar field $\phi$ in a perturbed Friedmann-Lemaître-Robertson-Walker universe (FLRW) is

$$ \ddot{\phi} + 3H\dot{\phi} - \frac{1}{a^2}\nabla^2\phi + (1+2\Phi)m^2\phi + \frac{dV_1}{d\phi} = 0, \quad (8) $$

where $a$ is the scale factor of the Universe, normalized to unity now. As we are interested in the classical behavior of the field $\phi$ in the nonrelativistic limit, it is convenient to decompose

$$ \phi = \frac{1}{\sqrt{2m}}(\psi e^{-imt} + \psi^* e^{imt}), \quad (9) $$

when the spatial and time variations of $\psi$ are small compared to $m$. This ansatz emphasizes the fact that the scalar field oscillates with a pulsation $m$ as the quadratic terms in the scalar field action (5) dominate, following (7). From this we can deduce the equation of motion of the nonrelativistic complex scalar field $\psi$,

$$ i\left(\dot{\psi} + \frac{3}{2}H\psi\right) = -\frac{\nabla^2\psi}{2ma^2} + m\Phi\psi + \frac{\partial V_1}{\partial \psi^*}, \quad (10) $$

which is a nonlinear version of the Schrödinger equation. Here we introduced the effective nonrelativistic self-interaction potential $V_1(\psi, \psi^*)$, which is obtained from $V_1$ by averaging over the leading oscillations $e^{\pm imt}$ of $\phi$. For polynomial self-interactions, or analytic potentials that can be defined by their Taylor expansion, with

$$ V_1(\phi) = \Lambda^4 \sum_{p \ge 3} \frac{\lambda_p}{p} \left(\frac{\phi}{\Lambda}\right)^p, \quad (11) $$

one obtains [33]

$$ V_1(\psi, \psi^*) = \Lambda^4 \sum_{p \ge 2} \frac{\lambda_{2p}}{2p} \frac{(2p)!}{(p!)^2} \left( \frac{\psi \psi^*}{2m\Lambda^2} \right)^p. \quad (12) $$

It is convenient to introduce the Madelung transform [49]

$$ \psi = \sqrt{\frac{\rho}{m}} e^{iS}. \quad (13) $$

This defines the effective density field $\rho$, which coincides with the scalar-field energy density in this nonrelativistic limit. The phase $S$ defines an effective curl-free velocity field $\vec{v}$,

$$ \vec{v} = \frac{\vec{\nabla} S}{ma}. \quad (14) $$

Then, the equations of motion take a familiar form, i.e., the one of hydrodynamics [39]. The real part of the nonlinear Schrödinger equation gives the continuity equation

$$ \dot{\rho} + 3H\rho + \frac{1}{a}\nabla \cdot (\rho\vec{v}) = 0. \quad (15) $$

We can see that the self-interactions due to $V_1$ do not modify this continuity equation. The imaginary part of the nonlinear Schrödinger equation becomes the Hamilton-Jacobi relation

$$ \dot{S} + \frac{(\nabla S)^2}{2ma^2} = -m\Phi - m\frac{dV_1}{d\rho} + \frac{1}{2ma^2}\frac{\nabla^2\sqrt{\rho}}{\sqrt{\rho}}, \quad (16) $$

where the nonrelativistic self-interaction potential $V_1(\rho)$ is directly obtained from $V_1(\psi, \psi^*)$ in Eq. (12) with $\psi\psi^* = \rho/m$,

$$ V_1(\rho) = \Lambda^4 \sum_{p \ge 2} \frac{\lambda_{2p}}{2p} \frac{(2p)!}{(p!)^2} \left( \frac{\rho}{2m^2\Lambda^2} \right)^p. \quad (17) $$

Then, taking the gradient of Eq. (16) gives the hydrodynamical Euler equation,

$$ \dot{\vec{v}} + H\vec{v} + \frac{1}{a}(\vec{v} \cdot \nabla)\vec{v} = -\frac{1}{a}\nabla(\Phi + \Phi_1 + \Phi_Q), \quad (18) $$

where we used $\nabla(\vec{v}^2) = 2(\vec{v} \cdot \nabla)\vec{v}$ as $\nabla \times \vec{v} = 0$. The self-interaction potential $\Phi_1(\rho)$ is defined by

$$ \Phi_1(\rho) = \frac{dV_1}{d\rho}, \quad (19) $$

and we have introduced the “quantum pressure” term

$$ \Phi_Q = - \frac{\nabla^2 \sqrt{\rho}}{2m^2 a^2 \sqrt{\rho}}. \quad (20) $$

The continuity equation and the Euler equation will show unstable solutions in the examples we consider in this article, because of attractive self-interactions $\Phi_1$ at low densities. This description is valid provided the nonlinear terms are small compared to the quadratic terms in the original action, as in (7). This translates into the conditions

$$ V_1 \ll \rho, \quad \text{hence} \quad \Phi_1 \ll 1. \quad (21) $$

### C. Cosmological background

#### 1. Real scalar field $\phi$

We now restrict our attention to the cosmological background, where the scalar field $\vec{\phi}$ only depends on time. The corresponding equation of motion is
$$ \ddot{\bar{\phi}} + 3H\dot{\bar{\phi}} + m^2\bar{\phi} + \frac{dV_1}{d\bar{\phi}} = 0, \quad (22) $$

whose solution can be written as a slowly varying deformation of the harmonic oscillator,

$$ \bar{\phi}(t) = \bar{\phi}(t) \cos(mt - \bar{S}(t)). \quad (23) $$

Notice the similarity with the ansatz (9) defining the complex scalar field $\psi$. The amplitude of the scalar field evolves in time and decreases with the scale factor

$$ \bar{\phi} = \bar{\phi}_0 a^{-3/2}, \quad (24) $$

while the phase evolves according to

$$ \bar{S}(t) = \bar{S}_0 - \int_{t_0}^{t} dt m \Phi_1 \left( \frac{m^2 \bar{\phi}_0^2}{2a^3} \right). \quad (25) $$

Hence, at the background level, the scalar field oscillates harmonically at the leading order, with the high frequency $m$ given by the scalar mass. The Hubble expansion and the self-interactions give rise to a slow decay of the amplitude and to a phase shift. The power-law decay $\bar{\phi} \propto a^{-3/2}$ shows that the scalar-field energy density $\bar{\rho}_{\phi} \simeq m^2 \dot{\bar{\phi}}^2 / 2$ decreases like $a^{-3}$ and plays the role of a nonrelativistic dark-matter component.

## 2. Nonrelativistic limit

Comparing the solution (23) with the nonrelativistic decomposition (9), we can see that, at the background level, the complex scalar field $\bar{\psi}$ is

$$ \bar{\psi}(t) = \bar{\psi}_0 a^{-3/2} e^{i\bar{S}}, \quad \text{with} \quad \bar{\psi}_0 = \sqrt{\frac{m}{2}} \bar{\phi}_0 = \sqrt{\frac{\bar{\rho}_0}{m}}. \quad (26) $$

<figure><img src="image_3.png" /></figure>

We can check that the solution defined by $\bar{\rho} = \bar{\rho}_0 / a^3$ and $\bar{S}$ given by Eq. (25), which also can be written as

$$ \dot{\bar{S}} = - \frac{m \Lambda^4 a^3}{2 \bar{\rho}_0} \sum_{n=2}^{\infty} \lambda_{2n} \frac{(2n)!}{(n!)^2} \left( \frac{\bar{\rho}_0}{2m^2 \Lambda^2 a^3} \right)^n, \quad (27) $$

is indeed the solution of the equations of motion derived from the hydrodynamical action, which read

$$ \dot{\bar{S}} = -m \frac{d\bar{V}_1}{d\bar{\rho}}, \quad (28) $$

$$ \dot{\bar{\rho}} + 3H\bar{\rho} = 0. \quad (29) $$

Hence, at the background level, the evolution of the scalar field given by the hydrodynamical equations reproduces the full solution to the scalar-field equation (22).

# III. TACHYONIC INSTABILITY FOR SMOOTH SELF-INTERACTIONS

## A. Polynomial self-interactions

In the first part of this paper, we consider the scenario illustrated in Fig. 1, associated with slowly-varying self-interaction potentials. For template, we take a low-order polynomial case where we directly define the model at the nonrelativistic level,

$$ \Phi_I = -c_1 \frac{\rho}{\rho_\Lambda} + c_2 \frac{\rho^2}{\rho_\Lambda^2}, \quad \mathcal{V}_I = -c_1 \frac{\rho^2}{2\rho_\Lambda} + c_2 \frac{\rho^3}{3\rho_\Lambda^2}, \quad (30) $$

with $c_i > 0$. This corresponds to

$$ V_I(\phi) = - \frac{c_1 m^4}{3 \rho_\Lambda} \phi^4 + \frac{2 c_2 m^6}{15 \rho_\Lambda^2} \phi^6. \quad (31) $$

<figure><img src="image_4.png" /><figcaption>FIG. 1. The main stages of the formation of scalar dark-matter clumps for the tachyonic scenario (31). Cosmic time grows from the left column to the right column, and from the upper panel to the lower panel within each column. See the main text for explanations.</figcaption></figure>
We focus on the choice of parameters

$$c_1 \sim c_2 \ll 1, \quad (32)$$

meaning that $\rho_\Lambda$ is the characteristic density that governs the shape of the nonlinear self-interaction potential $V_1$. As we will consider redshifts where $\bar{\rho} \sim \rho_\Lambda$, when the change of shape of the potential can be felt, we also require $c_i \ll 1$ to fulfil the nonrelativistic condition (21). In practice, at densities $\rho \gtrsim \rho_\Lambda$, we can expect higher-order terms to come into play, if (30) is understood as a Taylor expansion in powers of $\rho$, originating from the Taylor expansion in $\phi$ of $V_1(\phi)$. However, the physics will not change, as long as $d\Phi_1/d\rho$ shows one change of sign, being negative at low densities and positive at high densities.

To facilitate the reading of this section, we already present in Fig. 1 the formation process of the scalar-field clumps that will play the role of dark-matter particles at low redshifts. From the first to the fourth column, this goes as follows.

(1) The scalar field $\phi$ quickly oscillates in the potential $V(\phi)$, which is dominated by its quadratic component with a small correction $V_1$. In the nonrelativistic regime, we can integrate over the fast oscillations of $\phi$. The slow dynamics is then described by the complex scalar field $\psi$, or the hydrodynamics density and velocity fields $\{\rho, \vec{v}\}$, and the self-interaction potential $\Phi_1(\rho)$ defined by Eq. (19), i.e., Eq. (30) in our polynomial example.

(2) At early times, when $d\Phi_1/d\rho > 0$, the scalar-field density perturbations oscillate as acoustic waves. As the background density $\bar{\rho}$ decreases with time, it finally enters the regime where $d\Phi_1/d\rho < 0$. This quickly leads to a tachyonic instability ($c_s^2 < 0$) for some intermediate wave numbers $k$ and an exponential growth of the density contrast $\delta(\vec{k})$.

(3) The scalar density field then quickly reaches the nonlinear regime and the overdense regions collapse to form stable configurations (solitons).

(4) Because of their non-negligible velocities, these scalar clouds collide and grow by aggregation, relaxing toward more massive solitons. Next, the expansion of the Universe dilutes these scalar clumps, which behave as isolated CDM clumps. At lower redshifts, the amplification by gravitational instability of perturbations on much larger scales will form the cosmic web, galaxies, and clusters, as in the standard $\Lambda$CDM scenario.

We describe in the following sections these various stages in more detail.

## B. Cosmological perturbations

### 1. Linear theory

For small perturbations with respect to the FLRW background, we can linearize the equations of motion.

As explained in the previous sections, and as illustrated by the first column in Fig. 1, in the nonrelativistic regime it is convenient to work with the fluid approach, where the fast harmonic oscillations of the scalar field $\phi$ have been integrated out and we are left with the density-dependent self-interaction potential $\Phi_1(\rho)$. Then, defining the linear density contrast $\delta$ and the divergence $\theta$ of the fluid velocity,

$$\delta = \frac{\rho - \bar{\rho}}{\bar{\rho}}, \quad \theta = \frac{\nabla \cdot \vec{v}}{a}, \quad (33)$$

the continuity equation gives the familiar constraint between the density contrast and the divergence of the velocity field,

$$\theta = -\dot{\delta}, \quad (34)$$

while the Euler equation (18) implies

$$\dot{\theta} + 2H\theta = -\frac{1}{a^2}\nabla^2(\Phi + \Phi_1 + \Phi_Q). \quad (35)$$

Combining these two equations, and upon using the Poisson equation $\nabla^2\Phi = 4\pi G a^2 \bar{\rho} \delta$ and the expression (20) of the quantum potential, we obtain, in Fourier space, the modified growth equation [36]

$$\ddot{\delta} + 2H\dot{\delta} + \left( c_s^2 \frac{k^2}{a^2} - 4\pi G \bar{\rho} \right) \delta = 0, \quad (36)$$

where we introduced the speed of sound $c_s$ as

$$c_s^2 = \frac{k^2}{4a^2m^2} + \bar{\rho} \frac{d\bar{\Phi}_1}{d\bar{\rho}}. \quad (37)$$

The first term comes from the quantum potential and only plays a role at short distances. This will be crucial in what follows.

### 2. Exponential instability

As long as $d\bar{\Phi}_1/d\bar{\rho} > 0$, the only destabilizing force is gravity, which is negligible in the regime we consider here, i.e., at large density and on short distances, and is only important at very large scales. However, when $d\bar{\Phi}_1/d\bar{\rho} < 0$, the self-interactions lead to an attractive force that can dominate on intermediate scales, as also noticed in [50]. Indeed, the quantum pressure always becomes dominant on very small scales, which are thus stabilized. On very large scales, gravity plays a role too. In this section, we investigate scenarios where $d\Phi_1/d\rho > 0$ at high densities and $d\Phi_1/d\rho < 0$ at low densities. Then, at high redshifts with a large background density $\bar{\rho}$, the system is stable, apart from the slow gravitational instability on large scales, i.e., the Jean's instability, and the scalar field remains homogeneous. At lower densities, the self-interactions become attractive and destabilize the system, with a fast growth of perturbations on
intermediate scales. We denote by the subscript $c_s$ the scalar background density and the scale factor when $d\Phi_1/d\rho$ changes sign to become negative,

$$ \frac{d\Phi_1}{d\rho}(\rho_{c_s}) = 0, \quad a_{c_s} = \left( \frac{3\Omega_{\text{m}0} M_{\text{Pl}}^2 H_0^2}{\rho_{c_s}} \right)^{1/3}, \quad (38) $$

where we used $\bar{\rho} \propto a^{-3}$ from that period until today. For the simple polynomial case (30), this density is given by

$$ \rho_{c_s} = \frac{c_1}{2c_2} \rho_{\Lambda}. \quad (39) $$

More generally, this change of slope of $d\Phi_1/d\rho$ will occur at a characteristic density $\rho_{\Lambda}$ that governs the self-interaction potential $V_1$. Typically, as in the polynomial case (30), shortly after the time $t_{c_s}$, e.g., after the Universe has expanded by a factor 2, $d\Phi_1/d\rho$ will be nonzero and of the order of $\Phi_{1c_s}/\rho_{c_s}$,

$$ a \sim 2a_{c_s}: \bar{\rho} \frac{d\bar{\Phi}_1}{d\bar{\rho}} \sim -|\Phi_{1c_s}| \ll 1, \quad (40) $$

where the last constraint is the nonrelativistic condition (21). For the polynomial case (30), this reads $|\Phi_{1c_s}| \sim c_1 \ll 1$. When $a \gtrsim a_{c_s}$, and considering time scales that are short compared to the Hubble time and neglecting gravity, we obtain exponential growing and decaying modes on intermediate scales, $\delta_{\pm} \propto e^{\pm\gamma_k t}$, with

$$ q < q_{\text{up}}: \gamma_q = \frac{q}{2m} \sqrt{q_{\text{up}}^2 - q^2}, \quad (41) $$

where we have introduced the upper unstable wave number

$$ \frac{d\bar{\Phi}_1}{d\bar{\rho}} < 0: \quad q_{\text{up}} = 2m\sqrt{-\bar{\rho}\frac{d\bar{\Phi}_1}{d\bar{\rho}}}, \quad (42) $$

and we denote by $q = k/a$ the physical wave number. The maximum growth rate $\gamma_{\text{max}}$ is reached at the wave number $q_{\text{max}}$, with

$$ q_{\text{max}} = \frac{q_{\text{up}}}{\sqrt{2}}, \quad \gamma_{\text{max}} = \frac{q_{\text{up}}^2}{4m} = m \left| \bar{\rho} \frac{d\bar{\Phi}_1}{d\bar{\rho}} \right|. \quad (43) $$

Therefore, wave numbers around $q_{\text{max}}$ become nonlinear first, as long as the initial power spectrum is not too steep, and structures of physical size $r \sim 2\pi/q_{\text{max}}$ appear. This perturbative growth of the scalar density perturbations is illustrated by the second column in Fig. 1. Then, shortly after this time $t_{c_s}$, the system fragments into clumps of size

$$ r_{\text{NL}} \sim \frac{2\pi}{q_{\text{max}}} \sim \frac{1}{m\sqrt{|\Phi_{1c_s}|}} \gg \frac{1}{m}, \quad (44) $$

and typical density of the order of $\rho_{\text{NL}} \sim \rho_{c_s}$, with a mass

$$ M_{\text{NL}} \sim \frac{\rho_{c_s}}{m^3 |\Phi_{1c_s}|^{3/2}}, \quad \rho_{\text{NL}} \sim \rho_{c_s} \sim \rho_{\Lambda}. \quad (45) $$

Here the subscript “NL” refers to the fact that these are the first scalar-field structures to reach the nonlinear regime, in terms of the density contrast $\delta \sim 1$.

### 3. Constraints from the linear stage

At the redshift $z_{c_s}$, assuming a standard inflationary scenario with adiabatic initial conditions, the linear density contrast on subhorizon scales during the radiation era is [51]

$$ \delta \sim -9\Phi_i \ln \left( \frac{k\eta}{\sqrt{3}} \right) = -9\Phi_i \ln \left( \frac{q}{\sqrt{3}H} \right). \quad (46) $$

This holds before the onset of the exponential instability and beyond the quantum pressure scale, which stops the logarithmic growth. Here, $\eta$ is the conformal time, with $d\eta = dt/a$, and the initial amplitude is of the order of $\Phi_i \sim 10^{-5}$. Therefore, the exponential instability (43) reaches the nonlinear regime in less than a Hubble time provided we have

$$ e^{\gamma_{\text{max}}/H} > 10^5, \quad \text{hence } \gamma_{\text{max}} > 12H_{c_s}. \quad (47) $$

Thus, we obtain the constraint that the growth rate is much greater than the Hubble expansion rate, $\gamma_{\text{max}} \gg H$, which reads

$$ m |\bar{\rho} \frac{d\bar{\Phi}_1}{d\bar{\rho}}| \gg H, \quad \text{hence } m |\Phi_{1c_s}| \gg H_{c_s}. \quad (48) $$

This gives a constraint on the parameters $m$ and $\rho_{c_s}$,

$$ \frac{H_0}{m} \Omega_{\gamma_0}^{1/2} \left( \frac{\rho_{c_s}}{M_{\text{Pl}}^2 H_0^2} \right)^{2/3} \ll |\Phi_{1c_s}| \ll 1, \quad (49) $$

which also reads

$$ 10^{-13} \left( \frac{m}{1 \text{ GeV}} \right)^{-1} \left( \frac{\rho_{c_s}}{1 \text{ GeV}^4} \right)^{2/3} \ll |\Phi_{1c_s}| \ll 1. \quad (50) $$

The condition (48) also ensures that we could neglect the expansion of the Universe in the analysis above and that the scalar field $\phi$ had already started fast oscillations in the zeroth-order quadratic potential $m^2\phi^2/2$, as $m \gg H$ (i.e., the slow-roll regime governed by the Hubble friction is already finished).

We have neglected gravity in this analysis. This is valid provided $\Phi \ll \Phi_1$. The typical gravitational potential associated with these scalar-field clumps is

$$ \Phi \sim \frac{GM_{\text{NL}}}{r_{\text{NL}}} \sim \frac{\rho_{c_s}}{M_{\text{Pl}}^2 m^2 |\Phi_{1c_s}|}. \quad (51) $$
Therefore, gravity is negligible during the formation of these structures if we have

$$ \Phi \ll \Phi_I: \frac{\rho_{c_s}}{M_{\text{Pl}}^2 m^2} \ll |\Phi_{I_{c_s}}|^2, \quad (52) $$

which reads

$$ 10^{-37} \left( \frac{m}{1 \text{ GeV}} \right)^{-2} \frac{\rho_{c_s}}{1 \text{ GeV}^4} \ll |\Phi_{I_{c_s}}|^2 \ll 1. \quad (53) $$

### C. Scalar-field solitons

Shortly after the entry into the nonlinear regime, the collapse of the first structures builds scalar-field clumps that can grow through collisions. We will analyse this aggregation process below in Sec. III D. However, after the scale factor $a(t)$ has increased by a factor two or so, the expansion of the Universe dilutes these scalar-field clumps. Then, they behave like isolated compact objects, such as MACHOs, and play the role of CDM particles.

In this section, we describe the way clumps, which are formed by the linear instability studied previously, eventually settle to equilibrium configurations. Of course, we cannot describe analytically the full time-dependent evolution of the scalar field, from the initial instability to the stable configurations that we find below. This would require numerical simulations that go beyond the present treatment. However, we check that the scalar-field dynamics admit static configurations, often called “solitons,” which are solutions to the equations of motion and are natural candidates for the end-point of the scalar-field structure-formation process. In particular, they correspond to minima of the total energy at fixed mass, which ensures their dynamical stability with respect to small nonlinear perturbations.

Therefore, we expect that the collapse of the first nonlinear structures, illustrated by the third column in Fig. 1, will build halos that are not too far from these solitons. Moreover, as they are later diluted by the Hubble expansion, these isolated clouds should naturally relax toward these solitons, possibly radiating a small amount of scalar waves that can be accreted by those clumps. This picture is also corroborated by a thermodynamical analysis, which we present in the Appendix A.

#### 1. Hydrostatic equilibrium as a minimum of the total energy

Neglecting the expansion of the Universe and using the fact that the velocity field is curl-free, the continuity and Euler equations (15) and (18) conserve the total energy

$$ E = E_{\text{kin}} + E_{\text{grav}} + E_I + E_Q, \quad (54) $$

where the kinetic, gravitational, self-interaction and quantum-pressure energies are given by

$$ \begin{aligned} E_{\text{kin}} &= \int d\vec{r} \rho \frac{\vec{v}^2}{2}, & E_{\text{grav}} &= \frac{1}{2} \int d\vec{r} \rho \Phi, \\ E_I &= \int d\vec{r} \nu_I, & E_Q &= \int d\vec{r} \frac{(\nabla \rho)^2}{8m^2 \rho}. \end{aligned} \quad (55) $$

Following [39,50], we can obtain the properties of isolated scalar clumps from an energy principle. Indeed, the conservation of energy implies that local minima of $E$ are dynamically stable with respect to small perturbations. This variational analysis goes beyond linear stability and infinitesimal perturbations, and we can expect isolated clumps to follow such profiles. Local minima at fixed mass $M$ are given by the equation $\delta E - \alpha \delta M = 0$, where $\alpha$ is the Lagrangian multiplier associated with the constraint of fixed mass [39,50]. For the energy (55), the first variation with respect to $\rho$ and $\vec{v}$ gives

$$ \delta \rho \frac{\vec{v}^2}{2} + \rho \vec{v} \cdot \delta \vec{v} + \delta \rho (\Phi + \Phi_I + \Phi_Q - \alpha) = 0. \quad (56) $$

This implies that $\vec{v} = 0$ and

$$ \Phi + \Phi_I + \Phi_Q = \alpha. \quad (57) $$

Thus, we recover the hydrostatic equilibrium of the Euler equation (18), $\nabla(\Phi + \Phi_I + \Phi_Q) = 0$. In the following we analyze the solutions to this equation.

#### 2. Gaussian ansatz for the radial profile

It is not possible to obtain an explicit solution of Eq. (57), but we can understand the main features of the equilibrium by minimizing the energy over a class of trial functions. Thus, as in [39,50], let us consider static Gaussian spherical density profiles at constant mass $M$,

$$ \rho(r) = \rho_c e^{-(r/R)^2}, \quad \text{with} \quad \rho_c = \frac{M}{\pi^{3/2} R^3}. \quad (58) $$

For the polynomial case (30), their energies are

$$ \begin{aligned} E_{\text{grav}} &= -\frac{G}{\sqrt{2}} M^{5/3} \rho_c^{1/3}, & E_Q &= \frac{3\pi M^{1/3} \rho_c^{2/3}}{4m^2}, \\ E_I &= M \left[ -\frac{c_1}{2^{5/2}} \frac{\rho_c}{\rho_\Lambda} + \frac{c_2}{3^{5/2}} \frac{\rho_c^2}{\rho_\Lambda^2} \right]. \end{aligned} \quad (59) $$

Let us neglect the gravitational energy, in agreement with (52). If we only had the quadratic term in $E_I$, both $E_Q$ and $E_I$ would be increasing functions of $\rho$. Then, the minimum of the energy would be at $\rho_c = 0$. Indeed, both the quantum pressure and the self-interactions would be repulsive, so that there would be no stable state and the scalar cloud would keep expanding and diluting (until gravity comes into play). Therefore, for small stable clumps to exist, the linear attractive term in $E_I$ must balance the quantum
pressure before it is dominated by the quadratic repulsive term. For $c_1 \sim c_2$ the transition between the attractive and repulsive regimes occurs at $\rho \sim \rho_\Lambda \sim \rho_{c_s}$, as in Sec. III B 1. Therefore, stable solitons exist provided $E_I \gtrsim E_Q$ at $\rho \sim \rho_\Lambda$. This gives

$$ \text{solitons exist if } M > M_{\text{min}}, \quad M_{\text{min}} \sim \frac{\rho_{\Lambda}}{c_1^{3/2} m^3}. \quad (60) $$

With $|\Phi_{I_{c_s}}| \sim c_1$, we find that the initial clumps (45) formed by the linear instability are actually of the order of the smallest stable mass $M_{\text{min}}$. We have seen in the derivation of (60) that this threshold also corresponds to a core density $\rho_c \sim \rho_\Lambda$. At higher masses, the quantum pressure becomes negligible and the density is set by the minimum of $E_I$. This gives again $\rho_c \sim \rho_\Lambda$. Therefore, for all masses above the threshold $M_{\text{min}}$ we have

$$ M > M_{\text{min}} : \rho_c \sim \rho_\Lambda, \quad R \sim \left(\frac{M}{\rho_\Lambda}\right)^{1/3}. \quad (61) $$

This also means that the total energy $E$ of these solitons is dominated by their self-interaction energy and it scales linearly with their mass,

$$ M \gg M_{\text{min}} : E_Q \ll E_I, \quad E \approx E_I \sim c_1 M. \quad (62) $$

The solitons (61) correspond to the regime III-a in Fig. 5 of the dense axion stars studied in [39], where the results (60)–(62) were also derived.

### 3. Numerical computation of the radial profile

A numerical computation of the soliton profiles confirms the analysis of the previous section. Neglecting the gravitational energy, the equation of equilibrium (57) that describes minima of the total energy at fixed mass reads

$$ \frac{d^2 y}{dx^2} + \frac{2 dy}{x dx} = 2y \left[ -y^2 + \frac{c_2}{c_1} y^4 + \tilde{\alpha} \right], \quad (63) $$

where we introduced the dimensionless variables

$$ y = \sqrt{\frac{\rho}{\rho_{\Lambda}}}, \quad x = \sqrt{c_1} m r, \quad \tilde{\alpha} = -\frac{\alpha}{c_1}. \quad (64) $$

Then, the soliton mass reads

$$ \frac{M}{M_0} = 4\pi \int_0^\infty dx x^2 y^2, \quad M_0 = \frac{\rho_\Lambda}{c_1^{3/2} m^3}. \quad (65) $$

As expected, $M_0$ also sets the order of magnitude of the lower mass threshold $M_{\text{min}}$ of Eq. (60).

We solve the boundary-value problem (63) with a double-shooting method (integrating from both boundaries and matching at an intermediate point) for given values of

<figure><img src="image_5.png" /><figcaption>FIG. 2. Radial density profile for the equilibrium (57). We show the masses $M = 70, 190$ and $660 M_0$, from left to right.</figcaption></figure>

$\tilde{\alpha}$. Next, integrating the density over the radius gives the total mass $M$ as a function of $\tilde{\alpha}$. We show in Fig. 2 the density profiles that we obtain for the masses $M = 70, 190$ and $660 M_0$, when we take $c_1 = c_2$. We find that at large masses the core density stabilizes at values of the order of $\rho_\Lambda$ while the mass grows as $R^3$ with the characteristic radius $R$. At large radii, $r \gg R$, the density shows an exponential tail, which is governed by the quantum pressure. We also find a lower value for the mass $M_{\text{min}}$ of these equilibrium solutions, with $M_{\text{min}} \approx 67 M_0$, in agreement with the scalings of Eq. (60). Thus, the numerical computation confirms the analytical predictions (60) and (61). Because the self-interaction potential selects the unique density scale $\rho_\Lambda$, which sets the scale of both the minima of $\mathcal{V}_I$ and $\Phi_I$, the equilibrium profiles have very simple properties. They show a flat core at a density of the order of $\rho_\Lambda$ and higher masses are obtained by increasing the radius, with $M \propto R^3$.

We discuss in more details in the Appendix C 1 the properties of these solitonic profiles, interpreting the differential equation (63) as the damped motion of a particle $y(x)$ with time $x$ in a potential $U(y)$. This provides another simple explanation of the behaviors found in Fig. 2.

As already advocated, within the nonrelativistic approximation, these solitons are stable configurations minimizing the energy functional (54) for a given value of the mass $M$. This is a feature of the nonrelativistic approximation, where the oscillation pulsation of the background field is $m$. For models where the field probes higher harmonics of the scalar potential, for instance in axionic cases where the term in $-\phi^4$ becomes of the same order as the quadratic term for large field values, the pulsation can vary at high enough density, leading to an instability of the solitons for large densities. This instability implies that the solitons can have a lifetime which can be much shorter than the age of the Universe [52]. This is not the case here, as the scalar field always follows harmonic oscillations at the leading order. Numerical simulations of a related case to the one presented in this paper confirm this observation [47].
## D. Aggregation

The perturbative analysis of Sec. III B shows that the fragmentation process starts at the redshift $z_{cs}$ where the squared sound speed of Eq. (37) changes sign to become negative. Moreover, the typical size $r_{NL}$ of these nonlinear structures is initially of the order of $2\pi/q_{max}$ as given by Eq. (44), their density of the order of $\rho_{\Lambda}$ and their mass $M_{NL}$ given by Eq. (45). The comparison with Eqs. (60)–(61) shows that this also corresponds to the lowest mass associated with stable solitons, as well as with their core density,

$$M_{NL} \sim M_{\min}, \quad \rho_{NL} \sim \rho_{\Lambda}. \qquad (66)$$

Therefore, we could expect these structures to relax toward stable solitons of mass of the order of $M_{\min}$. However, after formation and before gravity comes into play, these halos can grow (or be destroyed) through direct collisions. The typical peculiar velocity $\vec{v}_{NL}$ at the formation time $t_{cs}$ can be estimated from the linear theory, at its limit of validity when $\delta \sim 1$. From the linearized continuity equation (34) and the growth rate $\delta(\vec{k}, t) \propto e^{\gamma_k t}$, we obtain

$$v \sim \frac{\gamma_k \delta}{q}. \qquad (67)$$

At the onset of the formation of the clumps, when the modes of physical wave number $q_{max}$ reach the nonlinear regime, we obtain from Eqs. (42)–(43)

$$v_{NL}^2 \sim |\Phi_{I_{cs}}|. \qquad (68)$$

If we assume that the halos aggregate after each collision, and relax to a more massive equilibrium soliton with the scalings (61), their number density decreases with time as

$$\frac{dn}{dt} + 3Hn = -n^2\sigma v, \qquad (69)$$

with a cross section $\sigma \sim 4\pi R^2$ and a typical relative velocity $v$. This relies on the hydrodynamical picture, where scalar-field solitons behave as spheres of a barotropic fluid with a large pressure. In the regime where quantum pressure dominates, the wavelike nature of the system as described by the Schrödinger equation could lead to true solitonic behaviors, where the solitons cross each other (as in the one-dimensional cubic Schrödinger equation). However, in this paper we focus on a different regime where the self-interactions dominate over the quantum pressure. Thus, the bulk of the solitons and the scalings (61) are only determined by the shape of the self-interactions, while the quantum pressure only governs the low-density tail of the solitons. Then, we can expect the system to behave like a fluid rather than a set of waves. Thanks to the linear scaling with mass of the total energy (62), this aggregation

model conserves the total energy and can proceed without radiating significant scalar-field waves.

The Hubble expansion rate decreases as $H(t) \propto a^{-2}$ in the radiation era while the velocity dispersion decays as $v \propto 1/a$ with the expansion of the Universe, if we neglect for simplicity the velocity changes due to collisions. Assuming the mass distribution of the solitons remains peaked around a characteristic mass $M(t)$, we have $M(t) \propto 1/(a^3n)$ by conservation of the effective scalar-field density $\rho$ in a comoving volume. This expresses the growth of the halos as they merge and the falloff of their comoving number density. Then, the cross section grows as $\sigma^2 \propto (a^3n)^{-2/3}$. This gives for the solution of Eq. (69)

$$n(t) = n_i \left(\frac{a}{a_i}\right)^{-3} \left[1 + \frac{n_i \sigma_i v_i}{6H_i} \left(1 - \left(\frac{a_i}{a}\right)^2\right)\right]^{-3}, \qquad (70)$$

where the subscript $i$ stands for the initial condition at the formation time, $t_{cs}$. The first factor corresponds to the dilution by the expansion of the Universe and the second factor to the mergings of the clumps. At late times the comoving number density $n_c$ goes to a finite value,

$$a \gg a_i: \quad n_c = n_{ci} \left(1 + \frac{n_i \sigma_i v_i}{6H_i}\right)^{-3}, \qquad (71)$$

which corresponds to a typical size and mass of the final solitons of the order of

$$R_{\infty} = R_i \left(1 + \frac{n_i \sigma_i v_i}{6H_i}\right), \quad M_{\infty} = M_i \left(1 + \frac{n_i \sigma_i v_i}{6H_i}\right)^3. \qquad (72)$$

At the initial time, of the order of $t_{cs}$, we have from Eqs. (44) and (68), in agreement with the analysis of Sec. III C and with the relationship (66),

$$R_i \sim \frac{1}{m \sqrt{|\Phi_{I_{cs}}|}}, \quad v_i \sim \sqrt{|\Phi_{I_{cs}}|}, \quad \sigma_i \sim R_i^2, \quad n_i \sim \frac{1}{R_i^3}. \qquad (73)$$

This gives

$$\frac{n_i \sigma_i v_i}{H_i} \sim \frac{m |\Phi_{I_{cs}}|}{H_{cs}} \gg 1, \qquad (74)$$

where we used the constraint (48) associated with the exponential growth of small perturbations at $z_{cs}$. Thus, we have a significant merging and growth of the scalar clouds. Then, from Eq. (72) the typical size and mass of the scalar clumps formed at the end of the aggregation process is

$$R_{\text{clump}} \sim \frac{v_i}{H_i} \sim \frac{\sqrt{|\Phi_{I_{cs}}|}}{H_{cs}}, \quad M_{\text{clump}} \sim \frac{\bar{\rho}_{cs}}{H_{cs}^3} |\Phi_{I_{cs}}|^{3/2}. \qquad (75)$$
This size corresponds to the distance that can be traveled by an initial soliton during a Hubble time, sweeping material along the way, before the expansion of the Universe dilutes the scalar clouds and halts collisions. This ballistic approximation follows from the fact that we did not include the change of velocity after collisions (but we included the growth of the cross section with the rise of the halo mass). This is clearly an upper bound and we can expect a broad distribution of halo sizes, $R_i \le R \le R_\infty$, with a typical size at a lower value associated with Brownian-like trajectories. Because of this significant aggregation process, the scalar-cloud masses grow much beyond the threshold (60). This implies that the quantum pressure is negligible and the radial profile of the solitonic solutions (57) is close to a top-hat, as in (61). We will check in Sec. III F below that gravity remains negligible despite this growth of the soliton mass. This aggregation process and the final dilution by the Hubble expansion toward a collection of isolated dark-matter solitons are illustrated by the fourth column in Fig. 1.

In this section we have discussed the merging of the initial solitons by aggregation using an effective description based on the master equation (69). This provides a phenomenological understanding of the complex processes which occur when scalar-field configurations collide. A more precise characterisation of the dynamics of multi-soliton states and their collisions would require numerical simulations and a quantitative comparison with our effective results based on (69). Numerical studies of soliton collision have been performed in the recent past, for instance with Fuzzy Dark Matter in mind [53,54]. In the self-interacting case of interest here, semianalytic methods combined with numerical studies have been used in the case of quartic interactions [55] and also for bounded potentials with an attractive quartic behaviour for small field values [47]. The latter case is the closest to the one presented in this section. Numerically, various types of phenomena have been observed. Merging of solitons occurs as well as orbiting solitons in a binary system and even bouncing. When merging happens, a certain amount of scalar energy has been seen to be radiated away. This phenomenon was also observed in the case of negative quartic interactions [56], where it has been obtained that up to thirty percent of the initial soliton mass can be radiated away. This quantitative result has been obtained in a different part of the soliton phase diagram, i.e., where gravity is responsible for the existence and stability of the solitons whereas in our case gravity is negligible. In our case we hope that the compact nature of the solitons could lead to a reduced rate of scalar wave emission. All in all, a better quantitative understanding of the effect of this possible radiation loss on the asymptotic number of clumps is certainly important. One analytic possibility would be to include loss terms in (69) which could be fitted with numerical results. Such an analysis requires new numerical

simulations and a comparison with modified master equations with loss terms. We plan to come back to this topic in forthcoming works.

### E. No collapse into black holes

We now check that the scalar-field clumps formed in this process do not collapse eventually into black holes. This is avoided if the gravitational potential $\Phi$ at the surface of the stable solitons obtained above is weak and far in the Newtonian regime, $|\Phi| \ll 1$. From Eq. (75) we have

$$|\Phi| \sim \frac{GM_{\text{clump}}}{R_{\text{clump}}} \sim \frac{\bar{\rho}_{c_s} |\Phi_{I_{c_s}}|}{M_{\text{Pl}}^2 H_{c_s}^2} \ll |\Phi_{I_{c_s}}| \ll 1. \quad (76)$$

Here we used the fact that the scalar-field energy density is subdominant in the radiation era, so that $\bar{\rho}_{c_s} \ll M_{\text{Pl}}^2 H_{c_s}^2$ from the Friedmann equation. Therefore, the clumps are far in the weak-gravity regime and do not form black holes. This is consistent with the fact that gravity is always subdominant with respect to the scalar-field self-interactions.

### F. Parameter space

The scenario described in the previous sections leads to the formation of scalar clouds, or solitons, at times shortly after $t_{c_s}$. This is due to an exponential instability, which leads to a fragmentation of the homogeneous background and the formation of clumps of initial size (44). This is followed within a Hubble time by a strong aggregation process, where these scalar clouds merge to reach sizes up to (75). The profiles of these halos should relax to the solitonic solutions (57), which for the large final masses, $M_{\text{clump}} \gg M_{\text{min}}$, are approximately top-hats at the density $\rho_\Lambda$, from (61). These scalar clouds form the dark matter “particles” that play the role of the WIMPs or primordial black holes of other dark matter models. In this section, we derive the parameter space of the model allowed by theoretical constraints. This is shown in Fig. 3 in the $(T_{c_s}, m)$ plane, for the choices $|\Phi_{I_{c_s}}| = 10^{-5}$ and $|\Phi_{I_{c_s}}| = 10^{-8}$ (upper and lower panels).

First, we require the size of the scalar clumps to be below $R_{\text{max}} = 1$ pc, so that they remain much below the size of small galaxies and can build realistic dark matter profiles in galactic halos. The typical size of the clumps formed at the end of the aggregation phase was obtained in Eq. (75), which also reads

$$R_{\text{clump}} = \frac{|\Phi_{I_{c_s}}|^{1/2} \bar{\rho}_{c_s} M_{\text{Pl}}}{T_{c_s}^2}. \quad (77)$$

This gives the constraint

$$R_{\text{clump}} < R_{\text{max}} : T_{c_s} > \frac{|\Phi_{I_{c_s}}|^{1/4} \bar{\rho}_{c_s} M_{\text{Pl}}^{1/2}}{R_{\text{max}}^{1/2}}, \quad (78)$$
<figure><img src="image_6.png" /><figcaption>FIG. 3. The yellow shaded area is the domain of validity, in the plane $(T_{cs}, m)$, of the scenario described in this paper associated with potentials of the form (6). The upper panel shows the case $|\Phi_{I_{cs}}| = 10^{-5}$ and the lower panel the case $|\Phi_{I_{cs}}| = 10^{-8}$. From the left and turning clockwise, the constraints that delimit the allowed domain are associated with the maximum size of the clumps, the classicality condition, and the condition that the instability rate is much greater than the Hubble rate. The left red dot-dashed line is the temperature $T_{eq}$ at matter-radiation equality. The lower black dashed line is the condition for gravity to be negligible during the formation process (the linear stage of the tachyonic instability) while the left black dotted line in the lower panel is the condition for gravity to be negligible in the final nonlinear solitons (it does not appear in the upper panel as it is slightly to the left of this panel boundary). Within the region delimited by the previous conditions, we automatically have $T_{cs} > T_{eq}$ and negligible gravity.</figcaption></figure>

which is shown by the blue solid line labeled “$R_{clump} < R_{max}$” on the left in Fig. 3, with the choice $R_{max} = 1$ pc.

Second, we require the formation of the scalar clumps to occur before the time of matter-radiation equality. This ensures that we recover the standard CDM scenario at lower redshifts. Therefore, we impose the lower bound

$$T_{cs} > T_{eq}, \quad \text{with} \quad T_{eq} \approx 1 \text{ eV}, \qquad (79)$$

which is shown by the red dot-dashed line labeled “$T_{eq}$” on the left in Fig. 3. We can see that for $|\Phi_{I_{cs}}| \gtrsim 10^{-9}$ this constraint is automatically satisfied once we verify the first constraint (78), $R_{clump} < R_{max}$.

Next, we also have three theoretical self-consistency conditions. First, the condition (48) for an exponential instability gives a lower bound on the scalar-field mass $m$,

$$m|\Phi_{I_{cs}}| \gg H_{cs}: m \gg \frac{T_{cs}^2}{\sqrt{3}|\Phi_{I_{cs}}| M_{Pl}}. \qquad (80)$$

This corresponds to the orange solid line labeled “$m|\Phi_{I}| \gg H$” in Fig. 3. Here, we take a factor $10^3$ to ensure the left and right-hand sides are separated by at least three orders of magnitude.

Second, the classicality condition (2) provides an upper bound on the scalar mass $m$,

$$m \ll \rho_{cs}^{1/4} |\Phi_{I_{cs}}|^{-3/8}, \qquad (81)$$

where we used Eq. (73) for $v_i$. This can be written in terms of the temperature $T_{cs}$ as

$$\frac{\rho}{m^4 v^3} \gg 1: m \ll \frac{M_{Pl}^{1/8} H_0^{1/8} T_{cs}^{3/4}}{(3 \Omega_{\gamma 0})^{3/16} |\Phi_{I_{cs}}|^{3/8}}. \qquad (82)$$

This is shown by the green solid line labeled “$\rho \gg m^4 v^3$” in Fig. 3. Here, we again take a factor $10^3$ to ensure the left- and right-hand sides are separated by at least three orders of magnitude.

Third, we assumed that the gravitational force is negligible during the formation process. This is given by the constraint (52), which also reads

$$|\Phi| \ll |\Phi_{I_{cs}}|: m \gg \frac{H_0^{1/4} T_{cs}^{3/2}}{(3 \Omega_{\gamma 0})^{3/8} |\Phi_{I_{cs}}| M_{Pl}^{3/4}}. \qquad (83)$$

This corresponds to the black dashed line labeled “$|\Phi| \ll |\Phi_{I_{cs}}|$” in Fig. 3. We can see that it is automatically verified when the previous conditions are satisfied.

We can check that gravity remains small in the final solitons that are built after the nonlinear collapse and the aggregation stage. This is satisfied provided we have $|E_{grav}| \ll |E_I|$, where $E_{grav}$ and $E_I$ are the gravitational and self-interaction energies of the final solitons. From Eq. (75) and with $E_I \sim Mc_1 \sim M|\Phi_{I_{cs}}|$, this gives the condition

$$|E_{grav}| < |E_I|: T_{cs} > \frac{3^{1/4} M_{Pl}^{1/2} H_0^{1/2}}{8 \pi \Omega_{\gamma 0}^{3/4}}. \qquad (84)$$
This is shown by the vertical black dotted line labeled “$E_{grav} = E_I$” on the left in the lower panel in Fig. 3. This line does not depend on the choice of $|\Phi_{I_{cs}}|$ and it does not appear in the upper panel because it is located slightly to the left of this panel boundary. We can see that in both cases it is located to the left of the maximum-radius boundary (78). Therefore, the condition (84) is automatically satisfied and the solitons always remain governed by the self-interactions.

Thus, as shown in Fig. 3, the parameter space of the model takes the form of a triangle in the ($T_{cs}, m$) plane. It is delimited by the maximum clump size (78), the exponential-instability condition (80), and the classicality condition (82). The requirements that the formation occurs before the matter-radiation equality, (79), and that gravity remains small, (83), are automatically satisfied. Gravity also automatically remains small in the final scalar clumps, (84). Thus, we can see that the scalar-field mass spans the range

$$10^{-26} \text{ GeV} \lesssim m \lesssim 10 \text{ GeV}, \quad (85)$$

while the background temperature at the redshift $z_{cs}$ covers the range

$$10 \text{ eV} \lesssim T_{cs} \lesssim 10^5 \text{ GeV}. \quad (86)$$

This gives a wide range of temperatures and masses in the allowed parameter space.

### G. Mass and size of the scalar clumps

The typical size and mass of the clumps formed at the end of the aggregation phase were obtained in Eq. (75). This gave Eq. (77) for the radius, and for the mass:

$$M_{\text{clump}} = \frac{|\Phi_{I_{cs}}|^{3/2} 3^{3/4} M_{\text{Pl}}^{7/2} H_0^{1/2}}{\Omega_{f0}^{3/4} T_{cs}^3}. \quad (87)$$

The clump mass and radius are independent of the scalar-field mass $m$ and only depend on the redshift $z_{cs}$ when the tachyonic instability appears. We show in Fig. 4 the clump mass and radius as a function of $T_{cs}$. We also display the Schwarzschild radius of the clumps,

$$R_{\text{Sch}} = 2GM. \quad (88)$$

It is much smaller than the radius of the clumps, in agreement with the result (76) that the clumps are in the weak-gravity regime and do not form black holes.

We can see that the clumps cover a huge range of masses and radii, from microscopic to subgalactic scales. Thus, their mass goes from $10^{-3}$ gram up to $10^{37}$ gram ~ $10^4 M_{\odot}$, and their radius from 0.01 angstrom to 1 parsec. At low mass, their core density is of the order of $\rho \sim 10^{27}$ gram/cm$^3$, much above that of neutron stars, while at large mass it is of the order of $\rho \sim 10^{-13}$ gram/cm$^3 \sim 10^{17} \bar{\rho}_0$, which remains

<figure><img src="image_7.png" /><figcaption>FIG. 4. Upper panel: mass of the clumps as a function of the background temperature $T_{cs}$ at the onset of the tachyonic instability, for $|\Phi_{I_{cs}}| = 10^{-5}$ (upper blue solid line) and $|\Phi_{I_{cs}}| = 10^{-8}$ (lower red dashed line). Lower panel: radius of the clumps. The lower dotted lines show the Schwarzschild radius $R_{\text{Sch}}$.</figcaption></figure>

much greater than the current mean density $\bar{\rho}_0$ of the Universe. At the large-mass end, these clumps are thousand times more massive than the Sun, like the most massive stars, but have much greater radii, up to the parsec. Thus, they are similar to galactic molecular clouds and do not correspond to the standard stellar-mass MACHOs (massive compact halo objects), which are strongly constrained by microlensing observations.

### H. Evading microlensing constraints

Massive compact halo objects, such as primordial black holes, can be constrained by microlensing observations. Indeed, such MACHOs located in the Milky Way halo would cause a time-varying amplification of background stars when then cross their line of sight. Monitoring the Andromeda galaxy (M31) with the Subaru Hyper Suprime-Cam (HSC), the number of observed microlensing events has provided strong upper bounds on the abundance of primordial BH in the mass range $10^{-11} < M_{\text{BH}} < 10^{-6} M_{\odot}$ [57]. At low BH mass, the microlensing sensitivity is
<figure><img src="image_8.png" /><figcaption>FIG. 5. Characteristic radii in the lens plane. We show the clump radius $R_{\text{clump}}$ (blue solid line), the Einstein radius $R_E$ (black dashed line), and the outer impact parameter $\frac{d_L}{d_s} R_s$ of a source of one solar radius aligned with the lens (red dotted line). We take $d_L = 1$ kpc and $d_s = d_{M31} \approx 770$ kpc.</figcaption></figure>

strongly decreased by finite-source-size and wave-optics effects [57–59]. In this section, we show that the scalar-field solitons produced in our scenario only produce very small magnifications of distant stars, much below observational thresholds.

For pointlike lenses, the relevant scale in the lens plane is the Einstein radius $R_E$,

$$R_E = \sqrt{4GMd_L(1 - d_L/d_s)}, \quad (89)$$

where $d_L$ and $d_s$ are the distances from the observer to the lens and to the source [60,61]. For $d_L = 1$ kpc and $d_s = d_{M31} \approx 770$ kpc, this gives

$$R_E \approx 10^{-8} \sqrt{\frac{M_{\text{clump}}}{1 \text{ g}} \text{meter}}, \quad (90)$$

which is shown by the black dashed line in Fig. 5. We can see that at large masses $R_E$ becomes smaller than the radius $R_{\text{clump}}$ of the clumps. Therefore, in contrast with the case of primordial black holes, we must take into account the

effects associated with the finite size of these lenses and this will make high-mass solitons evade detection by microlensing. The flat red dotted line in Fig. 5 shows the impact parameter $\frac{d_L}{d_s} R_s$ of a circular source of one solar radius, $R_s = R_\odot$, which is aligned with the lens and the observer. At low clump mass, $\frac{d_L}{d_s} R_s$ is much greater than the clump size and the Einstein radius, which means that the finite size of the source plays a significant role. For the case of primordial black holes, this finite-source effect significantly decreases the lensing magnification. This implies that microlensing observations cannot constrain small black holes below $10^{22}$ gram [57–59]. This will also prevent the detection of low-mass solitons in our case. For small lenses, wave-optics effects also decrease the magnification as compared with the geometrical-optics prediction that neglects finite-lens effects. However, these wave-optics effects are subdominant and smeared out by the finite-size effects of the lens [58]. Therefore, in this paper we do not consider the subdominant wave-optics effects and focus on the dominant finite-size effects, which already reduce the microlensing magnification to a very small level.

To simplify the computation, we approximate the lens by a disk of constant surface density $\Sigma$. This should be a good approximation as the solitons have a flat core and a shallow envelope that shows a fast exponential decrease, see Fig. 2. Then, with the optical axis centered on the lens disk, we define the normalized radius $x_0$ of the lens, in the lens plane, as

$$x_0 = \frac{R_{\text{clump}}}{R_E}, \quad (91)$$

and the normalized impact parameter $y$ of a source at radius $r$ in the source plane, as

$$y = \frac{d_L r}{d_s R_E}. \quad (92)$$

In particular, the outer normalized impact parameter $y_s$, for a circular source of radius $R_s$ in the source plane that is aligned with the lens and the observer, is

$$y_s = \frac{d_L R_s}{d_s R_E}. \quad (93)$$

For such axially symmetric lenses, the lens equation is [60,61].

$$y = x - \frac{m(x)}{x}, \quad (94)$$

where the dimensionless lens mass within radius $x$ is
$$m(x) = 2 \int_{0}^{x} dx' x' \kappa(x'), \qquad (95)$$

with $\kappa$ the lens convergence. For a constant surface density disk, we have $\kappa = \frac{R_E^2}{R_{\text{clump}}^2} = 1/x_0^2$ inside the disk, and $\kappa = 0$ outside of the disk. This gives

$$x < x_0: m(x) = \frac{x^2}{x_0^2}, \quad x > x_0: m(x) = 1, \qquad (96)$$

and the lens mapping becomes [60]

$$\begin{aligned} |x| < x_0: & \quad y = x \left(1 - \frac{1}{x_0^2}\right), \\ |x| > x_0: & \quad y = x - \frac{1}{x}. \end{aligned} \qquad (97)$$

It is useful to define the quantity $y_0$ by

$$y_0 = \left| x_0 - \frac{1}{x_0} \right|. \qquad (98)$$

We show in Fig. 6 the normalized radii $x_0$, $y_0$, and $y_s$ in the lens plane, as a function of the clump mass. The inversion of the lens mapping (97) provides the position $x(y)$ of the image as a function of the position $y$ of the source. By axial symmetry, we can take $y \ge 0$. If there are several solutions $x_i(y)$, the lensing of the distant source gives rise to several images on the sky.

For $x_0 < 1$, the size of the lens is small and there can be strong lensing effects for small impact parameter. Thus, there are three images at small impact parameter [60],

$$x_0 < 1, \quad y < y_0: x_{\pm} = \frac{y \pm \sqrt{y^2 + 4}}{2}, \quad x_c = \frac{x_0^2}{x_0^2 - 1} y. \qquad (99)$$

The images $x_{\pm}$ are outside of the lens disk, the image $x_c$ is inside the disk. As light can propagate through the scalar cloud, the central image $x_c$ is a true solution. The magnifications associated with these images are

$$\begin{aligned} \mu_{\pm} &= \pm \frac{1}{4} \left[ \frac{y}{\sqrt{y^2 + 4}} + \frac{\sqrt{y^2 + 4}}{y} \pm 2 \right], \\ \mu_c &= \left( 1 - \frac{1}{x_0^2} \right)^{-2}, \end{aligned} \qquad (100)$$

and the total magnification is their sum

$$\mu = |\mu_+| + |\mu_-| + |\mu_c| = \frac{y^2 + 2}{y \sqrt{y^2 + 4}} + \left(1 - \frac{1}{x_0^2}\right)^{-2}. \qquad (101)$$

For large impact parameter, only the image $x_+$ exists

<figure><img src="image_9.png" /><figcaption>FIG. 6. Lensing radii normalized to the Einstein radius, in the lens plane. We show the normalized scalar-clump radius $x_0$ (blue solid line), the source radius $y_s$ (red dotted line) for a star of one solar radius, and the characteristic quantity $y_0 = |x_0 - 1/x_0|$. We take $d_L = 1$ kpc and $d_s = d_{\text{M31}} \approx 770$ kpc.</figcaption></figure>

$$x_0 < 1, \quad y > y_0: x_+ = \frac{y + \sqrt{y^2 + 4}}{2}. \qquad (102)$$

For $y \to \infty$ we recover $x_+ \approx y$, as the deflection angle decreases at large impact parameter. The magnification is

$$\mu = \mu_+ = \frac{1}{4} \left[ \frac{y}{\sqrt{y^2 + 4}} + \frac{\sqrt{y^2 + 4}}{y} + 2 \right], \qquad (103)$$

which goes to unity (no magnification) at large distance $y$.

For $x_0 > 1$, the size of the lens is large and there is always only one image. At small impact parameter it is inside the disk,

$$x_0 > 1, \quad y < y_0: x_c = \frac{x_0^2}{x_0^2 - 1} y, \qquad (104)$$
while at large impact parameter it is outside of the disk,

$$x_0 > 1, \quad y > y_0: \quad x_+ = \frac{y + \sqrt{y^2 + 4}}{2}. \qquad (105)$$

The total magnification is then either $\mu = \mu_c$ or $\mu = \mu_+$, with these quantities already given in (100).

As seen in Fig. 6, at small scalar-clump mass the size of the source is large, $y_s \gg y_0$, which means that finite-source effects must be taken into account. Therefore, we integrate the magnification over the surface of the source. Approximating the source as a disk of radius $y_s$ of constant surface brightness, the maximum total magnification is obtained when the source is centered on the optical axis, that is, it is aligned with the lens and the observer. This gives for the maximum total magnification

$$\bar{\mu}_0 = \frac{1}{\pi y_s^2} \int_0^{y_s} d^2 \bar{y} \mu(y) = \frac{2}{y_s^2} \int_0^{y_s} dy y \sum_i |\mu_i(y)|. \qquad (106)$$

In the last expression we used the axial symmetry and we explicitly wrote the sum over the images $i$ of the source, to include the case when there are several images. On the other hand, the magnification $\mu$ is also obtained from the Jacobian matrix $A_{ij} = \frac{\partial y_i}{\partial x_j}$ as $\mu = \frac{1}{\det(A)} = \frac{x}{y} \frac{dx}{dy}$, where we used the axial symmetry in the last expression. Therefore, the total magnification (106) also reads

$$\bar{\mu}_0 = \frac{2}{y_s^2} \sum_i \left| \int dx_i x_i \right|. \qquad (107)$$

As is well known, because gravitational light deflection does not involve emission, absorption, or frequency shift, the specific intensity and the surface brightness are not modified. Then, the magnification is the ratio of the solid angles subtended by the image and the source in the absence of lensing [60], $d^2x/d^2y$. We recover this ratio in Eq. (106), which we sum over the number of images.

From the above analysis, we have four cases associated with $x_0 \ge 1$ and $y_s \ge y_0$. We obtain

$$\begin{align*}
& x_0 < 1, \quad y_s < y_0: \quad \bar{\mu}_0 = \frac{x_{+s}^2 - x_{-s}^2 + x_{cs}^2}{y_s^2}, \\
& x_0 < 1, \quad y_s > y_0: \quad \bar{\mu}_0 = \frac{x_{+s}^2}{y_s^2}, \\
& x_0 > 1, \quad y_s < y_0: \quad \bar{\mu}_0 = \frac{x_{cs}^2}{y_s^2}, \\
& x_0 > 1, \quad y_s > y_0: \quad \bar{\mu}_0 = \frac{x_{+s}^2}{y_s^2}, \tag{108}
\end{align*}$$

where $x_{+s}$, $x_{-s}$, and $x_{cs}$ are the positions of the images associated with a source at position $y_s$. We can see in Fig. 6 that for low clump mass, where $x_0 < 1$, we have $y_s \gg y_0$. Thus, we obtain

$$\begin{align}
x_0(M_{\text{clump}}) < 1: \quad \bar{\mu}_0 &= \frac{x_{+s}^2}{y_s^2} = \left( \frac{1 + \sqrt{1 + 4/y_s^2}}{2} \right)^2 \nonumber \\
&\approx 1 + \frac{2}{y_s^2} \approx 1, \tag{109}
\end{align}$$

which is very close to unity as $y_s \gg 1$. For intermediate clump mass, where $x_0 > 1$ and $y_s > y_0$, we obtain again

$$\begin{align}
x_0(M_{\text{clump}}) > 1 \quad \text{and} \quad y_s > y_0: \nonumber \\
\bar{\mu}_0 = \frac{x_{+s}^2}{y_s^2} \approx 1 + \frac{2}{y_s^2} \approx 1, \tag{110}
\end{align}$$

which is again very close to unity as we still have $y_s \gg 1$. Finally, for large clump mass, where $x_0 > 1$ and $y_s < y_0$, we obtain

$$\begin{align}
x_0(M_{\text{clump}}) > 1 \quad \text{and} \quad y_s < y_0: \nonumber \\
\bar{\mu}_0 = \frac{x_{cs}^2}{y_s^2} = \left( \frac{x_0^2}{x_0^2 - 1} \right)^2 \approx 1 + \frac{2}{x_0^2} \approx 1, \tag{111}
\end{align}$$

which is very close to unity as $x_0 \gg 1$. Numerically, we find that $\bar{\mu}_0 - 1 < 10^{-3}$ over all clump masses. This is much below the observational threshold $\mu_T = 1.34$ [57]. Therefore, microlensing observations do not constrain the models studied in this paper. At low clump masses, this is because the finite-source effects decrease the lensing magnification. The same effect prevents the detection of small primordial black holes. At large masses, the microlensing inefficiency is due to the finite-lens effect, because the scalar-clump radius is much greater than the Einstein radius, $x_0 \gg 1$. This is different from primordial black hole scenarios, where large masses can be constrained by microlensing because the Schwarzschild radius is much smaller than the Einstein radius, as $R_E = \sqrt{2R_{\text{Sch}}d_L(1 - d_L/d_s)} \gg R_{\text{Sch}}$. Instead, our massive scalar clumps have a very large radius and are similar to galactic molecular clouds, rather than compact objects. This leads to small gravitational potential wells, hence to very small deflection angles and lensing magnifications. As noticed in Sec. III D, in a more realistic computation the clumps are expected to have a finite range of masses and radii below the peak values (75). However, the very small magnification $\bar{\mu}_0 - 1 < 10^{-3}$ ensures that our result should not change once we take into account the finite width of the clump mass function.

## IV. AXION MONODROMY

In this section we present another mechanism for the formation of clumps. In this case parametric resonance plays the main role. We consider this effect in the context of axion monodromy potentials as it is illustrated schematically in Fig. 7.
<figure><img src="image_10.png" /><figcaption>FIG. 7. The main stages of the formation of scalar dark-matter clumps for the parametric-resonance scenario (112).</figcaption></figure>

A. Cosine potential

Axions have been a long standing candidate for dark matter. In the case of the QCD axion, which arises after the breaking of the Peccei-Quinn symmetry, the potential term arises nonperturbatively and is periodic. This is a generic feature of axions or axionlike particles, where the axion field can be seen as a Goldstone mode of a globally broken symmetry. The potential terms for these scalar (or pseudoscalar) fields possesses the remaining discrete symmetry $\phi \to \phi + f$ which is associated to the original U(1) symmetry and arises either from nonperturbative effects or from soft breaking terms prior to the symmetry breaking phenomenon. Generically, this gives rise to cosine potentials. In the misalignment mechanism, the axion field is constant until the Hubble rate drops below the axion mass. Then oscillations start and the axion becomes a good dark matter candidate. In certain string scenarios, where the axions come from the integrated forms of string theory along closed cycles of the compactification manifold, the coupling of these fields to internal fluxes can give rise to additional polynomial interactions in the axion field [62]. This has been called axion monodromy as the potential does not remain periodic due to these fluxes but is shifted when $\phi \to \phi + f$. In the following we will focus on the potential [63]

$$V(\phi) = \frac{m_0^2}{2}\phi^2 + M_I^4[1 - \cos(\phi/f)], \quad \frac{M_I^4}{f^2} \ll m_0^2. \quad (112)$$

We can absorb the quadratic part of the cosine into the mass term and write $V(\phi) = \frac{m_0^2}{2}\phi^2 + V_1(\phi)$, with

$$m^2 = m_0^2 + \frac{M_I^4}{f^2} \approx m_0^2, \quad (113)$$

$$V_1(\phi) = M_I^4 \left[ 1 - \cos(\phi/f) - \frac{\phi^2}{2f^2} \right]. \quad (114)$$

For $\phi \ll f$ we recover an attractive quartic potential, with $\lambda_4 = -M_I^4/(6f^4) < 0$. In the following, we consider that the mass term dominates over the cosine interaction. Notice that this is different from [63], where the two terms have the same order of magnitude. In the nonrelativistic regime, where we average over the fast oscillations of the scalar field, the effective interaction potential becomes [33]

$$\Phi_1(\rho) = \frac{8\rho_b}{\rho_a} \left[ \frac{2J_1(\sqrt{\rho/\rho_b})}{\sqrt{\rho/\rho_b}} - 1 \right], \quad (115)$$

with

$$\rho_a = \frac{8m^4 f^4}{M_I^4}, \quad \rho_b = \frac{m^2 f^2}{2}, \quad \rho_b \ll \rho_a. \quad (116)$$

This corresponds to the integrated potential

$$V_1(\rho) = \frac{8\rho_b}{\rho_a} [-\rho + 4\rho_b - 4\rho_b J_0(\sqrt{\rho/\rho_b})], \quad (117)$$

obtained by averaging the potential $V_1(\phi)$ over the period of the fast leading-order oscillations. This gives for the squared-sound speed associated with the quantum pressure and the self-interactions

$$c_s^2 = \frac{k^2}{4a^2m^2} - \frac{8\rho_b}{\rho_a} J_2(\sqrt{\bar{\rho}/\rho_b}). \quad (118)$$

Thus, at large background densities the self-interaction contribution to the squared sound speed oscillates around zero, with increasing amplitude as time goes on. From the asymptotic behavior of the second-order Bessel function, we obtain at large background densities for this contribution

$$\bar{\rho} \gg \rho_b : c_s^2|_1 \approx \frac{8\rho_b}{\rho_a} \sqrt{\frac{2}{\pi}} \left(\frac{\bar{\rho}}{\rho_b}\right)^{-1/4} \cos\left(\sqrt{\bar{\rho}/\rho_b} - \frac{\pi}{4}\right). \quad (119)$$
As the speed of sound squared becomes negative quasiperiodically, there will be instabilities that we will spell out in the following section.

As for the tachyonic case, to facilitate the reading of the next sections, we already present in Fig. 7 the main stages of the formation of the scalar dark-matter clumps.

(1) The scalar field $\phi$ again quickly oscillates in the potential $V(\phi)$, dominated by its quadratic component with a small correction $V_1$. This self-interaction contribution now shows fast oscillations, such as the cosine in Eq. (112). Integrating out the fast leading-order oscillations of $\phi$, the nonrelativistic self-interaction potential $\Phi_1(\rho)$ defined by Eq. (19), i.e., Eq. (115) in our example, now shows oscillations with a decaying amplitude at large densities.

(2) At early times, the scalar-field density perturbations again oscillate as acoustic waves. As the background density $\bar{\rho}$ decreases with time, the amplitude of the self-interactions grows and they finally become relevant. As in the scenario presented in the first part of this paper and illustrated in Fig. 1, a tachyonic instability would develop at late times, associated with the first region connected to the origin where $d\Phi_1/d\rho < 0$ (the first significant drop of $\Phi_1$ seen in the lower left panel in Fig. 7). However, at much earlier times, still in the region where $\Phi_1(\rho)$ shows many oscillations, a parametric resonance triggered by these oscillatory features develops and amplifies the scalar density perturbations.

(3) The scalar density field then quickly reaches the nonlinear regime. Because of the intricate properties and time-dependent nature of the parametric resonance, the length scales and densities that first become nonlinear do not correspond to those associated with stable isolated structures. This suggests that the system will undergo a significant redistribution toward greater structures that can form stable solitons.

(4) The relative velocities are now rather modest and we do not expect significant collisional aggregation. Finally, the expansion of the Universe again dilutes the scalar clumps, which then behave as isolated CDM particles. At much lower redshifts, gravitational instability will again build the cosmic web and galaxies as in the standard ΛCDM scenario.

We describe in the following sections these various stages in more detail.

## B. Dynamics of the scalar density field

### 1. Acoustic oscillations of the density contrast

From Eqs. (118) and (119), the squared sound-speed $c_s^2$ becomes negative on subhorizon scales for the first time at the redshift $z_{cs}$, when $\bar{\rho} = \rho_{cs}$ with

$$ \frac{H_{cs}^2}{4m^2} = \frac{8\rho_b}{\rho_a} \sqrt{\frac{2}{\pi} \left(\frac{\rho_{cs}}{\rho_b}\right)^{-1/4}}, \qquad (120) $$

where we assumed that we are in the large-density regime (119),

$$ \rho_{cs} \gg \rho_b. \qquad (121) $$

This gives the useful relationship

$$ \frac{\rho_{cs}}{\rho_b} \sim \left(\frac{\rho_b}{\rho_a}\right)^4 \left(\frac{m}{H_{cs}}\right)^8. \qquad (122) $$

Since we have $\rho_b/\rho_a \ll 1$, Eq. (120) also implies

$$ z \le z_{cs} : \frac{H}{m} \ll 1. \qquad (123) $$

Thus, the slow-roll stage of the evolution of the scalar field $\phi$, when it was governed by the Hubble friction, finished long before $z_{cs}$ and the scalar field shows fast oscillations in its mainly quadratic potential $m^2\phi^2/2$. This justifies the effective description in terms of the hydrodynamical variables $\{\rho, \bar{v}\}$ and of the self-interaction potential $\Phi_1(\rho)$, as illustrated in the first column in Fig. 7. Equation (120) also reads

$$ H_{cs} \sim 10^{-2} \left( \frac{m}{1 \text{ GeV}} \right)^{16/19} \left( \frac{\rho_b}{\rho_a} \right)^{8/19} \left( \frac{\rho_b}{1 \text{ GeV}^4} \right)^{2/19} \text{ GeV}. \qquad (124) $$

Using the approximation (119), the evolution equation (36) of the linear density contrast reads

$$ \ddot{\delta} + \frac{1}{t} \dot{\delta} + \frac{H_{cs}^4}{4m^2} \frac{k^2}{k_{cs}^2} \left\{ \frac{k^2}{k_{cs}^2} \left(\frac{t}{t_{cs}}\right)^{-2} + \left(\frac{t}{t_{cs}}\right)^{-5/8} \times \cos \left[ \sqrt{\frac{\rho_{cs}}{\rho_b}} \left(\frac{t}{t_{cs}}\right)^{-3/4} \right] \right\} \delta = 0, \qquad (125) $$

where we introduced $k_{cs} = a_{cs}H_{cs}$. Here we used $H = 1/(2t)$ in the radiation era and in the cosine term we discarded the constant phase $-\pi/4$, which can be absorbed in a small change of $t_{cs}$ or of the origin of time. Making the change of time coordinate

$$ \eta = -\ln \left[ \frac{1}{2} \sqrt{\frac{\rho_{cs}}{\rho_b}} \left( \frac{t}{t_{cs}} \right)^{-3/4} \right] \ll -1, \qquad (126) $$

we obtain

$$ \frac{d^2\delta}{d\eta^2} + \frac{H_{cs}^2}{9m^2} \frac{k^2}{k_{cs}^2} \left[ \frac{k^2}{k_{cs}^2} + e^{11(\eta-\eta_{cs})/6} \cos(2e^{-\eta}) \right] \delta = 0. \qquad (127) $$
The time coordinate $\eta$ grows with cosmic time but it is restricted to large negative values, as the asymptotic form (119) of the Bessel function only applies when the argument of the cosine is large. As shown in Appendix B, at early times the density contrast shows acoustic oscillations with a constant amplitude, driven by the quantum pressure term. For moderate wave numbers, we obtain

$$ \frac{k^2}{k_{c_s}^2} \ll \frac{m}{H_{c_s}} e^{-\eta}, \quad \frac{H_{c_s} k}{mk_{c_s}} e^{\eta+11(\eta-\eta_{c_s})/12} \ll 1: \\ \delta(\eta) \approx \delta_i \cos \left( \frac{H_{c_s} k^2}{3mk_{c_s}^2} (\eta-\eta_i) \right) \left[ 1 + \frac{H_{c_s}^2 k^2}{36m^2 k_{c_s}^2} \times e^{2\eta+11(\eta-\eta_{c_s})/6} \cos(2e^{-\eta}) \right], \quad (128) $$

whereas for large wave numbers we have

$$ \frac{k^2}{k_{c_s}^2} \gg \frac{m}{H_{c_s}} e^{-\eta}, \quad \frac{H_{c_s} e^{\eta+11(\eta-\eta_{c_s})/6}}{m} \ll 1: \\ \delta(\eta) \approx \delta_i \cos \left( \frac{H_{c_s} k^2}{3mk_{c_s}^2} (\eta-\eta_i) \right) + \delta_i \sin \left( \frac{H_{c_s} k^2}{3mk_{c_s}^2} (\eta-\eta_i) \right) \\ \times \sin(2e^{-\eta}) \frac{H_{c_s}}{12m} e^{\eta+11(\eta-\eta_{c_s})/6}. \quad (129) $$

Because $H_{c_s} \ll m$ and $e^\eta \ll 1$, the density contrast only starts growing beyond its initial value $\delta_i \sim 10^{-5}$ long after the redshift $z_{c_s}$, at the time $t_g(k)$ with

$$ k < k_g: \quad t_g(k) = t_{g\infty} \left(\frac{k}{k_g}\right)^{-16/23}, \\ k > k_g: \quad t_g(k) = t_{g\infty}, \quad (130) $$

where we define

$$ k_g = k_{c_s} \left( \frac{\rho_{c_s}}{\rho_b} \right)^{11/68} \left( \frac{m}{H_{c_s}} \right)^{11/34} \gg k_{c_s}, \quad (131) $$

and

$$ t_{g\infty} = t_{c_s} \left( \frac{\rho_{c_s}}{\rho_b} \right)^{4/17} \left( \frac{m}{H_{c_s}} \right)^{8/17} \gg t_{c_s}. \quad (132) $$

Thus, the time $t_g(k)$ decreases at higher wave numbers, up to $k_g$. At greater wave numbers, $t_g(k) = t_{g\infty}$ is constant and fluctuations on these very small scales start growing simultaneously at $t_{g\infty}$. At that time, the argument of the cosine is of the order of

$$ |\tau_{g\infty}| \sim \sqrt{\frac{\rho_{c_s}}{\rho_b} \left(\frac{t_{g\infty}}{t_{c_s}}\right)^{-3/4}} \sim \left(\frac{\rho_{c_s}}{\rho_b}\right)^{11/34} \left(\frac{m}{H_{c_s}}\right)^{-6/17} \\ \sim \left(\frac{\rho_b}{\rho_a}\right)^{22/17} \left(\frac{m}{H_{c_s}}\right)^{38/17}, \quad (133) $$

where in the last expression we used the relation (122). This is still a large value if $\rho_b/\rho_a$ is not too small. Then, the squared-sound speed (119) can still show many oscillations as the background density decreases. If $|\tau_{g\infty}| \lesssim 1$, we have $\rho/\rho_b \lesssim 1$ and we are in the low-density regime of the self-interactions, where we approximate the Bessel functions by their low-order Taylor expansion. Then, we recover the polynomial case (30)–(31) with a tachyonic instability, associated with the negative value $\lambda_4 < 0$ of the quartic term of the potential $V_I(\phi)$ for $\phi \ll f$. Thus, we are back to the physics analyzed in the previous Secs. III B. Therefore, in the following we consider in more details the case $|\tau_{g\infty}| \gg 1$. We will see that in this scenario a parametric resonance takes place before the tachyonic instability can set in,

parametric resonance for $|\tau_{g\infty}| \gg 1$. \quad (134)

## 2. Mathieu-equation resonances

Changing time coordinate to

$$ \tau = -e^{-\eta} = -\frac{1}{2} \sqrt{\frac{\rho_{c_s}}{\rho_b} \left(\frac{t}{t_{c_s}}\right)^{-3/4}} \ll -1, \quad (135) $$

and writing $\delta(\tau)$ as

$$ \delta(\tau) = (-\tau)^{-1/2} y(\tau), \quad (136) $$

the evolution equation (127) becomes

$$ \frac{d^2 y}{d\tau^2} + [A(\tau) - 2q(\tau) \cos(2\tau)]y = 0, \quad (137) $$

where we shifted the argument of the cosine by a phase $\pi$ (corresponding to a negligible shift of $\tau$) to recover the standard sign of the Mathieu equation, and

$$ \begin{aligned} A(\tau) &= \left(\frac{1}{4} + \frac{H_{c_s}^2 k^4}{9m^2 k_{c_s}^4}\right) \frac{1}{\tau^2} \\ &= \frac{1}{4\tau^2} + \frac{4}{9} \left(\frac{\tau}{\tau_{g\infty}}\right)^{-2} \left(\frac{k}{k_g}\right)^4, \end{aligned} \quad (138) $$

$$ q(\tau) = \frac{H_{c_s}^2 k^2 (-\tau_{c_s})^{11/6}}{18m^2 k_{c_s}^2 (-\tau)^{23/6}} = \frac{2}{9} \left(\frac{\tau}{\tau_{g\infty}}\right)^{-23/6} \left(\frac{k}{k_g}\right)^2. \quad (139) $$

Here $\tau_{g\infty}$ is the value of $\tau$ at the time $t_{g\infty}$ introduced in (132), and it is of the order of (133). For wave numbers smaller than $k_g$, we can also write in terms of $\tau_g(k)$, associated with the time $t_g(k)$ of (130),
$$k < k_g: \begin{aligned} A(k, \tau) &= \frac{1}{4\tau^2} + \frac{4}{9} \left(\frac{\tau}{\tau_g(k)}\right)^{-2} \left(\frac{k}{k_g}\right)^{68/23}, \\ q(k, \tau) &= \frac{2}{9} \left(\frac{\tau}{\tau_g(k)}\right)^{-23/6}. \end{aligned} \quad (140)$$

Equation (137) has the form of a Mathieu equation with slowly-varying coefficients. The coefficients $A(t)$ and $q(t)$ grow with cosmic time as $|\tau|$ decreases. For constant coefficients $A$ and $q$, Floquet theory shows that the Mathieu equation has solutions of the form $e^{\pm i\nu\tau}P(\pm\tau)$, where $P(\tau)$ is periodic of period $\pi$ and $\nu$ is the characteristic Mathieu exponent [64,65]. When $\nu$ has a nonzero imaginary part, $\mu = |\mathrm{Im}(\nu)| \neq 0$, there is a growing and a decaying solution, $y_{\pm} \sim e^{\pm\mu\tau}$, up to an oscillating prefactor. This gives instability bands in the plane $(q, A)$ of the parameters, see [65]. These instability bands touch the A-axis, at $q=0$, at the discrete values $A_n = n^2$, where $n=1, 2, 3, ...$. Their width $\Delta A$ grows with $q$ for $q > 0$.

### a. Intermediate wave numbers

Let us first consider wave numbers of the order of $k_g$. For $t \lesssim t_{g\infty}$ we have $q \lesssim 1$ and $A \lesssim 1$. Therefore, we are in the

first stability region of the Mathieu-equation stability chart, see [65]. This agrees with the perturbative analysis of the Appendix B and the results (128)–(129). After that time, $A(t)$ and $q(t)$ grow beyond unity along a line $A \sim q^{12/23}$ and enter a first unstable region at $q \approx 0.4$. Then, $y(\tau)$ grows exponentially as $e^{\mu\tau}$, where $\mu = |\mathrm{Im}(\nu)|$ is the imaginary part of the characteristic Mathieu exponent [64,65]. Indeed, $A(\tau)$ and $q(\tau)$ evolve on a timescale given by $|\tau|$ whereas the $\cos(2\tau)$ term oscillates on the much shorter time scale $\pi \ll |\tau|$. Then, it takes a time $\Delta\tau$ for the density contrast $\delta$ to grow from $\delta_i \sim 10^{-5}$ to unity, with

$$\Delta\tau \sim \frac{5 \ln 10}{\mu}. \quad (141)$$

This is much shorter than $|\tau|$ as soon as $\mu$ is of the order of unity. Therefore, the density contrast reaches the nonlinear regime soon after the time $t_{g\infty}$, on comoving scales $x_g \sim 1/k_g$.

This case is shown by the upper right panel in Fig. 8. It clearly shows the oscillations with almost constant amplitude until the time $t_{g\infty}$ and the exponential rise shortly after $t_{g\infty}$. The density contrast quickly grows by a factor of $10^5$.

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 8. Evolution of the modulus of the density contrast $\delta(\vec{k}, \tau)$, for different wave numbers $k$, from Eq. (137). For illustration, the characteristic time $\tau_{g\infty}$ is set to $-10^3$ and the initial condition to $|\delta_i| = 10^{-5}$. For $k = 1.5k_g$ and $k = 1.9k_g$ (lower panels), the oscillations of $\delta(\tau)$ are so fast that they cannot be distinguished in the figure and the curve fills all the area under its upper envelope.</figcaption></figure>
### b. Low wave numbers

Let us now consider low wave numbers, $k \ll k_g$. From Eq. (140), we can see that they remain in the stability region $\{q \lesssim 1, A \lesssim 1\}$ until the time $t_g(k) > t_{g\infty}$, in agreement with the perturbative result (128). They enter the first unstable region with $q \approx 0.9$ and $A \approx 0$, along the $q$-axis. Since dark-matter clumps have already formed at the latest at the time $t_{g\infty}$, on the scale $x_g$, this is no longer relevant.

This case is shown by the upper left panel in Fig. 8. It clearly shows the oscillations with constant amplitude until the time $t_g(k) > t_{g\infty}$ and the exponential rise at $t_g(k)$. In agreement with the analysis above, the oscillation frequency is lower than for the case $k = k_g$ and the instability appears later.

### c. High wave numbers

Higher wave numbers have larger values of $q$ and $A$, along a line $A \sim q^2$. This mostly goes through the stable regions of the Mathieu equation, except for very narrow instability bands that are missed in the first-order perturbative result (129), apart from the first one. Indeed, the first signs of these higher instability bands appear as secular terms in higher orders of perturbation theory [66]. Nevertheless, we can check that the instability rate does not diverge at high wave numbers.

First, we can obtain a conservative lower bound on the time $t_g^{\min}$ when the density contrast becomes of order unity. Indeed, Eq. (B8) is valid at all orders (B7) of the perturbative expansion (B4) and provides an upper bound on secular terms. This gives $|\delta| < 3|\delta_i|$ until the time $t_g^{\min}$, with

$$t_g^{\min} = t_{c_s} \left( \frac{m}{H_{c_s}} \right)^{8/11} \gg t_{c_s}. \quad (142)$$

This time does not depend on the wave number. This implies that there is no ultraviolet divergence; the time when the density contrast becomes of the order of unity does not go to zero at high $k$ and remains above the finite value (142).

We can check that this agrees with estimates obtained from the stability chart of the constant-coefficients Mathieu equation. First, let us consider the behavior of large wave numbers, $k \gg k_g$, when they cross high-order instability bands. Let us recall that for large $n$, not too far from the $A$-axis, the $n$th instability band occurs at $A_n$ with an exponentially small width $\Delta A_n$ [67,68],

$$A_n \approx n^2, \quad \Delta A_n \sim \frac{8(q/4)^n}{[(n-1)!]^2} \left[ 1 - \frac{q^2}{4n^3} + \dots \right]. \quad (143)$$

At time $t$, we have for wave numbers greater than $k_g$, $A \approx \frac{4\tau_{g\infty}k^4}{9\tau^2k_g^4}$ and $q = \frac{2(-\tau_{g\infty})^{23/6}k^2}{9(-\tau)^{23/6}k_g^2}$. This gives $n \approx \frac{2\tau_{g\infty}k^2}{3\tau k_g^2}$ and we are

inside an instability band when $n$ is very close to an integer. We can check that the corrective term in the bracket in (143) is negligible for $|\tau| \ge |\tau_{g\infty}|$ and $k \gg k_g$. Then, we obtain for the width $\Delta A$ the asymptotic upper bound

$$n \gg 1: \quad \frac{\Delta A}{A} \lesssim \exp \left[ - \frac{4\tau_{g\infty}k^2}{3\tau k_g^2} \ln \frac{k}{k_g} \right]. \quad (144)$$

The time spent inside the instability region is $\frac{\Delta\tau}{|\tau|} = \frac{1}{2}\frac{\Delta A}{A}$. Therefore, with a growth exponent $\mu_n$, the density contrast grows during the time spent in the $n$th instability band by a factor

$$n \gg 1: \quad e^{\mu_n \Delta \tau_n} \lesssim \exp \left[ \mu_n |\tau_n| e^{-\frac{4\tau_{g\infty}k^2}{3\tau n k_g^2} \ln(k/k_g)} \right]. \quad (145)$$

As $\mu_n$ decreases at high $n$, we can see that the growth becomes negligible at high $k$. We can resum the cumulative growth due to the crossing of successive instability bands by a given wave number $k$. From $n \approx \frac{2\tau_{g\infty}k^2}{3\tau k_g^2}$, we obtain the crossing time $\tau_n$ of the $n$th band, $\tau_n \approx \frac{2\tau_{g\infty}k^2}{3nk_g^2}$. As $|\tau_n| \gg |\tau_{g\infty}|$ for all $n$, we can apply Eq. (145) for all $n \gg 1$. Neglecting the decrease of $\mu_n$ with $n$, we obtain the conservative estimate of the cumulative growth factor $G$ by the time $t_{g\infty}$,

$$G_{n_0,N} = \prod_{n=n_0}^{N} e^{\mu_n \Delta \tau_n} = e^{S_{n_0,N}}, \quad (146)$$

where $N = \frac{2k^2}{3k_g^2}$ is the final band reached at the time $t_{g\infty}$, $n_0 \gg 1$ is the lowest value where we can use Eq. (145), and

$$n_0 \gg 1: \quad S_{n_0,N} \lesssim \sum_{n=n_0}^{N} \mu |\tau_n| e^{-\frac{4\tau_{g\infty}k^2}{3\tau n k_g^2} \ln(k/k_g)} \\ \sim \mu |\tau_{g\infty}| \left(\frac{k}{k_g}\right)^{2-2n_0} \sum_{\ell=0}^{N-n_0} \frac{1}{\ell + n_0} \left(\frac{k}{k_g}\right)^{-2\ell}. \quad (147)$$

The sum over $\ell$ converges and the limit $N \to \infty$ provides an upper bound. This also shows that the cumulative growth is dominated by the lower bands, $n \sim n_0$. This gives

$$n_0 \gg 1: \quad G_{n_0,N} \lesssim \exp \left[ |\tau_{g\infty}| \left( \frac{k}{k_g} \right)^{2-2n_0} \right], \quad (148)$$

where we take the upper bound $\mu \lesssim 1$. Thus, the cumulative growth due to the crossing of high-order bands, for instance $n \ge 10$, decreases at high wave numbers. Therefore, there is no ultraviolet divergence due to the crossing of high-order instability bands by high wave numbers.

To estimate the growth associated with the crossing of the first few instability bands, we evaluate the growth
obtained for the first band $n = 1$, which should be the largest one. From Eqs. (138)–(139), we can see that high wave numbers, $k \gg k_g$, cross the first instability band, $n = 1$, at time $t_1 \sim t_{g\infty} (k/k_g)^{-8/3}$ with $A \simeq 1$ and $q \sim (k/k_g)^{-17/3} \ll 1$. At low $q$, the width of the first instability band is $\Delta A_1 \sim q$, with a growth rate $\mu_1 \sim q$ [69]. This gives a growth factor

$$n = 1, \quad k \gg k_g: e^{\mu_1 \Delta \tau_1} \sim \exp \left[ |\tau_{g\infty}| \left( \frac{k}{k_g} \right)^{-28/3} \right], \quad (149)$$

which again goes to unity at large $k$. Moreover, we can infer that Eq. (149) provides the extension down to $n_0 = 1$ of Eq. (148), which was only valid for large $n_0$ and neglected the decrease of $\mu$ at low $q$.

This high-wave number case is shown by the two lower panels in Fig. 8. In agreement with the analysis above, the instability appears earlier for higher $k$ but the amount of growth decreases as the instability bands are narrower with lower growth rates. Moreover, higher-order instability bands crossed at later times do not significantly change the amplitude of the density contrast. In the case $k = 1.9k_g$, shown in the lower right panel, the instability band crossed at $-\log(-\tau) \simeq -3.1$ actually leads to a small decrease of $|\delta|$. This can happen depending on the phase of the density contrast at the entry of the narrow instability band, if it starts with a greater weight on the decaying mode. At higher $k$, there is no significant change from the initial amplitude $|\delta_i| = 10^{-5}$ of the oscillations.

This analysis shows that the growth factor decreases at high wave numbers. Therefore, only a finite range of wave numbers above $k_g$ has been able to show a significant growth of the density contrast by the time $t_{g\infty}$. This agrees with the finiteness of the bound (142), which provides a lower bound for the earliest instability time of the fastest-growing mode $k$.

This linear growth of the scalar density perturbations by a parametric resonance is illustrated by the second column in Fig. 7.

### 3. Initial nonlinear scalar structures

Thus, we can conclude that the density contrast becomes of the order of unity at a time $t_{NL}$ with

$$t_{g \min} \le t_{NL} \le t_{g\infty}, \quad (150)$$

for wave numbers $k_{NL}$ somewhat greater than $k_g$. We can obtain an upper bound for the highest unstable wave number from Eq. (149), which as we explained above is not modified by the crossing of higher-order instability bands. This gives for the wave numbers where the Mathieu-equation instability bands can have some significant effect,

$$k \le k_g \left( \frac{\rho_{c_s}}{\rho_b} \right)^{33/952} \left( \frac{m}{H_{c_s}} \right)^{-9/238}. \quad (151)$$

This provides an upper bound for the wave numbers where the density contrast first becomes of the order of unity. The small exponents show that this upper bound is not many orders of magnitude greater than $k_g$.

Thus, we can consider that the density contrast reaches the nonlinear regime at times of the order of $t_{g\infty}$, on comoving scales $x_g \sim 1/k_g$. This gives a typical size for the first nonlinear structures in physical coordinates $r_{NL} \sim a_{g\infty}/k_g$, which yields

$$r_{\text{NL}} \sim \frac{1}{H_{c_s}} \left( \frac{\rho_{c_s}}{\rho_b} \right)^{-3/68} \left( \frac{m}{H_{c_s}} \right)^{-3/34} < \frac{1}{H_{c_s}}, \quad (152)$$

and a typical mass

$$M_{\text{NL}} \sim \frac{\rho_{c_s}}{H_{c_s}^3} \left( \frac{\rho_{c_s}}{\rho_b} \right)^{-33/68} \left( \frac{m}{H_{c_s}} \right)^{-33/34}. \quad (153)$$

Thus, we obtain a typical size that is somewhat smaller than $1/H_{c_s}$, but not by a great factor as the exponents in Eq. (152) are rather small. Using the relation (122) we can also write $r_{NL}$ and $M_{NL}$ as

$$r_{\text{NL}} \sim \frac{1}{m} \left( \frac{\rho_b}{\rho_a} \right)^{-3/17} \left( \frac{m}{H_{c_s}} \right)^{19/34}, \\ M_{\text{NL}} \sim \frac{\rho_b}{m^3} \left( \frac{\rho_b}{\rho_a} \right)^{35/17} \left( \frac{m}{H_{c_s}} \right)^{209/34}, \quad (154)$$

and the typical density as

$$\rho_{\text{NL}} \sim \bar{\rho}_{g\infty} \sim \rho_b \left( \frac{\rho_b}{\rho_a} \right)^{44/17} \left( \frac{m}{H_{c_s}} \right)^{76/17}. \quad (155)$$

Comparing with Eq. (133) we find

$$\rho_{\text{NL}} \sim \rho_b |\tau_{g\infty}|^2 \gg \rho_b, \quad (156)$$

as we assumed $|\tau_{g\infty}| \gg 1$ following (134).

In contrast with the polynomial case (45), the typical density $\rho_{NL}$ at the entry into the nonlinear density contrast regime is not only set by a characteristic density scale of the self-interaction potential, such as $\rho_b$ or $\rho_a$ that would play the role of $\rho_\Lambda$ in Eq. (45). This is clearly shown by the new factor $m/H_{c_s}$ that involves both the scalar-field mass $m$ and the Hubble expansion rate, which could be seen as an external parameter. This is due to the importance of the quantum pressure.

In the polynomial case, the instability was triggered by the change of sign of the self-interactions contribution to the squared speed of sound, see Eq. (38). The quantum pressure then only determined the lower bound for the
scales where the instability can develop, see Eqs. (41)–(43), introducing in this manner the length scale $1/m$. Thus, the quantum pressure only played a secondary role. This is also seen in the solitons found in the polynomial case, see Eqs. (60)–(61). The quantum pressure only sets the minimal mass $M_{min}$ and radius $R_{min}$ of the solitons, but their typical core density $\rho_{\Lambda}$ and their scaling law (61) do not involve the quantum pressure, which plays a negligible role at high masses (and only governs the low-density tails of these scalar clouds).

In contrast, in the case of the Bessel-type self-interaction potential (117), which decays at large densities in an oscillatory manner, the instability is triggered by a parametric resonance between the oscillations of the potential and the harmonic oscillator built by the combination of the scalar-field kinetic term and its quantum pressure. This harmonic oscillator corresponds to the term $\ddot{\delta} + c_s^2 \frac{k^2}{a^2} \delta$ in the equation of motion (36), where we only include the quantum pressure contribution to $c_s^2$, or to the term $y'' + Ay$ in the generalized Mathieu equation (137). Thus, instead of a tachyonic instability we have a parametric resonance. It clearly involves the interplay between the scalar-field kinetic terms, its quantum pressure, and its self-interactions. This explains the appearance of the new factor $m/H_{cs}$ in Eq. (155), as compared with Eq. (45). As could be expected, this multiplicative factor can be expressed in terms of $|\tau_{g\infty}|$ in (156), which measures the possibility for the parametric resonance to take place, and its advance before the tachyonic instability that would be found at low densities, see the discussion above (134).

### C. Scalar-field solitons

#### 1. Equilibrium profiles

As in Sec. III C, we now look for the equilibrium profiles of isolated scalar-field halos. Local minima of the energy $E$ at fixed mass $M$ are still given by Eqs. (56)–(57), which coincide with the hydrostatic equilibrium of the Euler equation (18). For the trial Gaussian density profile (58), the gravitational and quantum-pressure energies $E_{grav}$ and $E_Q$ are still given by Eq. (59), while the self-interaction energy $E_I$ reads

$$E_I^G = -\frac{8\rho_b}{\rho_a} M + \frac{8\rho_b}{\rho_a} M \mathcal{E}_I^G(\rho_c/\rho_b), \quad (157)$$

where $\rho_c$ is the density at the center of the halo and we introduced the function

$$\mathcal{E}_I^G(w) = \frac{16}{\sqrt{\pi w}} \int_0^\infty du u^2 [1 - J_0(\sqrt{w}e^{-u^2/2})]. \quad (158)$$

On the other hand, for a top-hat profile we obtain the same form (157) but with a scaling function $\mathcal{E}_I^{\text{T.H.}}(w)$ given by

$$\mathcal{E}_I^{\text{T.H.}}(w) = \frac{4}{w} [1 - J_0(\sqrt{w})]. \quad (159)$$

We display the functions $\mathcal{E}_I^G(w)$ and $\mathcal{E}_I^{\text{T.H.}}(w)$ in Fig. 9. The top-hat profile shows the decaying oscillations arising from the self-interaction potential (117). The regular Gaussian profile erases these small oscillations, through the smooth radial integration, and only shows a smooth decay. At large densities, $w \gg 1$, we have $\mathcal{E}_I^G(w) \sim (\ln w)^{3/2}/w$ and $\mathcal{E}_I^{\text{T.H.}}(w) \sim 1/w$, while at low densities we have $\mathcal{E}_I(0) = 1$.

The first term proportional to $M$ in Eq. (157) plays no role, as we consider minima of the total energy at constant mass. This had to be the case, because it originates from the linear term in $\rho$ in the self-interaction potential (117), which could be absorbed as a small correction to the quadratic term $\phi^2$ of the potential $V(\phi)$, see also Eq. (113). This corresponds to a small change of the scalar-field mass and should not alter the physics.

Neglecting logarithmic corrections, we write

$$E_I|_M \sim M \frac{\rho_b^2}{\rho_a(\rho_b + \rho_c)}, \quad (160)$$

which gives the correct asymptotes at both low and high core densities, except for numerical prefactors. Here, the subscript $|_M$ means that we have removed the irrelevant constant contribution $-\frac{8\rho_b}{\rho_a} M$. Then, looking for a minimum with respect to $\rho_c$ of the sum of the self-interaction and quantum-pressure energies, $E_I|_M + E_Q$, we obtain from Eqs. (59) and (160) that the minimum $\rho_c$ is nonzero for masses above a lower threshold $M_{min}$, with

<figure><img src="image_13.png" /><figcaption>FIG. 9. Self-interaction energy $\mathcal{E}_I(\rho_c/\rho_b)$ for a Gaussian profile (upper dashed line) and a top-hat profile (lower solid line). We can see the negligible influence of the oscillations compared to the overall decrease. The absence of minimum as seen for the Gaussian profile entails that a continuous distribution of densities can be present for clumps in axion monodromy models.</figcaption></figure>
$$M_{\min} \sim \left( \frac{\rho_b}{\rho_a} \right)^{-3/2} \frac{\rho_b}{m^3}, \qquad (161)$$

and for higher masses it scales as

$$M \gg M_{\min}: \rho_c \sim \rho_b \left( \frac{M}{M_{\min}} \right)^{2/5}. \qquad (162)$$

Below the mass $M_{\min}$ the self-interactions are not strong enough to resist the quantum pressure and the halo keeps on extending with a density that goes to zero. The existence of a critical mass $M_{\min}$ is thus common to the cosine potential (112) studied in this section and to potentials such as the polynomial case (30). This is because small halo masses require small radii for the density to be large enough for the self-interactions to become important, but small radii further increase the impact of the quantum pressure, as it involves gradients of the density.

On the other hand, the cosine potential (112) does not select a unique density $\rho_\Lambda$, up to factors of unity. This could be seen from the analysis of linear perturbations in Sec. IV B, where we obtained instabilities for a range of densities, which peak at a density $\rho_g$ that can be many orders of magnitude above the potential scale $\rho_b$. In the context of static isolated solitons, this is also seen from the self-interaction energy $E_I$ shown in Fig. 9, which does not display a unique minimum but keeps decreasing at large densities, possibly showing an infinite series of local minima along the way. Then, Eq. (162) shows that the interplay between the self-interactions and the quantum pressure select a mass-dependent typical density $\rho_c$ for the equilibrium profile. This core density grows with $M$ as $\rho_c \propto M^{2/5}$. Thus, while in the polynomial case the quantum pressure played no significant role in the soliton profiles (apart from setting their minimum mass) and their scaling law (61), becoming negligible at high masses, for the Bessel potential (117) the quantum pressure plays a key role at all soliton masses. There, the soliton profile is always set by the balance between the self-interactions and the quantum pressure.

The self-interactions and quantum-pressure energies of these solitons scale as

$$E_I \sim -\frac{\rho_b}{\rho_a} M + \frac{\rho_b}{\rho_a} M_{\min} \left( \frac{M}{M_{\min}} \right)^{3/5}, \\ E_Q \sim \frac{\rho_b}{\rho_a} M_{\min} \left( \frac{M}{M_{\min}} \right)^{3/5}. \qquad (163)$$

As in the polynomial scenario (62), the total energy $E \approx E_I$ is dominated by the self-interactions energy. However, this leading term $-\frac{\rho_b}{\rho_a} M$ does not play any role in the determination of the equilibrium profile, which is set by the balance between the self-interactions and quantum-pressure contributions associated with the subleading terms

$\propto M^{3/5}$ [but note that both contributions are positive in Eq. (163)].

Thus, the importance of the quantum pressure term for scenarios with Bessel-type potentials like (117), in contrast with the polynomial scenario (30) mostly governed by its self-interactions, appears both for the parametric resonance studied in Sec. IV B and for the isolated soliton profiles studied in this section.

## 2. Numerical computation of the radial profile

As for the polynomial case studied in Sec. III C 3, we confirm the analytical results with a numerical computation of the soliton profiles. Neglecting the gravitational energy, the equation of equilibrium (57) that describes minima of the total energy at fixed mass now reads

$$\frac{d^2 y}{dx^2} + \frac{2}{x} \frac{dy}{dx} = J_1(y) + \tilde{\alpha} y, \qquad (164)$$

where we introduced the dimensionless variables

$$y = \sqrt{\frac{\rho}{\rho_b}}, \quad x = \sqrt{\frac{32\rho_b}{\rho_a}}mr, \quad \tilde{\alpha} = -\frac{1}{2} - \alpha \frac{\rho_a}{16\rho_b}. \qquad (165)$$

The soliton mass is also given by the integral (65), but $M_0$ is now given by

$$M_0 = \left( \frac{32\rho_b}{\rho_a} \right)^{-3/2} \frac{\rho_b}{m^3}, \qquad (166)$$

which also sets the order of magnitude of the lower mass threshold $M_{\min}$ of Eq. (161).

We again solve the boundary-value problem (164) with a double-shooting method. We first show in Fig. 10 the mass—density relation of the equilibrium profiles that we obtain in this fashion. As expected we recover a mass-dependent core-density, with a slope that agrees with the analytical prediction (162). We show in Fig. 11 the soliton profiles

<figure><img src="image_14.png" /><figcaption>FIG. 10. Mass-core density relation for the spherical-equilibrium soliton profiles obtained from Eq. (164).</figcaption></figure>
obtained for three masses M. Even though the Bessel function $J_1(y)$ in Eq. (164) is not scale-free, its cosinelike oscillations are mostly erased by the smooth density profiles, as was the case for the self-interaction energy $E_I^G$ shown in Fig. 9 for the Gaussian profile ansatz. Then, the profiles obtained at these vastly different masses and densities are quite similar and regular, without significant oscillations. As for the polynomial case studied in Fig. 2, they show an exponential tail at large radii. These smooth behaviors explain why we recover the simple analytic prediction (162).

We discuss in more detail in the Appendix C 2 the properties of these solitonic profiles, interpreting again the differential equation (164) as the damped motion of a particle $y(x)$ with time $x$ in a potential $U(y)$. This provides another simple explanation for the behaviors found in Fig. 11. In particular, it clearly explains why the solitonic profiles obtained for the polynomial case in Fig. 2 and those obtained for the cosine model in Fig. 11 show different behaviors.

## D. Mergings in the nonlinear regime

### 1. Initial relaxation onto the soliton scaling law

We have seen in Sec. IV B 3 that, at the entry of the density contrast into the nonlinear regime, the first structures have a mass $M_{NL}$ and a density $\rho_{NL}$ given by Eqs. (154)–(155). Comparing with the minimum mass $M_{min}$ and the minimum density $\rho_{min} \sim \rho_b$ of the solitons found in Eqs. (161)–(162), we obtain

$$ \frac{M_{NL}}{M_{min}} \sim |\tau_{g\infty}|^{11/4} \gg 1, \quad \frac{\rho_{NL}}{\rho_{min}} \sim |\tau_{g\infty}|^2 \gg 1, \quad (167) $$

where we used Eq. (133). Thus, in contrast with the polynomial scenario (66), the first nonlinear structures are much greater than the smallest stable solitons and also have a greater density. As for the discussions below (156) and (162), this mismatch and the appearance of the factors $|\tau_{g\infty}|$ is due to the interplay between the self-interactions, the kinetic terms and the quantum pressure, which cannot be neglected in this scenario. Moreover, we find that these

initial structures {$M_{NL}, \rho_{NL}$} deviate from the soliton scaling law (162), since we obtain

$$ \rho_{NL} \sim \rho_b \left( \frac{M_{NL}}{M_{min}} \right)^{8/11} \gg \rho_b \left( \frac{M_{NL}}{M_{min}} \right)^{2/5}. \quad (168) $$

In other words, these initial structures are too dense as compared with the soliton equilibrium profiles. Therefore, they cannot relax to stable solitons without significant changes. In particular, if we consider an aggregation mechanism as in Sec. III D, we can no longer assume that they constitute the first steps of an aggregation process that evolves along the soliton scaling law (162), since this starting point itself deviates from this scaling law. The mismatch (168) is due to the fact that the static solitons are governed by the balance between the quantum pressure and the self-interactions, whereas the instability that gives rise to the first nonlinear structures (154) also involves the kinetic energy, associated with the time derivatives in the equations of motion (127) or (137).

The structures of masses $M_{NL}$ cannot expand within one Hubble time to lower their density so as to fall onto the scaling law (162), because of the conservation of mass within large comoving volumes (there is no outer space to expand into). Therefore, it is more natural to assume that they evolve toward the scaling law (162) by merging while keeping a density of the order of $\rho_{NL}$. From (167) we find that this target mass $M_i$ is

$$ M_i = M_{NL} |\tau_{g\infty}|^{9/4}, \quad \text{so that } \rho_{NL} \sim \rho_b \left( \frac{M_i}{M_{min}} \right)^{2/5}. \quad (169) $$

This also means that the radius of these clumps has grown to $R_i$ with

$$ R_i = r_{NL} |\tau_{g\infty}|^{3/4}. \quad (170) $$

We can compare this size with the initial velocity $v_{NL}$ of the structures that enter the nonlinear regime. At the time $t_{g\infty}$, we again estimate the typical velocity from the continuity equation (34), $v \sim r \frac{\partial \delta}{\partial t}$. With a growth rate $\delta \sim e^{\mu \tau}$, we obtain when $\delta \sim 1$ and with $\mu \sim 1$,

<figure><img src="image_15.png" /><figcaption>FIG. 11. Density profiles of the spherical-equilibrium solitons obtained from Eq. (164). We show the cases of soliton mass $M = 7 \times 10^6 M_0$, $2 \times 10^{10} M_0$ and $10^{15} M_0$ from left to right.</figcaption></figure>
$$v_{\text{NL}} \sim r_{\text{NL}} \frac{|\tau_{g\infty}|}{t_{g\infty}}. \qquad (171)$$

Comparing with (170), we can see that it takes less than a Hubble time for a disturbance to travel from $r_{\text{NL}}$ to $R_i$, as $R_i - r_{\text{NL}} \ll v_{\text{NL}}t_{g\infty}$. This suggests that it is indeed possible for the scalar-field structures to reach the mass $M_i$ within a Hubble time and to relax on the soliton scaling law (162). This stage is illustrated by the third column in Fig. 7. It seems difficult however to obtain a more rigorous description of this process by analytical means. A more detailed study of this nonlinear stage is left for future works using numerical simulations.

## 2. Lack of significant aggregation

As in the polynomial scenario discussed in Sec. III D, we could expect the scalar-field clumps $M_i$ formed by this relaxation process to grow further through collisions. From the analysis above, we start with the initial mass $M_i$ of Eq. (169), density $\rho_i \sim \rho_{\text{NL}}$, radius $R_i \sim (M_i/\rho_i)^{1/3}$ and the velocity $v_i \sim R_i/t_{g\infty}$. This velocity estimate is somewhat uncertain. It is somewhat lower than the initial velocity (171) at the entry into the nonlinear regime, but this is expected as the velocity of the clouds should decrease as they merge. We estimate the aggregation of the scalar clumps following the approach presented in Sec. III D. We still have $M(t) \propto 1/(a^3n)$ but the radius of the clouds no longer grows as $M^{1/3}$, because equilibrium profiles no longer remain at a constant density $\rho_\Lambda$. Instead, their characteristic density grows as $M^{2/5}$ with their mass, which means that their radius only grows as $R \propto M^{1/5}$ and their cross section as $M^{2/5}$. Then, the solution of the aggregation equation (69) becomes

$$n(t) = n_i \left(\frac{a}{a_i}\right)^{-3} \left[1 + \frac{3n_i \sigma_i v_i}{10H_i} \left(1 - \left(\frac{a_i}{a}\right)^2\right)\right]^{-5/3}, \quad (172)$$

where the subscript i stands for the initial condition at the time $t_{g\infty}$. At late times the comoving number density goes to the finite value

$$n_c = n_{c_i} \left( 1 + \frac{3n_i \sigma_i v_i}{10H_i} \right)^{-5/3}, \qquad (173)$$

which corresponds to a typical size and mass of the final halos of the order of

$$\begin{aligned} R_{\infty} &= R_i \left( 1 + \frac{3n_i \sigma_i v_i}{10H_i} \right)^{1/3}, \\ M_{\infty} &= M_i \left( 1 + \frac{3n_i \sigma_i v_i}{10H_i} \right)^{5/3}. \end{aligned} \qquad (174)$$

As compared with Eq. (72), we can see that the slower growth of the radius and cross section as the clumps merge

significantly damps the efficiency of the aggregation process. Moreover, with $n_i \sim 1/R_i^3$, $\sigma_i \sim R_i^2$ and $v_i \sim R_i/t_i$ we obtain

$$\frac{n_i \sigma_i v_i}{H_i} \sim 1, \quad \text{hence } R_{\infty} \sim R_i, \quad M_{\infty} \sim M_i. \qquad (175)$$

Therefore, there should be no significant aggregation through collisions. This is quite different from the strong aggregation process found for the polynomial case in Sec. III D. This is due to the much slower velocity, which we took as $v_i \sim R_i/t_{g\infty}$. This is related to the different formation process of the first nonlinear structures. Whereas in the polynomial scenario we had a tachyonic instability, leading to the fast formation of virialized scalar clouds with a typical velocity set by the strength of the self-interactions, $v \sim \sqrt{|\Phi_1|}$ as in Eq. (73), for the Bessel-type self-interaction potential we have a parametric resonance that is not directly set by the strength of the self-interactions, $\rho_b/\rho_a$, but by the interplay between the kinetic terms, the quantum pressure and the self-interactions, leading to a resonance between the oscillatory behavior of the self-interaction potential and the oscillations of the scalar field due to its wavelike properties (the quantum pressure term combined with the kinetic term). This leads to very different scalings, as seen by the comparison of $v_{\text{NL}}$ in Eq. (171), which explicitly involves the cosmic time $t_{g\infty}$, with $v_i$ in Eq. (73), which only involves the self-interactions strength $\Phi_1$.

On the other hand, if we take the larger initial value $v_{\text{NL}}$ of Eq. (171) instead of $R_i/t_{g\infty}$, which is greater by a factor $|\tau_{g\infty}|^{1/4}$, we obtain a more significant aggregation process with $R_{\infty} \sim R_i|\tau_{g\infty}|^{1/12}$ and $M_{\infty} \sim M_i|\tau_{g\infty}|^{5/12}$. The relatively small exponents show that these values are not so much larger than the previous estimates (175), unless $|\tau_{g\infty}|$ is huge.

Another difference from the polynomial scenario of Sec. III D is that energy is no longer conserved along the soliton scaling law. Indeed, from Eq. (163) we find for the total energy per unit mass (apart from kinetic energy)

$$\frac{E}{M} \sim -\frac{\rho_b}{\rho_a} + \frac{\rho_b}{\rho_a} \left( \frac{M}{M_{\text{min}}} \right)^{-2/5}. \qquad (176)$$

This means that the internal specific energy decreases as the solitons merge. This favors the mergings toward more massive halos but also suggests that some energy is radiated away as low-mass scalar waves. These may later form smaller objects or a continuous component, that could be accreted at later times by the solitons.

Thus, the estimate (175) is more uncertain than for the polynomial case (75). We can expect a broad range of halo masses and more complex nonlinear dynamics than for the polynomial case studied in Sec. III D. A more detailed investigation is left for future numerical simulations.
As in the case of the tachyonic instability, the solitons that are studied here evade the large-excursion instability which can happen for dense configurations, when the amplitude of the background field probes anharmonic parts of the scalar potential. In the case of axionic potentials, this instability could have interesting consequences such as the implosion of the solitonic configurations and potential detectable effects in the form of gravitational waves [40]. Here, we avoid these phenomena as the field never violates harmonicity at leading order. On the other hand, as the effective potential $\Phi_1(\rho)$ for the axion monodromy case oscillates at large-enough density, in the nonrelativistic regime that we have considered, it is plausible that some of the extrema of the energy functional that we have found are not in fact bona fide minima but local maxima of the energy. In this case, and similarly to the large-excursion case, there could be metastable or unstable solitonic configurations. We could for instance envisage that a maximum evolves toward a minimum simply by rearrangement of its field configuration or explodes under the destabilising effect of the quantum pressure. We leave these questions to future investigations.

### E. Solitons dominated by gravity

In contrast with the tachyonic scenario presented in Sec. III, the relatively weaker strength of the self-interactions in this parametric-resonance scenario implies that, for certain values of the model parameters, the solitons formed during the nonlinear stage become dominated by gravity rather than by the self-interactions.

First, if we consider the structures of mass $M_{NL}$ and density $\rho_{g\infty}$, at the entry into the nonlinear regime, we obtain from (167) and (59) that $E_{grav} \ll E_Q$ provided we have

$$m \gg \left(\frac{\rho_b}{\rho_a}\right)^{-1} \frac{H_0^{3/8} T_{g\infty}^{5/4}}{3^{5/16} \Omega_{\gamma 0}^{9/16} M_{Pl}^{5/8}}. \quad (177)$$

We will check in Sec. IV G that this condition is always satisfied, for the range of parameters that we consider.

However, as the structures merge to reach the greater mass $M_i$ of Eq. (169), their self-gravity also grows and can dominate over the self-interactions. Then, we find that gravity remains small as compared with the quantum pressure and the self-interactions for these masses $M_i$ provided we have

$$\text{negligible gravity: } m < \left(\frac{\rho_b}{\rho_a}\right)^{-1} \frac{\Omega_{\gamma 0}^{9/8} T_{g\infty}^{7/2}}{3^{7/8} M_{Pl}^{7/4} H_0^{3/4}}. \quad (178)$$

This boundary is shown by the black dotted line in Fig. 12 below and it is not satisfied for low values of the parameter $T_{g\infty}$. In this case, before they reach the mass $M_i$ of Eq. (169), the nonlinear scalar-field clumps become dominated by gravity rather than by the self-interactions. This

<figure><img src="image_16.png" /><figcaption>FIG. 12. The shaded area is the domain of validity, in the plane $(T_{g\infty}, m)$, of the scenario described in this paper associated with potentials of the form (112). The upper panel shows the case $\rho_b/\rho_a = 10^{-5}$ and the lower panel the case $\rho_b/\rho_a = 10^{-8}$. From the left and turning clockwise, the constraints that delimit the allowed domain are associated with the condition of formation before the matter-radiation equality, the classicality condition, and the parametric-resonance condition. The lower black dashed line is the condition for gravity to be negligible during the parametric-resonance stage, until the density fluctuations reach the nonlinear regime. The black dotted line labeled “$E_{grav} = E_Q$” is the boundary condition (178). Thus, the triangular parameter space is split into two allowed regions, separated by this transition line. Nonlinear scalar-field solitons in the yellow domain to the right of this line are governed by the balance between the quantum pressure and the self-interactions, whereas solitons in the magenta domain to the left of this line are governed by the balance between the quantum pressure and gravity.</figcaption></figure>

leads to a different scaling law from (162) for the resulting solitons. Thus, the balance between self-gravity and the quantum pressure yields the new scaling law

$$\rho_c \sim \frac{m^6 M^4}{M_{Pl}^6}, \quad (179)$$
for the relationship between the core density and the soliton mass. If the nonlinear structures now merge at the characteristic density $\rho_{g\infty}$ until they reach this new scaling law (179), the clump mass $M_i$ of Eq. (169) is replaced by

$$M_{i, grav} = \frac{\rho_{g\infty}^{1/4} M_{Pl}^{3/2}}{m^{3/2}}. \qquad (180)$$

Assuming that there is no significant aggregation afterwards, as the clumps are diluted by the expansion of the Universe, this also gives the order of magnitude of the final scalar-field clumps that play the role of the dark matter particles.

These final relaxation and dilution phases are illustrated by the fourth column in Fig. 7.

## F. No collapse into black holes

### 1. Solitons dominated by self-interactions

As for the tachyonic scenario, we again check that the scalar-field clumps do not collapse into black holes. We first consider the case where the condition (178) is satisfied: the solitons are dominated by the balance between the quantum pressure and the self-interactions, while gravity is negligible. Then, from Eq. (175) the gravitational potential at the surface of these stable solitons reads

$$|\Phi| \sim \frac{GM_{\infty}}{R_{\infty}} \sim \frac{\bar{\rho}_{c_s}}{M_{Pl}^2 H_{c_s}^2} \left(\frac{\rho_b}{\rho_a}\right)^{3/17} \left(\frac{m}{H_{c_s}}\right)^{-18/17} \ll 1, \quad (181)$$

as all factors in the last expression are much smaller than unity. Thus, these clumps are far in the weak-gravity regime and do not form black holes. This is again consistent with the fact that gravity is subdominant with respect to the scalar-field self-interactions, which are already weak.

### 2. Solitons dominated by gravity

However, in the regime studied in Sec. IV E when the condition (178) is violated, gravity dominates over the self-interactions and the scaling law of the solitons is changed to Eq. (179). Together with Eq. (180), this gives for the gravitational potential

$$\Phi \sim \frac{\rho_{g\infty}^{1/2}}{mM_{Pl}}. \qquad (182)$$

We will check in Sec. IV G and in Fig. 15 below that $|\Phi| \ll 1$ over the allowed parameter space delimited by other constraints (parametric-resonance condition, classically condition, ...). Therefore, in this case again, the clumps remain far in the Newtonian-gravity regime and do not form black holes.

## G. Parameter space

We now study the parameter space of this parametric-resonance model for the formation of dark matter as scalar clouds. This is shown in Fig. 12 over the $(T_{g\infty}, m)$ plane, for the choices $\frac{\rho_b}{\rho_a} = 10^{-5}$ and $\frac{\rho_b}{\rho_a} = 10^{-8}$. As compared with the tachyonic scenario considered in Sec. III F, the background temperature $T_{g\infty}$ at the formation time (when the parametric resonance is in full swing) plays the same role as $T_{c_s}$ (when the tachyonic instability appeared). The ratio $\frac{\rho_b}{\rho_a}$ [which sets the magnitude of the self-interactions as compared with the quadratic term in the scalar-field potential $V(\phi)$, see Eqs. (112) and (116)], plays the same role as $\Phi_{l_{c_s}}$ (which also set the relative magnitude of the self-interactions in the tachyonic model).

As in the tachyonic scenario, we require that the scalar clouds form before the time of matter-radiation equality, in order to recover a standard CDM scenario at low redshifts. Therefore, we impose the lower bound

$$T_{g\infty} > T_{eq}, \quad \text{with} \quad T_{eq} \approx 1 \text{ eV}, \qquad (183)$$

which is shown by the red dot-dashed line labeled “$T_{eq}$” on the left in Fig. 12.

We also have further theoretical self-consistency conditions. Again, we must satisfy the condition $m \gg H$, so that the slow-roll stage ends much before the formation of the scalar clouds and our nonrelativistic analysis is valid, far inside the oscillatory stage of the scalar field $\phi$ at the bottom of its mainly quadratic potential. Using Eqs. (132), (122), and (133), we obtain the useful relations

$$\frac{H_{c_s}}{m} \sim \left(\frac{H_{g\infty}}{m}\right)^{17/57} \left(\frac{\rho_b}{\rho_a}\right)^{16/57}, \qquad (184)$$

and

$$|\tau_{g\infty}| \sim \left(\frac{\rho_b}{\rho_a}\right)^{2/3} \left(\frac{m}{H_{g\infty}}\right)^{2/3} \sim 3^{1/3} \left(\frac{\rho_b}{\rho_a}\right)^{2/3} \frac{m^{2/3} M_{Pl}^{2/3}}{T_{g\infty}^{4/3}}. \qquad (185)$$

The relations (184) and (185) show that the conditions $|\tau_{g\infty}| \gg 1$ and $\frac{\rho_b}{\rho_a} \ll 1$ automatically ensure $m \gg H_{g\infty}$ and $m \gg H_{c_s}$. Therefore, the condition $m \gg H$ is automatically satisfied, once the parametric-resonance condition (134), $|\tau_{g\infty}| \gg 1$, is verified. Using Eq. (185), this gives the condition

$$|\tau_{g\infty}| \gg 1: m \gg \left(\frac{\rho_b}{\rho_a}\right)^{-1} \frac{T_{g\infty}^2}{\sqrt{3}M_{Pl}}. \qquad (186)$$

This is shown by the orange solid line labeled “$|\tau_{g\infty}| \gg 1$” in Fig. 12. Here, we take a factor $10^3$ to ensure that the left
and right-hand sides in Eq. (186) are separated by at least three orders of magnitude.

Second, the classicality condition (2) provides an upper bound on the scalar mass $m$,

$$m \ll \rho_{g\infty}^{1/4} v_{\text{NL}}^{-3/4}. \qquad (187)$$

From Eq. (171) we obtain

$$v_{\text{NL}} \sim \left(\frac{\rho_b}{\rho_a}\right)^{1/2} |\tau_{g\infty}|^{-1/4} \ll 1, \qquad (188)$$

which shows that velocities are indeed nonrelativistic. Then, Eq. (187) gives

$$\frac{\rho}{m^4 v^3} \gg 1: m \ll \left(\frac{\rho_b}{\rho_a}\right)^{-2/7} \frac{M_{\text{Pl}}^{2/7} H_0^{1/7} T_{g\infty}^{4/7}}{3^{1/7} \Omega_{\gamma 0}^{3/14}}. \qquad (189)$$

This is shown by the green solid line labeled “$\rho \gg m^4 v^3$” in Fig. 12. Here, we again take a factor $10^3$ to ensure the left- and right-hand sides are separated by at least three orders of magnitude.

Third, we assumed that gravity is negligible during the initial growth of the scalar-field fluctuations. The equation of motion (36) shows that this is satisfied if $4\pi G\bar{\rho} \ll c_s^2 k^2 / a^2$. From Eq. (119) the self-interaction contribution to the squared sound speed is $c_s^2|_{\text{I}} \sim \frac{\rho_b}{\rho_a} \left(\frac{\bar{\rho}}{\rho_b}\right)^{-1/4}$. At time $t_{g\infty}$, for density $\rho_{g\infty}$ and wave number $k_g$, this condition gives

$$|\Phi| \ll |\Phi_{1_{g\infty}}|: m \gg \left(\frac{\rho_b}{\rho_a}\right)^{-1} \frac{H_0^{3/8} T_{g\infty}^{5/4}}{3^{5/16} \Omega_{\gamma 0}^{9/16} M_{\text{Pl}}^{5/8}}. \qquad (190)$$

This corresponds to the black dashed line labeled “$|\Phi| \ll |\Phi_{1_{g\infty}}|$” in Fig. 12. We can see that it is automatically verified when the previous conditions are satisfied. As expected, Eq. (190) coincides with the condition (177) that ensures that gravity is still negligible at the entry into the nonlinear regime.

Thus, as shown in Fig. 12, the parameter space of the model takes the form of a triangle in the $(T_{g\infty}, m)$ plane. It is delimited by the background temperature $T_{\text{eq}}$ at matter-radiation equality, (183), the parametric-resonance condition (186), and the classicality condition (189). The requirement that gravity remains small during the formation process, (190), is automatically satisfied. Thus, we can see that the scalar-field mass spans the range

$$10^{-28} \text{ GeV} \lesssim m \lesssim 10 \text{ GeV}, \qquad (191)$$

while the background temperature at the redshift $z_{g\infty}$ covers the range

$$1 \text{ eV} \lesssim T_{g\infty} \lesssim 10^6 \text{ GeV}. \qquad (192)$$

As for the tachyonic scenarios, this gives a wide range of temperatures and masses in the allowed parameter space.

In contrast with the tachyonic scenarios, although gravity is always negligible during the parametric-resonance stage, where the density fluctuations grow until they reach the nonlinear regime, gravity can become dominant in the final solitons that form after the nonlinear collapse and the relaxation toward the soliton scaling laws. This is the new phenomenon studied in Sec. IV E: for scalar masses above the threshold (178) the clumps formed at the end of the nonlinear stage are dominated by gravity. The transition between the regimes where gravity is negligible or dominant with respect to the self-interactions in the final clumps is shown by the black dotted line labeled “$E_{\text{grav}} = E_Q$”, given by Eq. (178). This divides the triangle of the allowed parameter space in the $(T_{g\infty}, m)$ plane in two parts. In the right part, shown by the yellow shaded area, the final solitons are governed by the balance between the quantum pressure and the self-interactions. In the left part, shown by the magenta shaded area, the final solitons are governed by the balance between the quantum pressure and their self-gravity.

We can also check that the scalar-field clumps do not form black holes. We have seen in Sec. IV F 1 that this is guaranteed by Eq. (181) when the solitons are governed by the balance between the self-interactions and the quantum pressure, i.e., to the right of the black dotted line “$E_{\text{grav}} = E_Q$” in Fig. 12. For models to the left of this transition line, the soliton self-gravity dominates over the scalar-field self-interactions and their gravitational potential is given by Eq. (182). The latter remains small provided we have:

$$\Phi \ll 1 \quad \text{for } m \gg \frac{H_0^{1/4} T_{g\infty}^{3/2}}{\Omega_{\gamma 0}^{3/8} M_{\text{Pl}}^{3/4}}. \qquad (193)$$

We again checked that this boundary line is much below the shaded area in Fig. 12. Therefore, over all the allowed parameter space the scalar-field clumps do not collapse into black holes.

## H. Mass and size of the scalar clumps

In the regime dominated by the self-interactions, the mass and the size of the solitons formed at the end of the nonlinear stage are given by (175). This yields

$$\text{negligible gravity: } M_{\text{clump}} \sim \left(\frac{\rho_b}{\rho_a}\right)^{1/2} \frac{3^{1/4} M_{\text{Pl}}^{5/2} H_0^{1/2}}{\Omega_{\gamma 0}^{3/4} m T_{g\infty}}, \\ R_{\text{clump}} \sim \left(\frac{\rho_b}{\rho_a}\right)^{1/6} \frac{3^{1/3} M_{\text{Pl}}^{2/3}}{m^{1/3} T_{g\infty}^{4/3}}. \qquad (194)$$

In the regime dominated by the self-gravity, the mass and the size of the solitons formed at the end of the nonlinear stage are given by Eqs. (179)–(180). This yields
<figure><img src="image_17.png" /><figcaption>FIG. 13. Mass of the clumps as a function of the background temperature $T_{g\infty}$ at the peak of the parametric resonance, for $\frac{\rho_b}{\rho_a} = 10^{-5}$ (upper panel) and $\frac{\rho_b}{\rho_a} = 10^{-8}$ (lower panel). For a given $T_{g\infty}$ there is a wide range of possible clump masses $M_{\text{clump}}$. The yellow and magenta domains, on either side of the black dotted line “$E_{\text{grav}} = E_Q$”, correspond to the yellow and magenta domains shown in Fig. 12.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 14. Radius of the clumps as a function of the background temperature $T_{g\infty}$ at the peak of the parametric resonance, for $\frac{\rho_b}{\rho_a} = 10^{-5}$ (upper panel) and $\frac{\rho_b}{\rho_a} = 10^{-8}$ (lower panel). For a given $T_{g\infty}$ there is a wide range of possible clump radii $R_{\text{clump}}$. The yellow and magenta domains, on either side of the black dotted line “$E_{\text{grav}} = E_Q$”, correspond to the yellow and magenta domains shown in Figs. 12 and 13.</figcaption></figure>

$$ \begin{aligned} \text{negligible self-interactions: } M_{\text{clump}} &\sim \frac{M_{\text{Pl}}^{13/8} H_0^{1/8} T_{g\infty}^{3/4}}{(3\Omega_{\gamma 0})^{3/16} m^{3/2}}, \\ R_{\text{clump}} &\sim \frac{(3\Omega_{\gamma 0})^{3/16} M_{\text{Pl}}^{3/8}}{m^{1/2} H_0^{1/8} T_{g\infty}^{3/4}}. \end{aligned} \tag{195} $$

In contrast with the tachyonic case studied in Sec. III G, the mass and size of the clumps depend on the scalar mass $m$, in addition to the background temperature $T_{g\infty}$. Therefore, there is a finite range of clump mass and radius for a given $T_{g\infty}$, as displayed in Figs. 13 and 14. We obtain a deformed triangular domain, which corresponds to the domain of parameter space shown in Fig. 12. Its boundaries are again set by the background temperature $T_{\text{eq}}$ at matter-radiation

equality, (183), the parametric-resonance condition (186), and the classicality condition (189), as labeled in the figure. The black dotted line labeled “$E_{\text{grav}} = E_Q$” again divides the allowed domain into a region where self-gravity is negligible (to the right of this transition line) and a region where it is dominant (to the left). The slope of the upper and lower boundaries differs on either side of the transition because the clump mass and radius are either given by Eq. (194) or by Eq. (195).

As for the tachyonic case shown in Fig. 4, we find that the clumps cover a huge range of masses and radii, from microscopic to subgalactic scales. Thus, their mass goes from $10^{-6}$ gram up to $10^{36}$ gram ~ $10^3 M_{\odot}$, and their radius from 0.01 angstrom to 0.1 parsec. Again, the largest clumps are similar to galactic molecular clouds and do not correspond to the standard stellar-mass MACHOs (massive compact halo objects).
<figure><img src="image_19.png" /></figure>

FIG. 15. Characteristic radii in the lens plane. The range of allowed clump radii for a given mass $M_{\text{clump}}$ is the domain labeled “$R_{\text{clump}}$” delimited by the orange, green, and red dot-dashed curves. We also show the Schwarzschild radius $R_{\text{Sch}}$ (blue dot-dashed line), the Einstein radius $R_E$ (black dashed line), and the outer impact parameter $\frac{d_s}{d_e} R_s$ of a source of one solar radius aligned with the lens (red dotted line). We take $d_L = 1$ kpc and $d_s = d_{M31} \approx 770$ kpc.

<figure><img src="image_20.png" /></figure>

FIG. 16. Lensing radii normalized to the Einstein radius, in the lens plane. We show the normalized scalar-clump radius $x_0$ associated with the smallest clumps of a given mass in the allowed domain (orange solid line), and the associated quantity $y_0 = |x_0 - 1/x_0|$ (black dot-dashed line). The red dotted line is the source radius $y_s$ for a star of one solar radius. We take $d_L = 1$ kpc and $d_s = d_{M31} \approx 770$ kpc.

Because the largest radius obtained in Fig. 14 is slightly below 1 pc, the condition $R_{\text{clump}} < R_{\text{max}}$ is automatically satisfied for $R_{\text{max}} = 1$ pc. This is why we did not plot this condition in Fig. 12.

## I. Evading microlensing constraints

As for the tachyonic scenario, we now check whether these scalar clumps can be detected through microlensing observations. Considering again the lensing of a distant star of one solar radius in M31, at $d_s = M_{M31} \approx 770$ kpc, by a clump located in the Milky Way at $d_L = 1$ kpc, we show in Fig. 15 the characteristic radii in the lens plane as a function of the clump mass $M_{\text{clump}}$. The source projected radius, $\frac{d_L}{d_s} R_s \sim 10^6$ meter and the Einstein radius $R_E$ of Eq. (89) are the same as in Fig. 5. We also show the Schwarzschild radius $R_{\text{Sch}}$ of Eq. (88). As in Figs. 13 and 14, there is now a finite range of clump radii for a given clump mass. The finite allowed parameter space shown in Fig. 12 translates into the orange, green and red dot-dashed line that enclose the label “$R_{\text{clump}}$.” The breaks in the orange and green boundary lines, associated with the parametric-resonance and classicality conditions (186) and (189), are due to the transition from the self-interaction regime (194) to the self-gravity regime (195). The red dot-dashed curve, associated with the constraint (183) on the background temperature, does not show any break because it is fully in the self-gravity regime, see Fig. 12.

As compared with the tachyonic case displayed in Fig. 5, we obtain similar clump masses and radii, but with the line

$R_{\text{clump}}$ of Fig. 5 being thickened toward higher radii into a finite-size band.

First, we can see that the clump radii are always much greater than the Schwarzschild radius $R_{\text{Sch}}$. This confirms that the clumps do not form black holes, in agreement with the analysis of Sec. IV F and Eq. (193).

Second, the comparison with Fig. 5 shows that we have the same lensing properties as in the tachyonic case. The Einstein radius $R_E$ is always much smaller than either the projected source radius, $\frac{d_L}{d_s} R_s \sim 10^6$ meter, or the lens radius, $R_{\text{clump}}$. This implies that gravitational lensing effects are very small. The strongest lensing effects are obtained for clump radii along the lower boundary of the allowed domain, the orange curve associated with the parametric-resonance condition $|\tau_{g\infty}| \gg 1$ in Fig. 12. Indeed, this minimizes the decrease of the lensing magnification due to finite-lens effects. We show in Fig. 16 the normalized lensing radii $x_0$ and $y_0$ obtained along this lower boundary of the clump-radius domain. We also plot the normalized source radius $y_s$. We can see that we have the same configuration as in Fig. 6. At low clump masses, where $x_0 < 1$, we have $y_s \gg y_0$ and $y_s \gg 1$; using Eq. (109) this gives again $\bar{\mu}_0 \approx 1$. At intermediate clump masses, we have $x_0 > 1$, $y_s > y_0$ and $y_s \gg 1$; using Eq. (110) this also gives $\bar{\mu}_0 \approx 1$. At large clump masses, we have $x_0 \gg 1$ and $y_s < y_0$; using Eq. (111) this gives again $\bar{\mu}_0 \approx 1$. We found by a numerical computation that $\bar{\mu}_0 - 1 \ll 10^{-6}$ over all clump masses.

Thus, as for the clumps formed in the tachyonic scenario, the clumps formed in the parametric-resonance scenario
cannot be detected by microlensing. Again, at low clump masses this is due to the finite size of the source, which also prevents the detection of low-mass black holes, while at large masses this is due to the large size of the lens. In this regime, they are much bigger than both the Schwarzschild and the Einstein radii; these large clumps are similar to galactic molecular clouds, rather than compact objects, with shallow gravitational potential wells.

## J. Discussion

The model we have described in this section is formally equivalent to the one discussed in [63] for a different range of parameters. It is relevant to define

$$ \kappa = \frac{M_{\mathrm{I}}^4}{m_0^2 f^2} \approx 16 \frac{\rho_b}{\rho_a} \ll 1. \qquad (196) $$

The regime described in [63] corresponds to $\kappa \gtrsim 1$, where a rapid growth of the perturbations and the nonlinear evolution of the scalar field have been studied using numerical simulations. The formation of clumps has been observed and the consequences for structure formation analyzed. In this paper, we conduct a similar analysis in the $\kappa \ll 1$ regime. If the argument of the cosine interaction term were small, the model would reduce to the tachyonic instability case that we treated in the first part of the paper. On the contrary, as the argument of the cosine term $\tau_{g\infty} \gg 1$ is large, this regime is never attained and a parametric-resonance phase sets in first. In this case, the instability is slow initially and a long period of acoustic oscillations takes place before the onset of the parametric-resonance instability. Subsequently, we find that the result of this instability can only be the formation of solitons maintained in an equilibrium state by either the self-interactions or gravity. As our treatment is only analytical, we have no description of the intermediate steps, which we plan to investigate numerically in the future. Technically, we have obtained our description of the instabilities using the nonrelativistic approximation of the scalar-field dynamics. This should give an accurate picture as the velocities of the matter perturbations in the linear regime, up to its limit, are small. Similarly, the solitons are stable configurations where the fluid is at rest. In the intermediate regime, relativistic effects might be at play and a full numerical investigation needs to be performed. This is left for future work.

## V. CONCLUSION

We have studied the formation of clumps in scalar-field models of dark matter. These small clumps form at very high redshift, in the radiation era, and could be a candidate for the dark matter (in a manner similar to primordial black holes or small compact objects). We have explicitly worked in the nonrelativistic regime, where the homogeneous background density of dark matter is realised in the form of rapid

oscillations of the scalar field around the origin and self-interactions appear as small corrections to this background behavior. In the nonrelativistic regime, the dark-matter field can be described by a fluid with nontrivial pressure. The pressure comprises two terms. The first one originates from the kinetic terms of the scalar field and appears in the nonrelativistic description as a so-called quantum pressure. The second is due to the self-interactions and leads to a pressure term that is a function of the scalar-field energy density. We have shown that the fluid equations, in particular the Euler equation, develop unstable behaviors when the effective speed of sound squared becomes negative.

We have envisaged two scenarios. In the first one, the speed of sound squared becomes negative below a certain energy density, resulting in a tachyonic instability. At the field-theory level, this instability appears when the quartic term of the field potential is negative. This is similar to the case of the axions where the cosine potential changes convexity at large enough values of the field. For axions, this implies that perturbations of the scalar field have a tachyon instability at large values of the field along its background oscillations, leading to the formation of axions. Here, we show that a tachyonic instability due to the negative quartic interaction term in the potential is present in the nonrelativistic regime, where the oscillations of the scalar field are still almost harmonic. The resulting growth of the density contrast for the scalar energy density shows a fast exponential growth, which leads to a nonlinear regime where clumps with a nontrivial spherical profile emerge. These solitons have a well-defined density, which depends on the scalar potential of the scalar field. As a result, the clumps have a mass-radius relationship of the $M \sim R^3$ type.

A second scenario appears for axion monodromy models, where a dominant quadratic term for the scalar field is perturbed by cosine interactions. In this case, the density contrast shows a parametric-resonance instability and grows after a period of acoustic oscillations governed by the quantum pressure. In the nonlinear regime, the corresponding solitons have spherical profiles with a mass and a radius that are continuously distributed above a minimum mass threshold and obey a scaling law that follows from the balance between the self-interactions and the quantum pressure. This results in a mass-radius relationship $M \sim R^5$. This is reminiscent of the formation of oscillons in the relativistic regime, where a delayed formation occurs before parametric resonance takes place. Interestingly, despite gravity being always negligible during the formation mechanism, for small scalar mass and low formation redshift, gravity can eventually dominate the final relaxation toward the highly nonlinear solitons, and hence the properties of the scalar clumps after the aggregation phase. In this case, the mass-radius relationship is in $M \sim 1/R$. As a result, the final scalar clumps in the axion monodromy case can be governed by the balance of the quantum pressure with either the scalar self-interactions or the clump self-gravity.
We have been able to give an analytic description of the formation of scalar clumps using both linear and nonlinear arguments. The linear analysis shows that the instability due to the negative values of the speed of sound squared is always at the origin of the clumps considered here. We have also solved numerically for the nonlinear profiles of the final collapsed objects, which must satisfy the equations of hydrostatic equilibrium. In the case of the tachyonic instability, we also present in Appendix A a thermodynamical analysis that confirms the fragmentation of the system toward highly inhomogeneous configurations, with clumps at the characteristic density $\rho_{\Lambda}$. However, we have not followed the detailed relaxation from the entry into the nonlinear regime toward these stable spherical configurations. This would require numerical simulations which go beyond the present work and are left for future studies.

We have computed the allowed parameter space of these models and found that the formation redshift and the scalar-field mass span many orders of magnitude, $10^{-26} \text{ GeV} \lesssim m \lesssim 10 \text{ GeV}$. The dark-matter clumps formed by the scalar-field solitons also cover a huge range of scales, much beyond the usual MACHOs, as we find $10^{-3} \text{ gram} \lesssim M_{\text{clump}} \lesssim 10^3 M_{\odot}$ and 0.01 angstrom $\lesssim R_{\text{clump}} \lesssim 1$ parsec. Thus, they run from the size of atoms to that of galactic molecular clouds. Because of finite-source and finite-lens effects, we found that these dark-matter clumps are far below the detection thresholds of microlensing observations.

Scalar clumps are particularly interesting as they would be amenable to new tests of dark matter [40]. For instance, the creation of the clumps in the nonlinear regime could lead to the emission of gravitational waves [70]. Their existence could even be detected by the ultrasensitive detectors of gravitational wave experiments [71]. In the future, we intend to perform a more thorough investigation of the dynamics of nonrelativistic clump formation using numerical methods [47].

## ACKNOWLEDGMENTS

This work has been financially supported in part by the MINECO (Spain) Projects No. FIS2016-78859-P (AEI/FEDER) and No. PID2019-107394GB-I00.

## APPENDIX A: THERMODYNAMICS

In Secs. III B and III C in the main text, we have described the dynamics leading to the formation of scalar-field clumps using a three-pronged approach. We have first studied the tachyonic linear instabilities leading to the nonlinear regime, where the system can develop strong inhomogeneities. Second, we have obtained stable static equilibrium configurations. Third, we have described the aggregation process which yields the final masses and radii of the clumps. In this Appendix, we describe in this polynomial scenario a thermodynamic approach, where the transition from a smooth background to a strongly

inhomogeneous system, associated with the formation of clumps, can be seen as resulting from the thermodynamics of the dark-matter fluid and its interaction potential $\Phi_I$. A similar analysis may be envisaged for the case of the axion monodromy models. This is left for future work.

### 1. Phase diagram

We describe here how a thermodynamical analysis shows that structures at the characteristic density $\rho_{\Lambda}$ should form. We discard the expansion of the Universe and use the physical coordinate $\vec{r} = a\vec{x}$ as the spatial coordinate. The total energy $E$, conserved by the continuity and Euler equations, is given by Eq. (54), which reads

$$E = \int d\vec{r} \left[ \rho \frac{\vec{v}^2}{2} + \frac{1}{2} \rho \Phi + \nu_1 + \frac{(\nabla \rho)^2}{8m^2 \rho} \right]. \quad (A1)$$

To analyze the effect of the change of sign of the self-interactions with the density, we neglect the quantum pressure and gravity, which only come into play at very small and very large scales. In particular, we have seen in Sec. III C that they are negligible for the static equilibrium configurations of interest in the case of isolated scalar clouds (solitons). We consider static equilibrium configurations, with vanishing velocity field, within a given volume $V$. Then, each state is described by the density field $\rho(\vec{r})$. It is characterized by its mass $M$, energy $E$ and entropy $S$,

$$M = \int d\vec{r} \rho(\vec{r}), \quad E = \int d\vec{r} \nu_1[\rho(\vec{r})], \quad (A2)$$

$$S = - \int d\vec{r} \rho(\vec{r}) \ln \frac{\rho(\vec{r})}{\rho_{\Lambda}}. \quad (A3)$$

In the microcanonical ensemble, the most likely configurations are those that maximise the entropy at fixed values of the mass and the energy. This means that we look for the maximum of $S - E/T + (\mu/T)M$, where $1/T$ and $(\mu/T)$ are Lagrange multipliers associated with the energy and mass constraints. This is equivalent to the minimization of the grand potential $\Omega$ defined by

$$\Omega = E - TS - \mu M = \int d\vec{r} \left[ \nu_1 + T\rho \ln \frac{\rho}{\rho_{\Lambda}} - \mu\rho \right]. \quad (A4)$$

Formally, this also corresponds to the macrocanonical ensemble, where we also look for the minimum of the grand potential $\Omega$ where $T$ and $\mu$ are the temperature and the chemical potential. In our case, because there is no external thermal bath or reservoir of particles, we work in the microcanonical ensemble and consider a fixed volume with a given mass and energy. Then, $T$ and $\mu$ are only Lagrange multipliers. However, we will refer to $T$ and $\mu$ as the effective temperature and chemical potential in the
following, to simplify the terminology and facilitate the intuition of the behaviors that we obtain, which follow the standard properties of phase transitions.

The characteristic density $\rho_{\Lambda}$ is introduced in Eq. (A3) to make the argument of the logarithm dimensionless. Any other density could be used, as its value is irrelevant because it is degenerate with the Lagrange multiplier $\mu$ in Eq. (A4). The thermodynamical equilibrium is given by the minimum of the grand potential. This gives

$$ \frac{\delta\Omega}{\delta\rho} = 0: \Phi_1 + T[\ln(\rho/\rho_{\Lambda}) + 1] - \mu = 0, \quad (A5) $$

where we neglect gravity and the quantum pressure. In this approximation, we obtain a local equation in $\vec{r}$, except for the global constraints associated with the total mass and energy. It is convenient to introduce a reduced effective chemical potential $\hat{\mu}$ by

$$ \hat{\mu} = \mu/T - 1, \quad (A6) $$

so that the thermodynamical equilibrium reads

$$ \hat{\mu} = \ln(\rho/\rho_{\Lambda}) + \Phi_1/T. \quad (A7) $$

At fixed temperature, this implicit equation determines the density $\rho(\hat{\mu})$ as a function of $\hat{\mu}$. We show the curve $\hat{\mu}(\rho)$ for several values of $T$ in the upper panel in Fig. 17, for the polynomial case (30) with $c_1 = c_2$.

At large temperature, $T \to \infty$, the self-interactions $\Phi_1$ become negligible; the grand potential is governed by the entropy. Then, the function $\hat{\mu}(\rho)$ is monotonically increasing and Eq. (A7) gives the unique solution $\rho \approx \rho_{\Lambda}e^{\hat{\mu}}$. This implies a homogeneous system at this density. At low temperature, $T \to 0$, the self-interactions come into play. There is still a single solution to Eq. (A7), hence a homogeneous equilibrium, as long as $\hat{\mu}(\rho)$ remains a monotonic increasing function of $\rho$,

$$ \text{homogeneous: } 1 + \frac{1}{T} \rho \frac{d\Phi_1}{d\rho} > 0. \quad (A8) $$

Thus, we recover the condition (42), $\frac{d\Phi_1}{d\rho} < 0$, for instabilities and inhomogeneities to appear. More precisely, let us consider self-interaction potentials $\Phi_1(\rho)$ such that the derivative with respect to $\ln \rho$ has a finite minimum. Then, the curve $\hat{\mu}(\rho)$ becomes nonmonotonic below the critical temperature $T_c$ with

$$ T_c = \max \left( - \frac{d\Phi_1}{d\ln\rho} \right). \quad (A9) $$

For the polynomial case (30) this gives the critical temperature

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>FIG. 17. Thermodynamical diagrams for the polynomial case (30) with $c_1 = c_2$. Upper panel: curve $\hat{\mu}(\rho)$ at fixed temperature $T$ from Eq. (A7), for $T = 3T_c$, $T_c$ and $T_c/3$. Middle panel: grand potential, normalized by $\rho_{\Lambda}V$, from Eq. (A4), at the low temperature $T = T_c/3$. Lower panel: phase diagram in the plane $(T, \rho)$. At low temperature, $T < T_c$, the system splits in two phases of densities $\rho_+$ and $\rho_-$ if $\rho_- < \bar{\rho} < \rho_+$.</figcaption></figure>

$$ T_c = \frac{c_1^2}{8c_2}. \quad (A10) $$

As seen in the upper panel in Fig. 17, for $T < T_c$, there are three solutions $\rho(\hat{\mu})$ to the equation (A7), $\rho_- < \rho_m < \rho_+$, when $\hat{\mu}$ is in the range $\hat{\mu}_2 < \hat{\mu} < \hat{\mu}_1$, where $\hat{\mu}_1$ and $\hat{\mu}_2$ are the
local maximum and minimum of the curve $\hat{\rho}(\rho)$, at densities $\rho_1 < \rho_2$. From Eq. (A7), we obtain for the self-interaction potential (30)

$$ \frac{\rho_1}{\rho_\Lambda} = \frac{c_1 - \sqrt{c_1^2 - 8c_2T}}{4c_2}, \quad \frac{\rho_2}{\rho_\Lambda} = \frac{c_1 + \sqrt{c_1^2 - 8c_2T}}{4c_2}. \quad (A11) $$

Both $\rho_- < \rho_1$ and $\rho_+ > \rho_2$ are local minima of the grand potential $\Omega$ whereas $\rho_m$ is a local maximum, as seen in the middle panel in Fig. 17. The low-density minimum $\rho_-$ is not easily seen in the figure because the potential $\Omega$ is only slightly below zero. However, its presence is easily seen from the fact that the grand potential (A4) behaves as $\Omega \approx VT\rho \ln\frac{\rho}{\rho_\Lambda}$ at low densities, because $V_1(\rho)$ goes to zero as $\rho^2$. This means that $\Omega(\rho)$ is a decreasing function of $\rho$ at low densities, which implies that there is a local minimum $\rho_-$ in the middle panel in Fig. 17 at a density below $\rho_\Lambda/10$. The physical solution is the deepest minimum among $\{\rho_-, \rho_+\}$. For $\hat{\mu} \approx \hat{\mu}_2$ (upper dashed line in the middle panel in Fig. 17), close to the low-density monotonic branch, this is $\rho_-$, whereas for $\hat{\mu} \approx \hat{\mu}_1$ (lower dot-dashed line), close to the high-density monotonic branch, this is $\rho_+$. In-between these two regimes, there is a critical value $\hat{\mu}_s$ (solid line), $\hat{\mu}_2 < \hat{\mu}_s < \hat{\mu}_1$, where we make the transition from $\rho_-$ to $\rho_+$ as the values $\Omega_-$ and $\Omega_+$ of the grand potential cross each other. This gives a first-order phase transition, with a finite density jump at the critical chemical potential $\hat{\mu}_s$. In the limit of low temperature, we can obtain the analytic behavior of $\rho_-$ and $\rho_+$ at the critical chemical potential $\hat{\mu}_s$ as follows. From Eqs. (A4) and (A5), we find that the grand potential at equilibrium reads

$$ \Omega_{\text{eq}} = V[V_1 - \rho\Phi_1 - \rho T]. \quad (A12) $$

For the polynomial case (30), this gives

$$ \Omega_{\text{eq}} = V \left[ -\rho T + \frac{c_1}{2} \frac{\rho^2}{\rho_{\Lambda}} - \frac{2c_2}{3} \frac{\rho^3}{\rho_{\Lambda}^2} \right]. \quad (A13) $$

On the other hand, from Eq. (A11) we obtain the asymptotic behaviors

$$ T \to 0: \quad \frac{\rho_1}{\rho_\Lambda} \approx \frac{T}{c_1} \to 0, \quad \frac{\rho_2}{\rho_\Lambda} \to \frac{c_1}{2c_2}. \quad (A14) $$

Then, from $\rho_- < \rho_1$ we obtain $\rho_- \to 0$ and $\Omega_- \to 0$. Therefore, at the critical chemical potential $\hat{\mu}_s$, where $\Omega_- = \Omega_+$, we also have $\Omega_+ \to 0$ at low temperature. From Eq. (A13), together with the finite lower bound $\rho_+ > \rho_2$, this implies $\frac{c_1}{2} \frac{\rho_+}{\rho_\Lambda} - \frac{2c_2}{3} \frac{\rho_+^2}{\rho_\Lambda^2} \to 0$. Thus, we obtain the low-temperature asymptotes

$$ T \to 0: \quad \rho_- \to 0 \quad \text{and} \quad \rho_+ \to \rho_\infty \quad \text{at} \quad \hat{\mu}_s, \quad (A15) $$

with

$$ \rho_\infty = \frac{3c_1}{4c_2} \rho_\Lambda. \quad (A16) $$

We show the densities $\rho_-(T)$ and $\rho_+(T)$ of these two phases in the lower panel in Fig. 17. The curves agree with the asymptotic limits (A15).

Therefore, at high temperature, $T > T_c$, the system is homogeneous with the density $\rho = M/V$. At low temperature, $T < T_c$, the system shows a phase transition with a coexistence of two phases at densities $\rho_- < \rho_+$, with the chemical potential given by the critical value $\hat{\mu}_s(T)$ (to coexist the two phases must have the same value of $\Omega$). Thus, if the mean density $\bar{\rho}$ in the volume $V$ is below $\rho_-$ or above $\rho_+$, the system is homogeneous at the density $\bar{\rho}$. If we have $\rho_- < \bar{\rho} < \rho_+$, the system is inhomogeneous, with a coexistence of the two phases at densities $\rho_-$ and $\rho_+$. Their relative abundance is then given by the constraint on the total mass,

$$ M = \bar{\rho}V = \rho_-V_- + \rho_+V_+. \quad (A17) $$

Since $V_- < V$ is bounded and $\rho_- \to 0$, we find at low temperature

$$ T \to 0: \quad \rho_+V_+ \approx M, \quad V_+ \approx \frac{\bar{\rho}}{\rho_+} V \ll V. \quad (A18) $$

In other words, at low temperature and density below $\rho_\Lambda$, the system goes to a configuration where most of the volume is empty and a small fraction of the volume is at the characteristic density $\rho_\infty$. This characteristic density is slightly above the density $\rho_{c_s}$ of Eq. (39) where $d\Phi_1/d\rho$ vanishes. However, this thermodynamical analysis does not predict the size of the high-density clumps.

## 2. Evolution in the phase diagram

We now go back to the minimization problem (A4) within the context of the micro-canonical ensemble and of the cosmological scalar-field dynamics studied in the main text, in Sec. III B. Let us consider a constant scalar-field mass $M$ within a constant large comoving volume $V$, as is appropriate for cosmological dynamics. The system is homogeneous until the redshift $z_{c_s}$ where the tachyonic instability sets in and quickly leads to nonlinear density contrasts. The thermodynamical analysis above is then meant as a shortcut to predict the final state of the relaxation associated with the highly nonlinear dynamics that follow the entry into the nonlinear regime. To do so, we must find where the initial configurations and their subsequent evolution lies in the phase diagram shown by the lower panel in Fig. 17. As we are interested in times after $z_{c_s}$, and the expansion of the Universe dilutes the mean density $\bar{\rho}$ below the initial value $\rho_{c_s} \sim \rho_\Lambda$, see Eq. (39), we have that $\bar{\rho}$
becomes increasingly small as compared with $\rho_\Lambda$ and with the upper branch $\rho_+ \sim \rho_\Lambda$ of the phase diagram. To find out whether the system is in the strongly inhomogeneous region to the left or in the homogeneous region to the right of the boundary curve $\rho_-$ we need the energy of the system (indeed, the effective temperature $T$ is the Lagrange multiplier associated with the energy). From Eq. (A2), the energy that corresponds to the homogeneous configuration is

$$ \bar{\rho} \ll \rho_{\Lambda}: E(\bar{\rho}) = V V_1(\bar{\rho}) \approx -\frac{c_1}{2} \frac{\bar{\rho}}{\rho_{\Lambda}} M. \quad (A19) $$

On the other hand, the energy that corresponds to inhomogeneous configurations, with domains at $\rho_+ \simeq \rho_\infty$ from Eq. (A16) and at $\rho_- \ll \rho_\Lambda$, is

$$ E(\rho_+, \rho_-) = V_+ V_1(\bar{\rho}_+) + V_- V_1(\bar{\rho}_-) \approx -\frac{3}{16} \frac{c_1^2}{c_2} M, \quad (A20) $$

where we used Eq. (A18). Thus, we obtain $E(\rho_+, \rho_-) \ll E(\bar{\rho}) < 0$ and as expected the inhomogeneous configuration associated with low $T$ is also associated with a low energy, in our case a large negative energy. On the other hand, at the entry into the nonlinear regime at the redshift $z_{c_s}$ we have $\rho \approx \rho_{c_s} \sim \rho_\Lambda$ and the initial energy is $E_{c_s} \sim -c_1 M$. As long as gravity is negligible, that is, until gravitational clustering develops at redshifts $z \lesssim 10$, the local self-interactions associated with the potential $V_1$ conserve the energy within large comoving volumes, which are essentially independent of each other (it is simply the sum of the internal energies of the scalar-field solitons contained within each comoving volume). Therefore, we keep $E \sim -c_1 M$, which selects the inhomogeneous configuration (A20), whereas the homogeneous configuration (A19) corresponds to an increasingly far high-energy configuration, with $E(\bar{\rho}) \to 0^-$.

Thus, we can conclude that at the entry in the nonlinear regime, at $z_{c_s}$, the system is close to the upper-right point in the phase diagram shown by the lower panel in Fig. 17, where the curves $\rho_+$ and $\rho_-$ meet with $\rho_+ \sim \rho_- \sim \bar{\rho} \sim \rho_\Lambda$, and that at later times the system moves to the lower-left part of the diagram, increasingly far into the inhomogeneous region to the left of the boundary curve $\rho_-$. Hence this simple thermodynamical analysis suggests that after the tachyonic instability studied in Sec. III B 2 has reached the nonlinear regime the complex dynamics that follow will lead to a fragmentation of the system over domains of density of the order of $\rho_+ \sim \rho_\Lambda$, which contain most of the mass, and domains of density $\rho_- \ll \rho_\Lambda$, which make most of the volume. This agrees with a simple halo model where the scalar field is clustered into the stable solitons obtained in Sec. III C amidst empty space.

## APPENDIX B: TIME-DEPENDENT MATHIEU EQUATION

The evolution equation (127) reads

$$ \frac{d^2\delta}{d\eta^2} + \omega^2\delta + \epsilon e^{11(\eta-\eta_{cs})/6} \cos(2e^{-\eta})\delta = 0, \quad (B1) $$

with

$$ \omega = \frac{H_{cs} k^2}{3mk_{cs}^2}, \quad \epsilon = \frac{H_{cs}^2 k^2}{9m^2 k_{cs}^2}. \quad (B2) $$

By assumption, for the asymptotic behavior (119) of the Bessel function to be valid, we restrict ourselves to the range

$$ e^{-\eta} \gg 1. \quad (B3) $$

We typically have $\omega \ll 1$ and $\epsilon \ll 1$ as $H_{cs} \ll m$, except for very large wave numbers. We can look for a perturbative expansion in $\epsilon$ of the form

$$ \delta(\eta) = \sum_{n=0}^{\infty} \epsilon^n \delta^{(n)}(\eta). \quad (B4) $$

The zeroth-order solution is

$$ \delta^{(0)}(\eta) = \delta_i \cos[\omega(\eta - \eta_i)], \quad (B5) $$

with the initial conditions $\{\delta = \delta_i, \delta' = 0\}$ at the initial time $\eta_i$. Thus, when the self-interactions are negligible the density contrast shows acoustic oscillations of constant amplitude because of the quantum pressure term. In this regime, the density perturbations do not grow. Using for instance the method of variation of parameters or Green's function [66], we obtain the solution of Eq. (B1) up to order $n$ as

$$ \begin{aligned} \delta^{(n)}(\eta) = & \cos[(\omega(\eta - \eta_i)) \int_{\eta_i}^{\eta} d\eta' \frac{\sin[\omega(\eta' - \eta_i)]}{\omega} \\ & \times e^{11(\eta' - \eta_{cs})/6} \cos(2e^{-\eta'}) \delta^{(n-1)}(\eta') - \sin[(\omega(\eta - \eta_i))] \\ & \times \int_{\eta_i}^{\eta} d\eta' \frac{\cos[\omega(\eta' - \eta_i)]}{\omega} e^{11(\eta' - \eta_{cs})/6} \cos(2e^{-\eta'}) \\ & \times \delta^{(n-1)}(\eta'). \end{aligned} \quad (B6) $$

From this recursion it is easy to obtain the upper bound

$$ |\epsilon^n \delta^{(n)}(\eta)| \le \frac{|\delta_i|}{n!} \left( \frac{12\epsilon}{11\omega} e^{11(\eta-\eta_{cs})/6} \right)^n. \quad (B7) $$

Therefore, the perturbative expansion (B4) converges for all values of $\omega$, $\epsilon$ and $\eta$. Moreover, we have
$$|\delta| < 3|\delta_i| \quad \text{for} \quad \frac{12\epsilon}{11\omega} e^{11(\eta-\eta_{cs})/6} < 1. \qquad (B8)$$

Thus, for any wave number the solution is well described by the zeroth-order acoustic oscillations (B5) at sufficiently early times.

From Eq. (B6) we obtain the first-order correction $\delta^{(1)}$ in terms of incomplete Gamma functions. For moderate values of $\omega$, and large values of $e^{-\eta}$, this gives

$$\omega \ll e^{-\eta}: \delta^{(1)}(\eta) \approx \frac{\delta_i}{4} e^{11(\eta-\eta_{cs})/6} e^{2\eta} \cos[(\omega(\eta-\eta_i))] \times \cos(2e^{-\eta}), \qquad (B9)$$

whereas for large values of $\omega$ we obtain

$$\omega \gg e^{-\eta}: \delta^{(1)}(\eta) \approx \frac{\delta_i}{4\omega} e^{11(\eta-\eta_{cs})/6} e^{\eta} \sin[(\omega(\eta-\eta_i))] \times \sin(2e^{-\eta}). \qquad (B10)$$

We can directly check on the equation of motion (B1) that these are the first-order perturbative corrections associated with the zeroth-order term (B5) in these two regimes. The amplitudes (B9) and (B10) are smaller than the conservative upper bound (B7) by factors $e^{\eta} \ll 1$. This is due to the fast oscillating factor $\cos(2e^{-\eta})$ in the perturbative term of the equation of motion (B1), which damps its impact on the dynamics. Then, the density contrast is well described by the zeroth-order solution (B5) until the first order correction, given by either (B9) or (B10), becomes of the same order.

# APPENDIX C: SOLITON RADIAL PROFILE AS A DAMPED TRAJECTORY IN A POTENTIAL

## 1. Polynomial case

The differential equation (63) that determines the radial profile of the scalar-field solitons, defined as a minimum of the energy at fixed mass, can also be interpreted as the trajectory $y(x)$ of particle that rolls down a potential $U(y)$, with a time-dependent friction,

$$\frac{d^2 y}{dx^2} + \frac{2}{x} \frac{dy}{dx} = - \frac{dU}{dy}, \qquad (C1)$$

where $x$ is interpreted as a time coordinate, running from 0 to $+\infty$. In particular, we have

$$\frac{d}{dx} \left[ \frac{1}{2} \left( \frac{dy}{dx} \right)^2 + U(y) \right] = - \frac{2}{x} \left( \frac{dy}{dx} \right)^2 \le 0, \qquad (C2)$$

which shows that the total energy of the particle, defined as the sum of its kinetic and potential energies, decreases with the time $x$. The boundary conditions of the trajectory are $\frac{dy}{dx} = 0$ at $x = 0$, because we require a regular profile at the origin, and $y(x) = 0$ at $x \to \infty$, as the density must

<figure><img src="image_24.png" /><figcaption>FIG. 18. Trajectories in the potentials $U(y)$ for different soliton masses, corresponding to the soliton radial profiles of Fig. 2.</figcaption></figure>

decrease at large radii to obtain a finite mass. With $U(0) = 0$, this means that at late times the particle must settle to the point $y = 0$ and that it starts at $x = 0$ from a value $y_0 > 0$ with $U(y_0) > 0$ and a vanishing velocity.

For the polynomial scalar-field potential (30), associated with the differential equation (63), the effective particle potential $U(y)$ reads

$$U(y) = \frac{1}{2} y^4 - \frac{c_2}{3c_1} y^6 - \tilde{\alpha} y^2. \qquad (C3)$$

It depends on the unknown parameter $\tilde{\alpha}$, which is a function of the soliton mass $M$. This parameter $\tilde{\alpha}$ is strictly positive so that the density shows an exponential tail at large radii with $y \sim e^{-\sqrt{2\tilde{\alpha}x}}$.

We show in Fig. 18 the potentials $U(y)$ for the soliton profiles displayed in Fig. 2, over the range $0 \le y \le y_0$ covered by the particle as it rolls down its potential from the starting point $y_0$. As $\tilde{\alpha} > 0$, we can see from (C3) and Fig. 18 that the potential $U(y)$ first decreases as $-\tilde{\alpha}y^2$ at low $y$. This corresponds to the fact that the particle coming from the right must take an infinite time ($x \to \infty$) to reach the zero-density point $y = 0$, by slowly climbing upward the potential $U(y)$. Note that the friction becomes negligible at late times because of the factor $2/x$. Thanks to the attractive self-interaction term $y^4/2$, the potential $U(y)$ turns upward to positive values at larger $y$, and next turns downward because of the large-density repulsive self-interaction term $-\frac{c_2}{3c_1}y^6$. In particular, there is only one local minimum $y_-$ and a global maximum $y_+ > y_-$ over the range $0 \le y < +\infty$. The particle must start slightly to the left of the maximum $y_+$ to roll down to $y = 0$, which is reached at infinite time. Therefore, we can see that the density profile can only reach large masses by having the particle start very close to the maximum $y_+$, so that it stays there for a very long time of the order of $x_+$, until it rolls down the potential $U(y)$ to finally settle at $y = 0$. This
<figure><img src="image_25.png" /><figcaption>FIG. 19. Trajectories in the potentials $U(y)$ for different soliton masses, corresponding to the soliton radial profiles of Fig. 11.</figcaption></figure>

convergence of the starting point to the maximum $y_+$ is clearly seen in Fig. 18 as we increase the soliton mass. This in turns means that we have a constant density core at $\rho = \rho_{\Lambda}y_0^2$ up to an increasingly large core radius $x_+$, beyond which the density falls off to converge to its exponential tail. This agrees with the profiles found in Fig. 2. The position of the maximum $y_+$ is set by the balance between the attractive and repulsive self-interactions, at $y^2 \sim \frac{c_1}{c_2} \sim 1$, and it does not significantly depend on the parameter $\tilde{\alpha}$, and hence nor on the soliton mass $M$. This means that the core density remains of the order of $\rho_{\Lambda}$ and stabilizes to a finite value for large masses, in agreement with Fig. 2.

## 2. Cosine model

For the cosine scalar-field potential (112), the soliton radial density profiles are given by Eq. (164). As in Eq. (C1), this can be interpreted as the trajectory of a particle $y(x)$ over time $x$, rolling down a potential $U(y)$ with a time-dependent friction term $\frac{2}{x} \frac{dy}{dx}$. The potential $U(y)$ is now given by

$$U(y) = J_0(y) - 1 - \frac{\tilde{\alpha}}{2} y^2. \qquad (C4)$$

We show in Fig. 19 the trajectories in the potentials $U(y)$ associated with the density profiles of Fig. 11. The parameter $\tilde{\alpha}$ is now negative, in the range $-\frac{1}{2} < \tilde{\alpha} < 0$. The density shows an exponential tail at large radii, with $y \sim e^{-\sqrt{\tilde{\alpha}+1/2}x}$. This again gives a potential $U(y)$ that decreases with $y$ near the origin, as $U(y) \sim -(1+2\tilde{\alpha})y^2/4$, so that the particle climbs the potential upward at late time to reach the origin at infinite time. In contrast with the polynomial case studied in Fig. 18, the soliton no longer gains mass by starting increasingly close to a maximum $y_+$ but by starting at increasingly large values $y_0$. There, the Bessel function $J_0(y)$ and its oscillations are negligible and the potential is dominated by the quadratic term $-\tilde{\alpha}y^2/2$. Thus, for large masses, the particle slowly rolls down the quadratic potential

$-\tilde{\alpha}y^2/2$ from increasingly large initial values $y_0$, until $y$ becomes of order unity with $U \simeq -1$. It next takes an infinite time to climb up to $(y=0, U=0)$. This agrees with the radial density profiles of Fig. 11, with a core density and a core radius that grow with the soliton mass.

We can also see that the oscillations of the Bessel function $J_0(y)$ in $U(y)$, which arise from the Bessel function $J_1(\sqrt{\rho/\rho_b})$ in the self-interaction potential $\Phi_I(\rho)$ in Eq. (115), have a negligible impact at large mass. This explains why the soliton mass-density relation shown in Fig. 10 takes the power-law form (162) predicted by the simple analytic ansatz (158)–(159). The Bessel function is not scale free but it only brings negligible deviations from the power-law behaviors arising from the main quadratic component $U(y) \simeq -\tilde{\alpha}y^2/2$.

We can infer that a similar behavior will be found for other models where the self-interaction component $U_I(y)$ goes to a constant or grows more slowly than $y^2$ at large $y$, in contrast with the polynomial case of Sec. C1 where we had $U_I(y) \sim -\frac{c_2}{3c_1}y^6$ at large $y$. In the general case, $U(y)$ is related to the self-interaction potential $V_I$ by

$$U(y) = \frac{\alpha}{2} y^2 - \frac{1}{2\rho_{\Lambda}} V_I(\rho_{\Lambda} y^2), \qquad (C5)$$

where we defined $y$ and the dimensionless radius $x$ by

$$y = \sqrt{\frac{\rho}{\rho_{\Lambda}}}, \qquad x = \sqrt{2mr}, \qquad (C6)$$

and $\rho_{\Lambda}$ is a characteristic density. Thus, models where $V_I(\rho)$ is bounded or grows more slowly than $\rho$, will have solitonic density profiles with a core density that grows at high masses, as in Figs. 19 and 11. In contrast, models where $V_I(\rho)$ has a minimum at a characteristic density $\rho_{\Lambda}$ and grows faster than $\rho$ at larger densities will have solitonic density profiles with a core density that converges to a finite value of order $\rho_{\Lambda}$ and a radius that grows as $M^{1/3}$ at high masses, as in Figs. 18 and 2.
[1] L. Bergström, in *40th Saas-Fee Course: Astrophysics at Very-High Energies*, Saas-Fee Advanced Course (Springer, 2013), pp. 123–222.

[2] K. Kowalska and E. M. Sessolo, *Adv. High Energy Phys.* **2018**, *6828560* (2018).

[3] M. Dine, W. Fischler, and M. Srednicki, *Phys. Lett. 104B*, **199** (1981).

[4] L. Abbott and P. Sikivie, *Phys. Lett. 120B*, **133** (1983).

[5] J. Preskill, M. B. Wise, and F. Wilczek, *Phys. Lett. 120B*, **127** (1983).

[6] R. D. Peccei and H. R. Quinn, *Phys. Rev. Lett.* **38**, 1440 (1977).

[7] F. Wilczek, *Phys. Rev. Lett.* **40**, 279 (1978).

[8] S. Weinberg, *Phys. Rev. Lett.* **40**, 223 (1978).

[9] M. Vysotsky, Y. Zeldovich, M. Khlopov, and V. Chechetkin, *Pis'ma Zh. Eksp. Teor. Fiz.* **27**, 533 (1978).

[10] D. J. E. Marsh, *Phys. Rep.* **643**, 1 (2016).

[11] L. A. Ureña-López, *Front. Astron. Space Sci.* **6**, 47 (2019).

[12] W. Hu, R. Barkana, and A. Gruzinov, *Phys. Rev. Lett.* **85**, 1158 (2000).

[13] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, *Phys. Rev. D* **95**, 043541 (2017).

[14] A. Sakharov and M. Khlopov, *Phys. At. Nucl.* **57**, 485 (1994).

[15] A. Sakharov, D. Sokoloff, and M. Khlopov, *Phys. At. Nucl.* **59**, 1005 (1996).

[16] M. C. Johnson and M. Kamionkowski, *Phys. Rev. D* **78**, 063010 (2008).

[17] J.-c. Hwang and H. Noh, *Phys. Lett. B* **680**, 1 (2009).

[18] C.-G. Park, J.-c. Hwang, and H. Noh, *Phys. Rev. D* **86**, 083535 (2012).

[19] R. Hlozek, D. Grin, D. J. E. Marsh, and P. G. Ferreira, *Phys. Rev. D* **91**, 103512 (2015).

[20] J. A. R. Cembranos, A. L. Maroto, and S. J. Núñez Jareño, *J. High Energy Phys.* **03** (2016) 013.

[21] J. A. R. Cembranos, A. L. Maroto, and S. J. Núñez Jareño, *J. High Energy Phys.* **02** (2017) 064.

[22] H.-Y. Schive, T. Chiueh, and T. Broadhurst, *Nat. Phys.* **10**, 496 (2014).

[23] T. Broadhurst, H. N. Luu, and S. H. H. Tye, *Phys. Dark Universe* **30**, 100636 (2020).

[24] J. P. Ostriker and P. J. Steinhardt, *Science* **300**, 1909 (2003).

[25] J. A. R. Cembranos, J. L. Feng, A. Rajaraman, and F. Takayama, *Phys. Rev. Lett.* **95**, 181301 (2005).

[26] D. H. Weinberg, J. S. Bullock, F. Governato, R. Kuzio de Naray, and A. H. G. Peter, *Proc. Natl. Acad. Sci. U.S.A.* **112**, 12249 (2015); **112**, 12249 (2015).

[27] A. Pontzen and F. Governato, *Nature (London)* **506**, 171 (2014).

[28] M. Boylan-Kolchin, J. S. Bullock, and M. Kaplinghat, *Mon. Not. R. Astron. Soc.* **415**, L40 (2011).

[29] B. Moore, S. Ghigna, F. Governato, G. Lake, T. R. Quinn, J. Stadel, and P. Tozzi, *Astrophys. J.* **524**, L19 (1999).

[30] W. J. G. de Blok, *Adv. Astron.* **2010**, 789293 (2010).

[31] J. A. R. Cembranos, A. L. Maroto, S. J. Núñez Jareño, and H. Villarrubia-Rojo, *J. High Energy Phys.* **08** (2018) 073.

[32] E. Armengaud, N. Palanque-Delabrouille, C. Yèche, D. J. E. Marsh, and J. Baur, *Mon. Not. R. Astron. Soc.* **471**, 4606 (2017).

[33] P. Brax, J. A. R. Cembranos, and P. Valageas, *Phys. Rev. D* **100**, 023526 (2019).

[34] P. Brax, J. A. Cembranos, and P. Valageas, *Phys. Rev. D* **101**, 023521 (2020).

[35] P. Brax, J. A. Cembranos, and P. Valageas, *Phys. Rev. D* **101**, 063510 (2020).

[36] P.-H. Chavanis, *Astron. Astrophys.* **537**, A127 (2012).

[37] B. Carr and F. Kuhnel, *arXiv:2006.02838*.

[38] P. Tisserand, L. Le Guillou, C. Afonso, J. N. Albert, J. Andersen, R. Ansari, É. Aubourg, P. Bareyre, J. P. Beaulieu, X. Charlot *et al.*, *Astron. Astrophys.* **469**, 387 (2007).

[39] P.-H. Chavanis, *Phys. Rev. D* **98**, 023009 (2018).

[40] A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, *Phys. Rev. D* **101**, 083014 (2020).

[41] E. W. Kolb and I. I. Tkachev, *Phys. Rev. Lett.* **71**, 3051 (1993).

[42] E. D. Schiappacasse and M. P. Hertzberg, *J. Cosmol. Astropart. Phys.* **01** (2018) 037; 03 (2018) E01.

[43] M. A. Amin, R. Easther, and H. Finkel, *J. Cosmol. Astropart. Phys.* **12** (2010) 001.

[44] M. A. Amin, R. Easther, H. Finkel, R. Flauger, and M. P. Hertzberg, *Phys. Rev. Lett.* **108**, 241302 (2012).

[45] J. Ollé, O. Pujolàs, and F. Rompineve, *J. Cosmol. Astropart. Phys.* **02** (2020) 006.

[46] H.-Y. Zhang, M. A. Amin, E. J. Copeland, P. M. Saffin, and K. D. Lozanov, *J. Cosmol. Astropart. Phys.* **07** (2020) 055.

[47] M. A. Amin and P. Mocz, *Phys. Rev. D* **100**, 063507 (2019).

[48] A. H. Guth, M. P. Hertzberg, and C. Prescod-Weinstein, *Phys. Rev. D* **92**, 103513 (2015).

[49] E. Madelung, *Z. Phys.* **40**, 322 (1927).

[50] P.-H. Chavanis, *Phys. Rev. D* **84**, 043531 (2011).

[51] D. S. Gorbunov and V. A. Rubakov, *Introduction to the Theory of the Early Universe: Cosmological Perturbations and Inflationary Theory* (World Scientific, 2011).

[52] L. Visinelli, S. Baum, J. Redondo, K. Freese, and F. Wilczek, *Phys. Lett. B* **777**, 64 (2018).

[53] F. Guzmán and A. A. Avilez, *Phys. Rev. D* **97**, 116003 (2018).

[54] B. Schwabe, J. C. Niemeyer, and J. F. Engels, *Phys. Rev. D* **94**, 043513 (2016).

[55] E. Cotner, *Phys. Rev. D* **94**, 063503 (2016).

[56] M. P. Hertzberg, Y. Li, and E. D. Schiappacasse, *J. Cosmol. Astropart. Phys.* **07** (2020) 067.

[57] H. Niikura *et al.*, *Nat. Astron.* **3**, 524 (2019).

[58] S. Sugiyama, T. Kurita, and M. Takada, *Mon. Not. R. Astron. Soc.* **493**, 3632 (2020).

[59] N. Smyth, S. Profumo, S. English, T. Jeltema, K. McKinnon, and P. Guhathakurta, *Phys. Rev. D* **101**, 063005 (2020).

[60] P. Schneider, J. Ehlers, and E. E. Falco, *Gravitational Lenses* (Springer, Berlin, Heidelberg, 1992).

[61] M. Bartelmann, *Classical Quantum Gravity* **27**, 233001 (2010).

[62] L. McAllister, E. Silverstein, A. Westphal, and T. Wrase, *J. High Energy Phys.* **09** (2014) 123.

[63] J. Berges, A. Chatrchyan, and J. Jaeckel, *J. Cosmol. Astropart. Phys.* **08** (2019) 020.

[64] M. Abramowitz, *Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables* (Dover Publications, Inc., New York, 1974), ISBN 0486612724.
[65] N. W. McLachlan, *Theory and Application of Mathieu Functions* (Oxford University Press, New York, 1947).

[66] C. Bender and S. Orszag, *Advanced Mathematical Methods for Scientists and Engineers* (Springer-Verlag, New York, 1999).

[67] J. Avron and B. Simon, *Ann. Phys. (N.Y.)* **134**, 76 (1981).

[68] B. Anahtarci and P. Djakov, *J. Math. Anal. Appl.* **396**, 243 (2012).

[69] H. Fukunaga, N. Kitajima, and Y. Urakawa, *J. Cosmol. Astropart. Phys.* **06** (2019) 055.

[70] A. Chatrchyan and J. Jaeckel, arXiv:2004.07844.

[71] J. Jaeckel, S. Schenk, and M. Spannowsky, arXiv:2004.13724.