#  Oscillons of Axion-Like Particle:
Mass distribution and power spectrum

Masahiro Kawasaki  Wakutaka Nakano  Hiromasa Nakatsuka  and Eisuke Sonomoto

###### Abstract

In string theory, the simultaneous existence of many Axion-Like Particles
(ALPs) are suggested over a vast mass range, and a variety of potentials have
been developed in the context of inflation. In such potentials shallower than
quadratic, the prominent instability can produce localized dense objects,
oscillons. Because of the approximate conservation of their adiabatic
invariant, oscillons generally survive quite long, maybe up to the current age
of the universe in the case of ultra-light ALPs with $m\sim 10^{-22}\ {\rm
eV}$. Such oscillons can have significant effects on the evolution of the
recent universe.

In this paper, we investigate the oscillons of the pure-natural type potential
by classical lattice simulation to explore the key quantities necessary for
phenomenological application: the number density of oscillons, the oscillon
mass distribution, the energy ratio of oscillons to the ALP field, and the
power spectrum. Then, we evolve these values in consideration of the analytic
decay rate.

##  1 Introduction

Although the $\Lambda$-CDM cosmology excellently explains the various
observations, some exceptions are found in small scale physics, so-called
”small scale crisis” [1] such as the core-cusp problem [2, 3, 4, 5]. One of
the solutions to the small scale crisis is the fuzzy dark matter, the ultra-
light bosonic dark matter with mass $\sim 10^{-22}\ \mathrm{eV}$. Such
extremely light mass is naturally explained by the pseudo-scalar field whose
mass is prohibited by a shift symmetry. Since the string theory predicts a
large number of axion-like particles (ALPs) with a broad range of masses, the
string axion is one of the promising candidates of the ultralight bosonic
fields.

In the literature, the ALP potential is often simply assumed as the cosine
type in analogy with QCD axion. However, other types of potentials are
suggested in the context of inflation [6, 7, 8, 9, 10, 11, 12, 13]. These
potential are shallower than quadratic, which results in the prominent
instability when ALP oscillates around the potential minimum. Such instability
can produce the gravitational waves in the reach of future observations [14,
15].

The large instability also produces the localized soliton of ALP, oscillon
[16, 17, 18]. For the QCD axion, such soliton is sometimes called ”axiton”
[19, 20, 21]. The formation of oscillons is confirmed on pure-natural type
potential numerically [22, 23, 24]. When the sizable energy fraction of the
ALPs is confined inside the oscillons, the constraints on the homogeneous ALPs
[25, 26, 27] should be reconsidered. The lifetime of oscillons is intensively
studied [28, 29, 30] in the last few years, and it is found that oscillons
live quite long. The long lifetime opens new possibilities of oscillons
discussed in several literature [28, 31, 32, 33]. However, no research has
been done so far on the current abundance of ALP oscillons including the
effects of their long lifetime and decay.

In this paper, we study properties of oscillons produced from ALP by lattice
simulation. First, we calculate the mass distribution, the energy ratio, the
number density, and the average mass of oscillons after the formation. Using
these values, we compare the simulated density power spectrum to the analytic
formula extended from the Poisson distribution. Finally, we study their time
evolution by estimating the analytic decay rate of oscillons.

Sec. 2 explains the nature of oscillons in the pure natural type potential.
Sec. 3 describes the setup of lattice simulations, which derive the number
density of oscillons, the energy ratio of oscillons to ALP, initial field
dependence, mass distribution, and the power spectrum in Sec. 4. Sec. 5
discusses the decay of oscillons and our paper concludes in Sec. 6.

##  2 Oscillons of Axion-Like Particle

Various kinds of ALP potentials are derived when the instanton effect breaks
the shift symmetry of ALP. In this paper, we use the pure natural type
potential arising from the coupling between the ALP and the pure Yang-Mills
gauge fields [12],

$\displaystyle
V(\phi)=\frac{m^{2}F^{2}}{2p}\left[1-\left(1+\frac{\phi^{2}}{F^{2}}\right)^{-p}\right].$
(2.1)

This potential is flatter than quadratic for $p>-1$, which causes the
instability of the field and produces localized objects, oscillons [24, 23].

We focus on the ALP which starts to oscillate during the radiation dominated
era. In this potential, the relation between the Hubble parameter at the onset
of the oscillation and the initial field value $\phi_{i}$ is represented by
[34]

$\displaystyle\left(\frac{H}{m}\right)^{2}\simeq\frac{1}{m^{2}\phi}\frac{\partial
V(\phi)}{\partial\phi}~{}\bigg{|}_{\phi=\phi_{i}}=\left(1+\frac{\phi_{i}^{2}}{F^{2}}\right)^{-(p+1)}.$
(2.2)

Thus, the large $|\phi_{i}|$ and large $p$ delay the start of the oscillation,
and the ALP mass must be

$\displaystyle m\gtrsim 10^{-28}\ {\rm eV}\
\left(1+\frac{\phi_{i}^{2}}{F^{2}}\right)^{\frac{p+1}{2}}\left(\frac{H_{\rm
eq}}{10^{-28}\ {\rm eV}}\right).$ (2.3)

where $H_{\rm eq}$ is the Hubble parameter at the matter radiation equality.

###  2.1 Instability Analysis

In this subsection, we will show that the coherent oscillation in this
potential causes two kinds of instabilities [15], which affect oscillon
formation as we will see in Sec. 4.2.

We perform the linear analysis by decomposing the field into the background
field $\phi_{0}(t)$ and its perturbation $\delta\phi(x)$ as

$\displaystyle\phi(t,\bm{x})=\phi_{0}(t)+\delta\phi(t,\bm{x}).$ (2.4)

We assume that the background field coherently oscillates as
$\phi_{0}(t)=\phi_{i}(a_{i}/a)^{3/2}\cos mt$ where the lower index $i$
represents the beginning of the oscillation, and define the Fourier components
as $\phi_{k}$ where $k$ is the comoving wavenumber.

The first kind of instability is the parametric resonance which causes
oscillon formation. When $\phi\ll F$, we can approximate the potential as
$V_{,\phi\phi}/m^{2}\simeq 1+\lambda(\phi/F)^{2}/2$ with $\lambda=-6(1+p)$. In
this case, the time evolution of $\phi_{k}$ is described by the Mathieu
equation,

$\displaystyle\frac{d^{2}}{d(mt)^{2}}\phi_{k}$
$\displaystyle+[A_{k}-2q\cos(2mt)]\phi_{k}=0,$ (2.5) $\displaystyle q$
$\displaystyle=-\frac{\lambda}{8}\left(\frac{\phi_{i}}{F}\right)^{2}\left(\frac{a_{i}}{a}\right)^{3},$
(2.6) $\displaystyle A_{k}$
$\displaystyle=1+\left(\frac{k}{am}\right)^{2}-2q.$ (2.7)

Neglecting the cosmic expansion ($a=a_{i}$), the most efficient resonance band
is

$\displaystyle\frac{k}{a_{i}m}\sim\sqrt{2q}=\frac{\sqrt{|\lambda|}|\phi_{i}|}{2F}.$
(2.8)

The oscillon formation becomes effective when the perturbations with the
wavelength of the oscillon size are amplified by the parametric resonance. If
these instabilities amplify the perturbations up to $\mathcal{O}(1)$, the
localized field configuration, the oscillon, can be formed,

The second kind of instability is tachyonic instability, which hinders the
coherent oscillation needed for the oscillon formation. The tachyonic
instability is caused by the negative effective mass,

$\displaystyle
V,_{\phi\phi}<0\quad\Leftrightarrow\quad\frac{\phi}{F}>\frac{1}{\sqrt{1+2p}}\equiv\frac{\phi_{c}}{F}.$
(2.9)

Thus, the potential has the inflection point only when $p>-1/2$. The tachyonic
instability occurs in the small $k$ mode with

$\displaystyle\frac{k}{am}<\left|\frac{V_{,\phi\phi}}{m^{2}}\right|^{1/2}.$
(2.10)

###  2.2 Oscillon Profile

The oscillon configuration is quasi-stable because of the approximate
conservation of the adiabatic invariant [35, 36],

$\displaystyle
I\equiv\frac{1}{\omega}\int\text{d}^{3}\bm{x}\overline{\dot{\phi}^{2}},$
(2.11)

where $\omega$ is the frequency of the periodic motion and the overline
represents the time average over one period. When we assume that $I$ is
invariant and set $\phi(t,\bm{x})=\Phi(r)\cos(\omega t)$, the lowest energy
configuration of the scalar field $\Phi(r)$ is determined by the following
differential equation:

$\displaystyle\frac{\text{d}^{2}\Phi}{\text{d}r^{2}}+\frac{2}{r}\frac{\text{d}\Phi}{\text{d}r}+\omega^{2}\Phi-2\,\overline{\frac{\text{d}V}{\text{d}\Phi}(\Phi\cos(\omega
t))\cos(\omega t)}=0,$ (2.12)

with boundary conditions $\Phi|_{r\to\infty}=0$ and
$\text{d}\Phi/\text{d}r|_{r\to 0}=0$. The profile of produced oscillons
depends on the parameter $\omega$. When we set $\omega$ smaller, the initial
condition of $\Phi(r=0)$ becomes larger, i.e., we obtain a larger oscillon.
The oscillon is produced under the condition that the Eq. (2.12) has the
solution, and the potential Eq. (2.1) satisfies this condition only when
$p>-1$ as we mentioned earlier.

The range of $\omega$ is constrained by the stability condition

$\displaystyle\frac{dI}{d\omega}<0.$ (2.13)

This charge $I$ relates to the classical particle number conservation [37, 38,
39, 40].

###  2.3 Oscillon Decay

Because the adiabatic invariance is only approximate, oscillons gradually
decay by emitting the self-radiation. In this subsection, we calculate the
classical decay rate of the oscillon by perturbative method following refs.
[29, 30].

The classical decay rate is defined by

$\displaystyle\Gamma(M_{\text{osc}})\equiv\frac{|\overline{\dot{M}_{\text{osc}}}|}{M_{\text{osc}}},$
(2.14)

where $M_{\text{osc}}$ is the energy of the oscillon configuration. We take
the radiation term $\xi(t,r)$ as

$\displaystyle\phi(t,{\bf x})=\Phi(r)\cos(\omega t)+\xi(t,r).$ (2.15)

The linear approximated equation of motion for $\xi(t,r)$ is

$\displaystyle\left[\,\square+V^{\prime\prime}(\Phi\cos(\omega
t))\right]\xi=-V^{\prime}(\Phi\cos(\omega
t))+2\,\overline{V^{\prime}(\Phi\cos(\omega t))\cos(\omega t)}\cos(\omega t),$
(2.16)

where prime denotes the derivative of $\phi$.

It is very convenient to use the Fourier expansions in time for the
calculation of the decay rate. We expand
$\xi(t,r)=\sum_{j=3}\xi_{j}(r)\cos(j\omega t)$ and the potential terms. The
equation of motion for $\xi_{j}$ is written as

$\displaystyle\left[\nabla^{2}+\kappa_{j}^{2}\right]\xi_{j}=-\mathcal{S}_{j},$
(2.17)

where $\kappa_{j}=\sqrt{j^{2}\omega^{2}-m^{2}}$ and $\mathcal{S}_{j}(r)$ is
the source term which includes the combination of $\xi_{k}$’s. The solution of
$\xi_{j}$’s is derived under the same boundary condition as for the oscillon
profile calculation, but it includes not only inhomogeneous solution but also
homogeneous solution, i.e. free field solution. We compare the derived
solutions with the formal outgoing solutions at the origin to remove the
homogeneous solution,

$\displaystyle\xi_{j}(0)=\int_{0}^{\infty}dr^{\prime}\mathcal{S}_{j}(r^{\prime})\,r^{\prime}\cos(\kappa_{j}r^{\prime}).$
(2.18)

The radiative solution of $\xi$ for $r\to\infty$ is obtained as

$\displaystyle\xi_{j}^{\mathrm{rad}}(t,r)\approx$ $\displaystyle\frac{1}{4\pi
r}\sum_{j=3}\tilde{\mathcal{S}}_{j}(\kappa_{j})\cos(\kappa_{j}r-j\omega t),$
(2.19) $\displaystyle\tilde{\mathcal{S}}_{j}(k)=$ $\displaystyle\int d^{3}{\bf
r}^{\prime}\mathcal{S}_{j}(r^{\prime})\,e^{-i{\bf k\cdot r}^{\prime}}.$ (2.20)

The time averaged classical energy loss rate is given by

$\displaystyle\overline{\dot{M}_{\text{osc}}}=4\pi
r^{2}\,\overline{\partial_{t}\xi^{\mathrm{rad}}(t,r)\partial_{r}\xi^{\mathrm{rad}}(t,r)}=-\frac{1}{8\pi}\sum_{j=3}j\omega\kappa_{j}\tilde{\mathcal{S}}_{j}(\kappa_{j})^{2}.$
(2.21)

Therefore, we can write $\Gamma=\sum_{j=3}\Gamma_{j}$ and find $\Gamma_{j}\geq
0$. In the following, we just calculate $\Gamma_{3}$ and assume
$\Gamma\simeq\Gamma_{3}$ because there is relation of
$|\xi_{j}|\gg|\xi_{j+1}|$ (hence $\Gamma_{j}\gg\Gamma_{j+1}$) confirmed in
$\omega\gtrsim 0.8m$ except for poles of $\Gamma_{3}$ according to the
previous study [30].

The decay rates are shown in Fig. 1 as a function of $\omega$. The upper limit
of $\omega$ is determined by the relation $dI/d\omega=0$. In Fig. 1, there are
poles that the decay rate is extremely small. Generally, the position of poles
of the decay rate $\Gamma_{j}$ does not correspond to that of $\Gamma_{j+1}$
and therefore, the relation $\Gamma_{j+1}\gg\Gamma_{j}$ could hold at the
poles of $\Gamma_{j}$. This means that the lifetime at poles calculated here
is the maximal lifetime, i.e., the lifetime could be overestimated by the pole
effect.

Figure 1:  Decay rates of oscillons calculated based on Eqs. (2.14) and (2.21)
with $j=3$. The oscillon mass is related to $\omega$. The different lines
shows the different $p$.

The time dependence of oscillon mass is formally given by

$\displaystyle
M_{\text{osc}}(t)=M_{\text{osc}}(0)\exp\left(-\int^{t}_{0}\Gamma(M_{\text{osc}}(t^{\prime}))\text{d}t^{\prime}\right).$
(2.22)

We will discuss the time evolution of oscillons based on Eq. (2.22) in Sec. 5.

###  2.4 Extended Poisson Distribution

The dark matter oscillons produce large isocurvature fluctuations. The simple
approximation assuming the random distribution of point particles with a
single mass leads to Poisson distribution. In this section, extending it with
consideration of broad mass of oscillons and energy conservation inside the
horizon scale, we analytically derive the power spectrum of ALP oscillons. In
the following, we use the comoving coordinate $\bm{x}$ and momentum $\bm{k}$.

Assuming that the spatial profile of oscillons is approximated by the Gaussian
function and their velocity is almost zero due to the Hubble friction, the
energy density of a single oscillon $\rho_{\rm osc}$ with the total energy $M$
is represented by

$\displaystyle\rho_{\text{osc}}(t,\bm{x};\bm{x}_{0},M)$
$\displaystyle\simeq\frac{[\omega(M)\psi_{0}(M)]^{2}}{2}\exp\left(-\frac{a(t)^{2}(\bm{x}-\bm{x}_{0})^{2}}{R(M)^{2}}\right)$
(2.23)

where $\bm{x}_{0}$, $R(M)$, $\omega(M)$, and $\psi_{0}(M)$ are the comoving
position, the physical radius, the oscillation frequency, and the central
field value of the oscillon. Then, the total energy of the oscillon is written
as

$\displaystyle
M=a(t)^{3}\int\text{d}^{3}\bm{x}\rho_{\text{osc}}(t,\bm{x};\bm{x}_{0},M)=\frac{[\omega(M)\psi_{0}(M)]^{2}}{2}(\pi
R(M)^{2})^{3/2}.$ (2.24)

The density perturbation $\delta_{\text{tot}}$ is given by the contribution of
$N$ oscillons over the average energy density of ALP ${\rho}_{\phi}$ including
oscillons, i.e,

$\displaystyle\delta_{\rm tot}(\bm{x})$
$\displaystyle\equiv\sum_{i}^{N}\frac{\rho_{\rm
osc}(t,\bm{x};\bm{x}_{i},M_{i})}{\rho_{\phi}(t)},$ (2.25)

where and $\bm{x}_{i}$ and $M_{i}$ is the comoving coordinate and mass of the
$i$-th oscillon. Then, the Fourier component of $\delta_{\rm tot}$ is

$\displaystyle\delta_{\bm{k}}(t)$
$\displaystyle\equiv\int\text{d}\bm{x}^{3}\delta_{\rm
tot}(\bm{x})e^{-i\bm{x}\cdot\bm{k}},$ (2.26)
$\displaystyle=\frac{1}{a(t)^{3}\rho_{\phi}(t)}\sum_{i}^{N}e^{-i\bm{x}_{i}\cdot\bm{k}}M_{i}\exp\left(-\frac{R(M_{i})^{2}k^{2}}{4a(t)^{2}}\right),$
(2.27)

To calculate the sum $\sum_{i}^{N}$, let us consider the comoving volume $V$
including $N$ oscillons. When the volume $V$ is smaller than the horizon size
at oscillon formation, we expect that the oscillons are randomly distributed.
The power spectrum is given by the ensemble average of the oscillon
coordinates, 111 Here, the power spectrum is defined as
$\displaystyle\braket{\delta_{\bm{k}}\delta_{\bm{k}^{\prime}}}$
$\displaystyle=(2\pi)^{3}\delta^{3}(\bm{k}-\bm{k}^{\prime})P_{\text{osc},k}\quad,\quad\delta_{\bm{k}}\equiv\int\text{d}x^{3}\delta_{\rm
tot}e^{-i\bm{x}\cdot\bm{k}}.$ (2.28)

$\displaystyle P_{\text{osc},k}(t)$
$\displaystyle\equiv\left(\prod_{i}\int_{V}\frac{\text{d}\bm{x}_{i}^{3}}{V}\right)\frac{|\delta_{\bm{k}}(t)|^{2}}{V},$
$\displaystyle=\left(\frac{1}{a(t)^{3}\rho_{\phi}(t)}\right)^{2}\left(\prod_{i}\int_{V}\frac{\text{d}\bm{x}_{i}^{3}}{V}\right)\frac{1}{V}\sum_{l,n}^{N}e^{i(\bm{x}_{l}-\bm{x}_{n})\cdot\bm{k}}M_{l}M_{n}\exp\left(-\frac{(R(M_{l})^{2}+R(M_{n})^{2})k^{2}}{4a(t)^{2}}\right)$
$\displaystyle\simeq\left(\frac{1}{a(t)^{3}\rho_{\phi}(t)}\right)^{2}\frac{1}{V}\sum_{i}^{N}M_{i}^{2}\exp\left(-\frac{R(M_{i})^{2}k^{2}}{2a(t)^{2}}\right),$
$\displaystyle=\frac{r_{\text{osc}}^{2}}{n_{\text{osc}}(t)}\frac{\Braket{M^{2}\exp\left(-\frac{R(M)^{2}k^{2}}{2a(t)^{2}}\right)}}{\braket{M}^{2}}.$
(2.29)

where we define $\braket{}$ as the average over mass $M$, and $n_{\rm osc}$
and $r_{\rm osc}$ are number density and energy ratio of oscillons to the ALP,
respectively: 222 We have used physical energy density and comoving number
density.

$\displaystyle n_{\rm osc}\equiv\frac{N}{V}\quad,\quad
r_{\text{osc}}=\frac{n_{\text{osc}}(t)\braket{M}}{a(t)^{3}\rho_{\phi}(t)}\quad.$
(2.30)

On a scale larger than the horizon size at the oscillon formation, the
position of oscillons is not purely randomized because of the energy
conservation inside the horizon. More precisely, energy is conserved in a
volume inside which the energy transfer of the scalar field is not effective.
Thus, the power spectrum damps over the scale of energy transfer at the
oscillon formation $k_{c}$, which is smaller than the horizon size.
Considering the energy conservation, the average of oscillon coordinates in
Eq. (2.29) is modified, and we introduce the following damping factor (see
Appendix A for the derivation)

$\displaystyle P_{\text{osc},k}(t)$
$\displaystyle=\frac{r_{\text{osc}}^{2}}{n_{\text{osc}}(t)}\frac{\Braket{M^{2}\exp\left(-\frac{R(M)^{2}k^{2}}{2a(t)^{2}}\right)}}{\braket{M}^{2}}K(k/k_{c}),$
(2.31) $\displaystyle K(x)$
$\displaystyle\equiv\left[1-\left(\frac{2}{x}\right)^{2}\sin^{2}\left(\frac{x}{2}\right)\right].$
(2.32)

where $k_{c}$ is the cutoff scale of perturbations and should be determined by
simulation.

Note that the above formula reproduces the Poisson distribution when we
neglect the size of oscillons ($R(M)\to 0$), use the monochromatic mass
($\braket{M^{2}}=\braket{M}^{2}$) and ignore the energy conservation
($k_{c}/k\ll 1$).

We call the above formula as extended Poisson distribution in this paper. In
the following section, we will confirm that the above estimation is in good
agreement with the power spectrum derived from the numerical simulation.

##  3 Simulation Setup

In this section, we briefly explain the setup of the lattice simulations to
investigate the oscillon formation on various $p$’s and initial field values.
We use tilde to denote a simulation value normalized by a physical scale. The
normalized length and time in units of $m$ and the field value in the units of
$F$ are given by

$\displaystyle(\tilde{t},\tilde{\bm{x}})\equiv
m(t,\bm{x})\quad,\quad\tilde{\phi}\equiv\frac{\phi}{F}.$ (3.1)

In this normalization, the mass of oscillon $M$ is normalized as

$\displaystyle\tilde{M}=\frac{m}{F^{2}}M.$ (3.2)

###  3.1 Simulation Parameters

We take the following initial condition in simulations:

$\displaystyle\tilde{\phi}_{i}(\bm{x})$
$\displaystyle=\tilde{\phi}_{i}(1+\xi(\bm{x}))\quad,\quad\tilde{\phi}_{i}^{\prime}=0\quad,\quad
H_{i}=m^{-1}\quad,\quad\tilde{a}_{i}=1,$ (3.3)

where $\tilde{\phi}_{i}$ is order one, and $\xi(\bm{x})$ is the initial noise
defined by the scale-free power spectrum

$\displaystyle\braket{\xi_{\bm{k}}\xi_{\bm{k}^{\prime}}}=(2\pi)^{3}\delta^{3}(\bm{k}-\bm{k}^{\prime})\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{\xi},$
(3.4)

with a small constant $\mathcal{P}_{\xi}=2.1\times 10^{-9}$ as a reference
value333Note that this definition slightly changes in two dimension as
$\braket{\xi_{\bm{k}}\xi_{\bm{k}^{\prime}}}=(2\pi)^{2}\delta^{2}(\bm{k}-\bm{k}^{\prime})\frac{2\pi}{k^{2}}\mathcal{P}_{\xi}.$
. Because we assume the radiation-dominated era and use the conformal time
$\tilde{\tau}\equiv\sqrt{2\tilde{t}}=\tilde{a}$, the initial physical and
conformal times are given by $\tilde{t}_{i}=m/(2H_{i})=1/2$ and
$\tilde{\tau}_{i}=m\tau_{i}=1$.

The parameters of our simulations are summarized in Table 1 for three
dimensions and Table 2 for two dimensions. The small box setup ($L=4,8$) is
used to investigate the initial field dependence in Sec. 4.2, and the large
box setup ($L=16,32$) is used to investigate the oscillon mass distribution in
Sec. 4.3. Both setup satisfies the condition that all oscillons should be
resolved at the end of simulations. We use the threshold value $\rho_{\rm th}$
to identify the oscillons explained in the next subsection.

Table 1:  The setups of our three dimensional numerical simulations. We use the setup with the small box size when we need a large number of iterations and the large one when we calculate the oscillon mass distribution and the power spectrum.  $p$ | Box size | Grid size | Time step | Final time | Threshold
---|---|---|---|---|---
$L$ | $N$ | d$\tilde{\tau}$ | $\tilde{\tau}_{e}$ | $\tilde{\rho}_{\rm th}$
$-0.75$ | $8$ | $512^{3}$ | $4.0\times 10^{-3}$ | $81.0$ | $0.05$
$-0.5$ | $4$ | $256^{3}$ | $4.0\times 10^{-3}$ | $121.0$
$32$ | $1024^{3}$ | $8.0\times 10^{-3}$ | $81.0$
$0.5$ | $4$ | $512^{3}$ | $2.0\times 10^{-3}$ | $101.0$ | $0.04$
$2$ | $4$ | $256^{3}$ | $4.0\times 10^{-3}$ | $121.0$ | $0.01$
$16$ | $1024^{3}$ | $81.0$
$3$ | $4$ | $256^{3}$ | $121.0$
Table 2:  The setup of our two dimensional numerical simulations. We use this setup only for the calculation of the oscillon power spectrum.  $p$ | Box size | Grid size | Time step | Final time | Threshold
---|---|---|---|---|---
$L$ | $N$ | d$\tilde{\tau}$ | $\tilde{\tau}_{e}$ | $\tilde{\rho}_{\rm th}$
$-0.5$ | $1024$ | $32768^{2}$ | $8.0\times 10^{-3}$ | $81.0$ | 0.025

The lattice code adopts the fourth-order symplectic integration scheme for
time evolution and the fourth-order central differential scheme for spatial
derivatives with the periodic boundary condition [29, 41, 23]. We investigated
the effect of the boundary condition with different box sizes and confirmed
that the box size difference does not change the results.

###  3.2 Oscillon Identification

The numerical simulation of localized object are investigated in the early
works for Q-ball [42] and oscillon [43, 44, 45].

To estimate the contribution from the tail of the oscillon profile correctly,
we regard the connected regions with the energy density larger than the
threshold value $\tilde{\rho}_{\text{th}}$ and one more grid outside the
surface of them as oscillons. The threshold value $\tilde{\rho}_{\rm th}$ is
chosen so that $\tilde{\rho}_{\rm
th}=\alpha\tilde{\rho}_{\phi}(\tilde{\tau}_{\rm e})$ where $\alpha$ is
$\mathcal{O}(0.1)$ parameter and the time $\tilde{\tau}_{\rm e}$ is the end of
oscillon formation.

The mass of the oscillons is estimated by summing over the energy density
$\tilde{\rho}$ inside the identified region [42] as

$\displaystyle\tilde{M}_{\rm osc}=\sum_{\bm{x}\in
V_{\text{osc}}}\tilde{\rho}(\bm{x})\left(\frac{\tilde{a}L}{N}\right)^{d},$
(3.5)

where $d$ is the spatial dimension and the sum is over the grid $\bm{x}$
inside the oscillon volume $V_{\text{osc}}$. The number of oscillon inside the
box $N_{\rm osc}$ is given by the number of these disconnected regions.

We define the energy ratio of oscillons to the total scalar field as

$\displaystyle r_{\rm osc}(t)=\frac{\sum_{i}^{N_{\rm osc}}\tilde{M}_{{\rm
osc},i}}{\tilde{\rho}(t)(\tilde{a}(t)L)^{d}},$ (3.6)

and the normalized comoving number density of the oscillons $\tilde{n}_{\rm
osc}$ and the average mass of the oscillon $\langle\tilde{M}_{\rm osc}\rangle$
as

$\displaystyle\tilde{n}_{\rm
osc}=\frac{N_{\text{osc}}}{L^{d}},\quad\braket{\tilde{M}_{\text{osc}}}=\frac{\sum^{N_{\text{osc}}}_{i}\tilde{M}_{{\rm
osc},i}}{N_{\rm osc}}.$ (3.7)

Because of the contribution from the surface region, the threshold value and
the spatial resolution slightly change $M_{\rm osc}$ and hence $r_{\rm osc}$.
We found that the estimation of energy ratio defined in Eq. (3.6) differs less
than $5\%$ even if we change the spatial resolution and the threshold value up
to a factor of $2$.

##  4 Simulation Results

In this section, we exhibit the results of the classical lattice simulations.
In Sec. 4.1, we show the time dependence of the number of oscillons and
$r_{\text{osc}}$. In Sec. 4.2, we discuss the initial value dependence of
$r_{\rm osc}$ changing the initial field value $\tilde{\phi}_{i}$ with small
box setup for each $p$. In Sec. 4.3, we simulate the oscillon formation with
the large box setup to investigate the oscillon mass distribution. Finally in
Sec. 4.4, we numerically calculate the power spectrum of oscillons and compare
it with the analytic formula in Eq. (2.31).

###  4.1 Time Evolution of Oscillon Number and its Energy Ratio

---
Figure 2:  The time evolution of the number (left figure) and energy ratio
(right figure) of oscillons. The energy ratio converge to the stable value at
$\tilde{\tau}>40$, which we call the energy ratio of oscillons
$r_{\text{osc}}$.

First, we simulate the time evolution of the oscillon number $N_{\rm osc}$ and
the oscillon energy ratio $r_{\rm osc}$ changing the threshold value
$\tilde{\rho}_{\rm th}$ to observe their threshold dependence during the
oscillon formation. The result for $p=-0.5$ and $\tilde{\phi}_{i}=5.0\pi$ with
$L=4$ is shown in Fig. 2. The red, green, and blue lines represent
$\bar{\rho}_{\text{th}}=0.025,~{}0.05,~{}0.1$ respectively.

The left figure shows the number of oscillons. After the formation around
$\tilde{\tau}\sim 20$, $N_{\text{osc}}$ slightly decreases until
$\tilde{\tau}\sim 40$ because of the merger of the produced oscillons and the
damping of fluctuations larger than $\tilde{\rho}_{\rm th}$ by the gradient
term and the cosmic expansion. Then, $N_{\text{osc}}$ becomes stable
independent of the threshold values by the simulation end. In the following
calculation, we evaluate the energy of each oscillon when $N_{\text{osc}}$
becomes stable.

The right figure shows the time evolution of $r_{\text{osc}}(t)$. After the
oscillon number becomes stable, the energy ratio also becomes stable, which
hardly depends on the threshold values. Hereafter, we call this stable value
as the energy ratio of oscillon $r_{\text{osc}}$ 444 For some choice of
parameters, $r_{\text{osc}}(t)$ slightly decreases since the large oscillons
decay in the simulation time. We discuss the details of decay in the next
section. .

###  4.2 The Initial Field Value Dependence of Energy Ratio

Next, we examine the dependence of $r_{\text{osc}}$ on the initial field value
$\phi_{i}$. For each index $p$, we perform the simulations with the small box
setup in Table 1 changing initial field value $\phi_{i}$.

Figure 3:  The oscillon energy ratio with different index $p$ and initial
field value $\phi_{i}$. The different line corresponds to the different index
$p$. The error bar shows the variance over three times of simulations with
different random seeds.

We show $r_{\text{osc}}$ as a function of $\phi_{i}$ for each $p$ in Fig. 3.
For a small $\phi_{i}$, because the instability is too small for the oscillon
formation, $r_{\rm osc}$ ends up in zero. For a larger $\phi_{i}$, $r_{\rm
osc}$ also gets larger because the parametric resonance becomes effective.
Then, for $p=-0.75,-0.5$, $r_{\rm osc}$ converges on a plateau in Fig. 3, and
we define $\tilde{\phi}_{i,m}$ as the typical field value on the plateau. On
the other hand, when $p>-0.5$, the tachyonic instability amplifies the modes
different from the oscillon size, which hinders the subsequent oscillon
formation as discussed in Sec. 2.1. Therefore, the oscillon formation rate has
a peak at some value of $\phi_{i}$ for $p>-0.5$ and we define it as
$\tilde{\phi}_{i,m}$ 555 To confirm the effect of tachyonic instability, we
calculate $r_{\text{osc}}$ for $p=2$ with smaller initial noise
$\mathcal{P}_{\xi}$. We find that $r_{\text{osc}}$ increases for
$\phi_{i}=0.9\pi>\phi_{i,m}$ and decreases for $\phi_{i}=\phi_{i,m}$. For the
former case, the small initial noise delays the tachyonic instability. For the
latter case, $\phi_{i,m}$ is small enough to avoid the tachyonic instability
and the small noise just results in small oscillon formation. This result is
consistent with the effect of tachyonic instability. . For the later use, we
perform the same analysis in two dimensions with $p=-0.5$.

We comment on the initial noise dependence of the oscillon energy ratio.
Because the instability growth only logarithmically depends on the size of the
initial fluctuations, we expect that $r_{\text{osc}}$ does not depend on the
initial noise much. Actually, we confirmed that $r_{\text{osc}}$ typically
varies less than $10\%$ when we change initial noise $\mathcal{P}_{\xi}$ by a
factor of about $10$.

###  4.3 The Mass Distribution of Oscillons

Here, we fix the initial field value as $\tilde{\phi}_{i}=\tilde{\phi}_{i,m}$
derived in the previous subsection and perform simulations for $p=-0.5$ and
$2$ in the large box setup and $p=-0.75,~{}0.5$, and $~{}3$ in the small box
setup in Table 1. We have simulated in each setup three times with different
random seeds to reduce the stochastic noise.

The mass distributions of the comoving number density of oscillons for
$p=-0.5$ and $2$ at $\tilde{\tau}=81$ are shown in Fig. 4 by the blue
histograms. Note that the vertical axis shows the number of oscillons in each
bin normalized by $L^{3}$. The orange and green regions show the oscillon
distribution after the oscillon decay over $m\Delta t$ and will be discussed
in Sec. 5. Using these mass distributions, we calculate the comoving number
density of oscillon $\tilde{n}_{\text{osc}}$ and averaged mass
$\braket{\tilde{M}_{\text{osc}}}$ in Eq. (3.7) for each $p$, which are
summarized in Table 3.

---
Figure 4:  The histogram of oscillon mass with $p=-0.5$ (left figure) and $p=2$ (right figure). The horizontal axis shows the normalized oscillon mass $\tilde{M}_{\rm osc}$ and the vertical axis shows the comoving number density in each bin $\Delta\tilde{n}_{i}\equiv N_{i}/L^{3}$ where $N_{i}$ is the number of oscillons in the $i$-th bin. Each bar is binned by the log scale $[10^{0.01i},10^{0.01(i+1)}]$ and the error bar is given by $\sqrt{N_{i}}/L^{3}$. The blue region shows the oscillon distribution at $\tilde{\tau}=81$ and the orange and green regions show the oscillon distribution after the oscillon decay over $m\Delta t$ discussed in Sec. 5. The split form of $p=2$ with $m\Delta t=10^{5}$ comes from the pole structure of the decay rate, which has the uncertainty discussed in Sec. 2.3.  Table 3:  The results of simulations. $\tilde{\phi}_{i,m}$ is the initial field value for efficient oscillon formation discussed in Sec. 4.2. $r_{\text{osc}}$, $\tilde{n}_{\text{osc}}$ and $\braket{\tilde{M}_{\text{osc}}}$ are the oscillon energy ratio, number density, and averaged mass when $\tilde{\phi}_{i}=\tilde{\phi}_{i,m}$.  $p$ | $\tilde{\phi}_{i,m}$ | $r_{\text{osc}}$ | $\tilde{n}_{\text{osc}}$ | $\braket{\tilde{M}_{\text{osc}}}$
---|---|---|---|---
$-0.75$ | $12.0\pi$ | 0.73 | 0.98 | 2100
$-0.5$ | $5.0\pi$ | 0.69 | 0.83 | 570
$0.5$ | $1.5\pi$ | 0.59 | 0.78 | 130
$2$ | $0.7\pi$ | 0.60 | 0.79 | 47
$3$ | $0.55\pi$ | 0.62 | 0.85 | 32

###  4.4 Oscillon Power Spectrum

To validate the analytical estimation of the power spectrum Eq. (2.31), we
numerically calculate the power spectrum after oscillon formation. In this
subsection, we only show the results for $p=-0.5$ in the large box setup with
$\tilde{\phi}_{i}=5\pi$ in three dimensions and $\tilde{\phi}_{i}=2.5\pi$ in
two dimensions, but similar result is obtained for $p=-0.75$ and $2$.

Figure 5:  The normalized power spectrum with $p=-0.5$. The horizontal axis
represents the normalized comoving wavenumber and the vertical axis represents
the normalized power spectrum. We also plot the reproduced power spectrum
based on the extended Poisson distribution discussed in Sec. 2.4.

In Fig. 5, the blue and green lines show the snapshot of the power spectrum at
$\tilde{\tau}=81.0$ when the oscillon formation has already finished. The red
line shows the analytical power spectrum from Eq. (2.31) based on the oscillon
distribution in Fig. 4. When the cutoff scale $k_{c}$ in Eq. (2.31) is taken
as $k_{c}/(a(t_{m})m)=0.35$, we find that the numerical power spectrum (green
and blue lines) are well approximated by the analytical formula (red line).

The power spectrum consists of three parts. The center plateau region
corresponds to the Poisson distribution of the oscillons. On the scale smaller
than the oscillon size, the power spectrum decreases following the oscillon
profile given by Eq. (2.23). On the scale larger than $k_{c}^{-1}$, the power
spectrum decreases because of the energy conservation following Eq. (2.32).
This suppression, however, does not hold for the lower end (see blue line in
Fig. 5). One of the possibilities is that the oscillon energy ratio at each
horizon itself fluctuates, which results in another source of Poisson noise
with amplitude $\mathcal{O}(10^{2})$ smaller than the plateau region.

##  5 Time Evolution of Oscillons

As we introduced in Sec. 2, the decay rate of the oscillons is too small to
simulate numerically. Therefore, in this section, we include the decay process
of oscillons based on the analytic formula Eq. (2.22) to calculate the
oscillon energy ratio $r_{\rm osc}$, the average oscillon energy $\langle
M_{\rm osc}\rangle$, and the oscillon power spectrum.

First, we evolve the oscillon mass distribution derived in Sec. 4.3 by time
$m\Delta t$ and the results of $p=-0.5$ and $2$ are shown in Fig. 4. Since the
lighter oscillons are more stable for $p=-0.5$, oscillons converge on the
quasi-stable mass. On the other hand, because for $p=2$ the oscillon decay
rate has a pole at $M_{\text{osc}}\sim 50$ (see Fig. 1) where the decay rate
is extremely small, the oscillon distribution becomes spitted into two peaks
666 Note that the lifetime could be overestimated by poles for $p=0.5$, $2$
and $3$ as explained in Sec. 2.3. .

Based on this mass distribution, we calculate the time evolution of the
oscillon ratio $r_{\rm osc}$ and the average oscillon mass $\braket{M_{\rm
osc}}$ shown in Fig. 6. The time evolution of $\braket{M_{\rm osc}}$ is
different for each $p$. When $p=-0.75$ or $-0.5$, the average mass
monotonically decreases because the lightest oscillon is the most stable as
one can see from Fig. 1. On the other hand, when $p=0.5$, $2$ or $3$, the
lifetime of oscillons at the poles of the decay rate is the longest. Thus, the
total number of oscillons $N_{\rm osc}$ decreases by the decay of smaller
oscillons and it temporarily increases the average mass for $p=2$ and $3$ at
$\Delta\tilde{t}\sim 10^{5}$ which corresponds to the lifetime of the smallest
oscillon.

---
Figure 6:  The time dependence of the oscillon ratio and averaged mass when
oscillons decay. We use the typical initial condition summarized in Table. 3.
For $\Delta\tilde{t}\gtrsim 10^{5}$, the oscillons with $p=0.5,$ 2 and 3 are
stabilized by the small decay rate at the pole, which has the uncertainty
discussed in Sec. 2.3.

Based on these results, we consider the power spectrum after the partial decay
of oscillons. When the oscillon lifetime is long enough to survive until the
matter-dominated era, the power spectrum is gravitationally amplified. Then,
the realistic power spectrum during the matter-dominated era is represented by

$\displaystyle P_{\text{osc},k}(t)$
$\displaystyle=\left(\frac{3}{2}\frac{a(t)}{a_{\text{eq}}}\right)^{2}\frac{r_{\text{osc}}^{2}}{n_{\rm
osc}}\frac{\Braket{M^{2}}}{\braket{M}^{2}}\left[1-\left(\frac{2k_{c}}{k}\right)^{2}\sin^{2}\left(\frac{k}{2k_{c}}\right)\right].$
(5.1)

where $a_{\text{eq}}$ is the scale factor at the matter radiation equality.
Note that the factor $\Braket{M^{2}}/\braket{M}^{2}$ converges on 1 after the
partial decay of oscillons because the mass distribution becomes
monochromatic. The result is plotted in Fig. 7 for $p=-0.5$. Because the
gravitational growth only occurs on the scale which includes at least two
oscillons, we cut off the power spectrum on the scale which includes 10
oscillons $n_{\rm osc}(2\pi/k)^{3}>10$ for instance, shown as the dotted
lines.

Figure 7:  The power spectrum of oscillons for $p=-0.5$ based on Eq. (5.1).
The different lines describe the power spectrum at different time including
the effect of oscillon decay. We normalize the power spectrum so as to neglect
the growth during matter dominated era. The dotted line shows the scales which
includes less than 10 oscillons, $k>2\pi(n_{\text{osc}}/10)^{1/3}$.

Finally, we comment on the power spectrum after oscillons have completely
decayed. Suppose that the oscillon with radius $R$ and energy $M$ breaks down
radiating the ALPs with the kinetic energy $K\gtrsim m$ 777 $K\gtrsim m$ is
satisfied for the perturbative decay discussed in Sec. 2.3. However, it is
unclear whether this condition holds on the sudden decay at $dI/d\omega=0$,
which is out of the scope of this paper. , the escape velocity from the
oscillon surface $K_{\rm esc}\simeq GmM/R\sim m(F/M_{Pl})^{2}$ is smaller than
$K$. Thus, the decaying oscillon cannot gravitationally bind itself and the
Poisson fluctuation would damp after the oscillon evaporation.

##  6 Conclusion

In this paper, we have performed a comprehensive analysis of oscillons from
its formation to the decay and have established the basics of oscillon
phenomenology.

We firstly performed classical lattice simulation and calculated various
quantities: the dependence of energy ratio of oscillons to the ALP field on
the initial field value in Fig. 3, the oscillon mass distribution in Fig. 4,
and the power spectrum of oscillons in Fig. 5. We found that the oscillon
formation becomes efficient when the tachyonic instability which hinders the
instability mode necessary for oscillon formation is weak. Comparing the
resultant power spectrum and the analytical formula derived in Sec. 2.4, we
confirmed that the extended Poisson distribution is well approximation of the
oscillon power spectrum.

We analyzed the time evolution of oscillons in Sec. 5 including the effect of
their decay. With the simulation results, we calculated the current number
density and typical mass of the oscillons. Combining them with the oscillon
decay and the gravitational growth, we have estimated the power spectrum of
oscillons after oscillons partially decay shown in Fig. 5.

Finally, let us comment on the gravitational amplification of oscillon
fluctuations. Although we simply treated them as the matter fluctuations in
Sec. 5, the evolution of the fluctuations produced by oscillons in the matter-
dominated universe is still unclear. Thus, given our result of the matter
power spectrum as the initial condition, we need a more precise study
including the gravitational effect in future work.

## Acknowledgments

We would like to thank T. Hiramatsu for useful comments on numerical
calculations. This work was supported by JSPS KAKENHI Grant Nos. 17H01131
(M.K.), 17K05434 (M.K.), 19H05810 (W.N.), JP19J21974 (H.N.), 19J12936 (E.S.),
Advanced Leading Graduate Course for Photon Science (H.N.), and World Premier
International Research Center Initiative (WPI Initiative), MEXT, Japan (M.K.).

## References

  * [1] D. H. Weinberg, J. S. Bullock, F. Governato, R. Kuzio de Naray and A. H. G. Peter, _Cold dark matter: controversies on small scales_ , _Proc. Nat. Acad. Sci._ 112 (2015) 12249–12255, [1306.0913].
  * [2] B. Moore, T. R. Quinn, F. Governato, J. Stadel and G. Lake, _Cold collapse and the core catastrophe_ , _Mon. Not. Roy. Astron. Soc._ 310 (1999) 1147–1152, [astro-ph/9903164].
  * [3] W. de Blok, S. S. McGaugh, A. Bosma and V. C. Rubin, _Mass density profiles of LSB galaxies_ , _Astrophys. J. Lett._ 552 (2001) L23–L26, [astro-ph/0103102].
  * [4] S.-H. Oh, W. de Blok, F. Walter, E. Brinks and J. Kennicutt, Robert C., _High-resolution dark matter density profiles of THINGS dwarf galaxies: Correcting for non-circular motions_ , _Astron. J._ 136 (2008) 2761, [0810.2119].
  * [5] N. Amorisco and N. Evans, _Dark matter cores and cusps: the case of multiple stellar populations in dwarf spheroidals_ , _Monthly Notices of the Royal Astronomical Society_ 419 (2012) 184–196.
  * [6] X. Dong, B. Horn, E. Silverstein and A. Westphal, _Simple exercises to flatten your potential_ , _Phys. Rev. D_ 84 (2011) 026011, [1011.4521].
  * [7] R. Kallosh and A. Linde, _Universality Class in Conformal Inflation_ , _JCAP_ 07 (2013) 002, [1306.5220].
  * [8] R. Kallosh, A. Linde and D. Roest, _Universal Attractor for Inflation at Strong Coupling_ , _Phys. Rev. Lett._ 112 (2014) 011303, [1310.3950].
  * [9] R. Kallosh, A. Linde and D. Roest, _Superconformal Inflationary $\alpha$-Attractors_, _JHEP_ 11 (2013) 198, [1311.0472].
  * [10] E. Silverstein and A. Westphal, _Monodromy in the CMB: Gravity Waves and String Inflation_ , _Phys. Rev. D_ 78 (2008) 106003, [0803.3085].
  * [11] L. McAllister, E. Silverstein and A. Westphal, _Gravity Waves and Linear Inflation from Axion Monodromy_ , _Phys. Rev. D_ 82 (2010) 046003, [0808.0706].
  * [12] Y. Nomura, T. Watari and M. Yamazaki, _Pure Natural Inflation_ , _Phys. Lett. B_ 776 (2018) 227–230, [1706.08522].
  * [13] S. Dubovsky, A. Lawrence and M. M. Roberts, _Axion monodromy in a model of holographic gluodynamics_ , _JHEP_ 02 (2012) 053, [1105.3740].
  * [14] J. Soda and Y. Urakawa, _Cosmological imprints of string axions in plateau_ , _Eur. Phys. J. C_ 78 (2018) 779, [1710.00305].
  * [15] N. Kitajima, J. Soda and Y. Urakawa, _Gravitational wave forest from string axiverse_ , _JCAP_ 10 (2018) 008, [1807.07037].
  * [16] I. Bogolyubsky and V. Makhankov, _Lifetime of Pulsating Solitons in Some Classical Models_ , _Pisma Zh. Eksp. Teor. Fiz._ 24 (1976) 15–18.
  * [17] M. Gleiser, _Pseudostable bubbles_ , _Phys. Rev. D_ 49 (1994) 2978–2981, [hep-ph/9308279].
  * [18] E. J. Copeland, M. Gleiser and H.-R. Muller, _Oscillons: Resonant configurations during bubble collapse_ , _Phys. Rev. D_ 52 (1995) 1920–1933, [hep-ph/9503217].
  * [19] A. Vaquero, J. Redondo and J. Stadler, _Early seeds of axion miniclusters_ , _JCAP_ 04 (2019) 012, [1809.09241].
  * [20] E. W. Kolb and I. I. Tkachev, _Nonlinear axion dynamics and formation of cosmological pseudosolitons_ , _Phys. Rev. D_ 49 (1994) 5040–5051, [astro-ph/9311037].
  * [21] E. W. Kolb and I. I. Tkachev, _Axion miniclusters and Bose stars_ , _Phys. Rev. Lett._ 71 (1993) 3051–3054, [hep-ph/9303313].
  * [22] M. A. Amin, R. Easther, H. Finkel, R. Flauger and M. P. Hertzberg, _Oscillons After Inflation_ , _Phys. Rev. Lett._ 108 (2012) 241302, [1106.3335].
  * [23] M. Kawasaki, W. Nakano and E. Sonomoto, _Oscillon of Ultra-Light Axion-like Particle_ , _JCAP_ 2001 (2020) 047, [1909.10805].
  * [24] J.-P. Hong, M. Kawasaki and M. Yamazaki, _Oscillons from Pure Natural Inflation_ , _Phys. Rev._ D98 (2018) 043531, [1711.10496].
  * [25] B. Li, T. Rindler-Daller and P. R. Shapiro, _Cosmological Constraints on Bose-Einstein-Condensed Scalar Field Dark Matter_ , _Phys. Rev. D_ 89 (2014) 083536, [1310.6061].
  * [26] N. K. Porayko et al., _Parkes Pulsar Timing Array constraints on ultralight scalar-field dark matter_ , _Phys. Rev. D_ 98 (2018) 102002, [1810.03227].
  * [27] M. Nori, R. Murgia, V. Iršič, M. Baldi and M. Viel, _Lyman $\alpha$ forest and non-linear structure characterization in Fuzzy Dark Matter cosmologies_, _Mon. Not. Roy. Astron. Soc._ 482 (2019) 3227–3243, [1809.09619].
  * [28] J. Ollé, O. Pujolàs and F. Rompineve, _Oscillons and Dark Matter_ , _JCAP_ 02 (2020) 006, [1906.06352].
  * [29] M. Ibe, M. Kawasaki, W. Nakano and E. Sonomoto, _Decay of I-ball/Oscillon in Classical Field Theory_ , _JHEP_ 04 (2019) 030, [1901.06130].
  * [30] H.-Y. Zhang, M. A. Amin, E. J. Copeland, P. M. Saffin and K. D. Lozanov, _Classical Decay Rates of Oscillons_ , _JCAP_ 07 (2020) 055, [2004.01202].
  * [31] M. A. Amin and Z.-G. Mou, _Electromagnetic Bursts from Mergers of Oscillons in Axion-like Fields_ , 2009.11337.
  * [32] J. H. Buckley, P. B. Dev, F. Ferrer and F. P. Huang, _Fast radio bursts from axion stars moving through pulsar magnetospheres_ , 2004.06486.
  * [33] A. Prabhu, _Optical Lensing by Axion Stars: Observational Prospects with Radio Astrometry_ , 2006.10231.
  * [34] M. Kawasaki, T. Kobayashi and F. Takahashi, _Non-Gaussianity from Curvatons Revisited_ , _Phys. Rev. D_ 84 (2011) 123506, [1107.6011].
  * [35] S. Kasuya, M. Kawasaki and F. Takahashi, _I-balls_ , _Phys. Lett. B_ 559 (2003) 99–106, [hep-ph/0209358].
  * [36] M. Kawasaki, F. Takahashi and N. Takeda, _Adiabatic Invariance of Oscillons/I-balls_ , _Phys. Rev. D_ 92 (2015) 105024, [1508.01028].
  * [37] N. G. Vakhitov and A. A. Kolokolov, _Stationary solutions of the wave equation in a medium with nonlinearity saturation_ , _Radiophysics and Quantum Electronics_ 16 (1973) 783–789.
  * [38] T. D. Lee and Y. Pang, _Nontopological solitons_ , _Phys. Rept._ 221 (1992) 251–350.
  * [39] K. Mukaida, M. Takimoto and M. Yamada, _On Longevity of I-ball/Oscillon_ , _JHEP_ 03 (2017) 122, [1612.07750].
  * [40] K. Mukaida and M. Takimoto, _Correspondence of I- and Q-balls as Non-relativistic Condensates_ , _JCAP_ 1408 (2014) 051, [1405.3233].
  * [41] M. Ibe, M. Kawasaki, W. Nakano and E. Sonomoto, _Fragileness of Exact I-ball/Oscillon_ , _Phys. Rev. D_ 100 (2020) 125021, [1908.11103].
  * [42] T. Hiramatsu, M. Kawasaki and F. Takahashi, _Numerical study of Q-ball formation in gravity mediation_ , _JCAP_ 1006 (2010) 008, [1003.1779].
  * [43] M. A. Amin, R. Easther and H. Finkel, _Inflaton Fragmentation and Oscillon Formation in Three Dimensions_ , _JCAP_ 12 (2010) 001, [1009.2505].
  * [44] M. A. Amin and P. Mocz, _Formation, gravitational clustering, and interactions of nonrelativistic solitons in an expanding universe_ , _Phys. Rev. D_ 100 (2019) 063507, [1902.07261].
  * [45] K. D. Lozanov and M. A. Amin, _Gravitational perturbations from oscillons and transients after inflation_ , _Phys. Rev. D_ 99 (2019) 123504, [1902.06736].

##  Appendix A Suppression of Power Spectrum by Energy Conservation

We estimate the effect of energy conservation, which suppresses the power
spectrum on the large scale as Eq. (2.32). Here, we simply assume that all the
energy of ALP is confined inside oscillons and oscillons have the same energy
$M_{o}$. In this case, Eq.(2.27) contains the following sum of phases:

$\displaystyle\psi=\sum_{i}^{N}\exp(-i\bm{q}_{i}\cdot\bm{k}),$ (A.1)

where index $i$ refers to $i$-th oscillon and $N$ is the total oscillon
number.

At first, we revisit the calculation when the volume $V$ is smaller than the
horizon size at oscillon formation, that is, we expect that the oscillons are
randomly distributed. Then, each phase factor $\bm{x}_{i}\cdot\bm{k}$ is a
random variable, and the sum is estimated in the large number limit
($N\to\infty$) as

$\displaystyle\left|\sum_{i}^{N}\exp(-i\bm{x}_{i}\cdot\bm{k})\right|^{2}=N.$
(A.2)

Next, we include the effect of energy conservation. Let us take $R$ as the
scale for energy conservation where the energy transfer occurs only for scales
smaller that $R$ and oscillon distribution is almost similar on the scales
larger than $R$. The oscillons are randomly distributed inside boxes with the
volume $R^{3}$, that is, the position of a oscillon
$\bm{q}=(q_{x},q_{y},q_{z})$ has the uniform distribution
$q_{w}\in[j_{w}R,(j_{w}+1)R]$ with $w=x$, $y$ and $z$. The index
$\bm{j}=(j_{x},j_{y},j_{z})$ refers to the different box in which the energy
of oscillons is conserved independently. The energy conservation requires that
the number of oscillons in each box is universal constant $N_{o}$. We take the
total space with volume $V=L^{3}$ and calculate the power spectrum
$\mathcal{P}_{k}$ with $k=2\pi/L$, where $L$ is larger than the single box
size, $J\equiv L/R\gg 1$. Note that the total number of oscillons is now
$N=N_{o}J^{3}$.

The average over the position of oscillons is given by

$\displaystyle\braket{\cdots}_{o}=\prod_{\bm{j},n_{j}}\left(\int_{j_{x}R}^{(j_{x}+1)R}\frac{\text{d}q_{x,n_{j}}}{R}\int_{j_{y}R}^{(j_{y}+1)R}\frac{\text{d}q_{y,n_{j}}}{R}\int_{j_{z}R}^{(j_{z}+1)R}\frac{\text{d}q_{z,n_{j}}}{R}\right)(\cdots),$
(A.3)

where index $n_{j}$ refers to oscillons inside the $\bm{j}$-th box and $n_{j}$
runs from $1$ to $N_{o}$. We take $\bm{k}=(2\pi/L,0,0)$ without loss of
generality assuming the isotropic oscillon distribution. The average of $\psi$
is calculated as

$\displaystyle\braket{\psi}_{o}$
$\displaystyle=\prod_{j_{x}^{\prime},n_{j}^{\prime}}\left(\int_{j_{x}^{\prime}R}^{(j_{x}^{\prime}+1)R}\frac{\text{d}q_{n_{j}^{\prime}}}{R}\right)\sum_{j_{x},n_{j}}^{J,N_{o}}\exp(-2\pi
iq_{n_{j}}/L),$ (A.4)
$\displaystyle=\sum_{j_{x},n_{j}}^{J,N_{o}}\frac{1}{-ikR}\left[\exp(-2\pi
i(j_{x}+1)/J)-\exp(-2\pi ij_{x}/J)\right],$ (A.5)
$\displaystyle=N_{o}\frac{\exp(-2\pi i/J)-1}{-ikR}\sum_{j_{x}}^{J}\exp(-2\pi
ij_{x}/J)=0,$ (A.6)

where we use the orthogonal condition $\sum_{j=1}^{J}\exp(-i2\pi j/J)=0$. The
variance is calculated as

$\displaystyle\braket{}{\psi}{{}^{2}}_{o}$
$\displaystyle=\prod_{j_{x}^{\prime\prime},n_{j}^{\prime\prime}}\left(\int_{j_{x}^{\prime\prime}R}^{(j_{x}^{\prime\prime}+1)R}\frac{\text{d}q_{n_{j}^{\prime\prime}}}{R}\right)\left|\sum_{j_{x},n_{j}}^{J,N_{o}}\exp(-2\pi
iq_{n_{j}}/L)\right|^{2},$ (A.7)
$\displaystyle=\prod_{j^{\prime\prime},n_{j}^{\prime\prime}}\left(\int_{j^{\prime\prime}R}^{(j^{\prime\prime}+1)R}\frac{\text{d}q_{n_{j}^{\prime\prime}}}{R}\right)\left[\sum_{j,j^{\prime}=1}^{J}\sum_{n_{j},n_{j^{\prime}}^{\prime}}^{N_{o}}e^{-2\pi
i(q_{n_{j}}-q_{n_{j^{\prime}}^{\prime}})/L}\right],$ (A.8)

where we omit the subscript ”$x$” of $j$. The average over $\text{d}q_{n_{j}}$
is 1 for the terms with $j=j^{\prime}$ and $n_{j}=n_{j^{\prime}}^{\prime}$
while the other combination can be integrated separately for $q_{n_{j}}$ and
$q_{n^{\prime}_{j^{\prime}}}$ as

$\displaystyle\int_{j_{x}R}^{(j_{x}+1)R}\frac{\text{d}q_{n_{j}}}{R}e^{-2\pi
iq_{n_{j}}/L}=\frac{e^{-2\pi i/J}-1}{-ikR}e^{-2\pi ij_{x}/J}=\epsilon e^{-2\pi
ij_{x}/J},$ (A.9)

where we define the suppression factor:

$\displaystyle\epsilon=\frac{J}{-i2\pi}(e^{-2\pi i/J}-1).$ (A.10)

The integration is performed as

$\displaystyle\braket{}{\psi}{{}^{2}}_{o}$
$\displaystyle=N_{o}J+N_{o}(N_{o}-1)|\epsilon|^{2}J+N_{o}^{2}|\epsilon|^{2}\sum_{j,j^{\prime}(\neq
j)}^{J}e^{-2\pi i(j-j^{\prime})/J},$ (A.11)
$\displaystyle=N_{o}J+N_{o}(N_{o}-1)|\epsilon|^{2}J+N_{o}^{2}|\epsilon|^{2}\left(\left|\sum_{j=1}^{J}e^{-2\pi
ij/J}\right|^{2}-J\right),$ (A.12)
$\displaystyle=N_{o}J(1-|\epsilon|^{2})=N\left[1-\left(\frac{2}{Rk}\right)^{2}\sin^{2}\left(\frac{Rk}{2}\right)\right].$
(A.13)

For the modes smaller than the box size $Rk\gg 1$, the variance coincides with
the Poisson distribution, $\braket{}{\psi}{{}^{2}}_{o}\to N$. On the other
hand, for the modes larger than the box size $Rk\ll 1$, the perturbation is
suppressed by energy conservation.