eq

$\displaystyle\begin{split}\BODY\end{split}$ (1)

# Nonperturbative structure in coupled axion sectors
and implications for direct detection

David Cyncynates<EMAIL_ADDRESS>Department of Physics, University of
Washington, Seattle, Washington 98195, USA  Stanford Institute for Theoretical
Physics, Stanford University, Stanford, California 94305, USA  Olivier Simon
<EMAIL_ADDRESS>Stanford Institute for Theoretical Physics, Stanford
University, Stanford, California 94305, USA  Jedidiah O. Thompson
<EMAIL_ADDRESS>Stanford Institute for Theoretical Physics, Stanford
University, Stanford, California 94305, USA  Zachary J. Weiner<EMAIL_ADDRESS>Department of Physics, University of Washington, Seattle, Washington 98195,
USA

###### Abstract

Pairs of misalignment-produced axions with nearby masses can experience a
nonlinear resonance that leads to enhanced direct and astrophysical signatures
of axion dark matter. In much of the relevant parameter space, self-
interactions cause axion fluctuations to become nonperturbative and to
collapse in the early Universe. We investigate the observational consequences
of such nonperturbative structure in this “friendly axion” scenario with $3+1$
dimensional simulations. Critically, in a substantial fraction of parameter
space we find that nonlinear dynamics work to equilibrate the abundance of the
two axions, making it easier than previously expected to experimentally
confirm the existence of a resonant pair. We also compute the gravitational
wave emission from friendly axion dark matter; while the resulting stochastic
background is likely undetectable for axion masses above
$10^{-22}\,\text{eV}$, the polarization of the cosmic microwave background
does constrain possible hyperlight, friendly subcomponents. Finally, we
demonstrate that dense, self-interaction–bound oscillons formed during the
period of strong nonlinearity are driven by the homogeneous axion background,
enhancing their lifetime beyond the in-vacuum expectation.

##  I Introduction

Axions are some of the best-motivated extensions to the Standard Model (SM).
The simplest such extension, the QCD axion, was originally proposed to solve
the strong $CP$ problem [1, 2, 3, 4], but it has since been realized that
axions are common in many theories beyond the SM (BSM) [5, 6, 7]. One
particularly important example is string theory, which generically predicts a
large number of light axions coupled weakly to the SM [8]. The possibility of
such a “string axiverse” is of particular interest because it offers a
potential low-energy window into extremely high-energy physics.

The simplest nonthermal production mechanism for a cosmological abundance of
axions is the misalignment mechanism [9, 10, 11, 12, 13]. Any axion with a
mass lighter than the Hubble scale during inflation would be seeded in an
approximately homogeneous state displaced from the vacuum. It would then
remain frozen at this “misaligned” field value until the expansion rate drops
below its mass, at which point it begins to coherently oscillate about the
minimum of its potential. Barring substantial sources of isocurvature, axions
have large-scale density perturbations that track the adiabatic fluctuations
also seeded during inflation and are thus a viable candidate for the observed
dark matter (DM) or a subcomponent thereof.

An axion’s potential is generically nonlinear, but at late times all axions
with a mass $m$ much larger than the present-day Hubble rate ($m\gg H_{0}$)
oscillate near the bottom of their potential and may be treated as free,
massive fields. This is not, however, a valid assumption at early times, and
it has become increasingly apparent that nonlinearities in an axion’s
potential can have an outsized impact on many late-time observables (see,
e.g., Refs. [14, 15, 16, 17, 18, 19, 20, 21]). If the dark matter comprises a
single axion, these early-time dynamics can strongly enhance structure on
scales that enter the horizon when the Hubble rate $H$ is approximately the
axion mass $m$ [16].

More generally, a string axiverse may consist of many axions interacting with
each other through a joint potential, and recent work has shown that when any
two of these have similar masses (within a factor of roughly 2) a new type of
efficient, resonant energy transfer is possible [18]. This mechanism, dubbed
“friendship” due to the necessary mild coincidence of masses, transfers energy
from an axion with a high decay constant to one with a lower decay constant.
Since an axion’s couplings to the SM are generically inversely proportional to
its decay constant, the mechanism boosts the abundance of the more strongly
coupled axion. In other words, friendly axion dark matter can be significantly
more visible to direct detection experiments than would be expected for either
axion individually.

In this paper, we follow up on the work of Ref. [18] with a suite of $3+1$
dimensional numerical simulations, corroborating its findings and extending
the results to the strongly nonlinear regime. As anticipated in that work,
large spatial inhomogeneities significantly modify the results of a
homogeneous analysis. Nonperturbative fluctuations collapse into dense
oscillons, nontopological field configurations bound by self-interactions [22,
23, 24, 25, 26, 27, 28, 29, 30, 31]. The oscillons quench the resonant
amplification and mediate energy transfer between the friendly pair, leading
to approximate energy density equipartition over a broad range of parameters.
In contrast to expectations from a homogeneous analysis, the enhanced
visibility of one axion therefore does not come at the expense of the other’s
detectability. In sum, nonlinear dynamics make the friendly axion model both
more predictive (by being less parameter-dependent) and more identifiable
(because both axions would be detectable).

This paper is divided as follows. In Sec. II we review the friendly axion
model and results within the spatially homogeneous approximation. Section III
presents the extension of these results into the nonlinear regime using
numerical simulations, with a primary focus on the late-time abundance as
relevant to direct detection experiments. We also investigate gravitational
wave signatures in these scenarios, which, while not promising if the friendly
axions make up all of the dark matter, are relevant for hyperlight
subcomponents. Finally, we study a novel driving effect in which oscillons
resonantly siphon energy from the axion background, parametrically enhancing
their lifetime. We conclude in Sec. IV, putting this work into the broader
context of the landscape of nonlinear axion models. For completeness and ease
of readability, we relegate an extended discussion of methodology and
additional results to the appendices. Appendix A enumerates the system of
evolution equations and the details of our numerical implementation, and App.
B expands upon our discussion of bound axion states.

##  II Review of friendly axions

As a concrete and illustrative model, Ref. [18] focuses on a simple two-axion
potential with two instanton contributions:111 Throughout, we work in units
where $\hbar=c=1$. We also define the reduced Planck mass
$M_{\mathrm{pl}}=1/\sqrt{8\pi G}\approx 2.44\times 10^{18}\,\mathrm{GeV}$.

$\displaystyle\begin{split}V(\phi_{S},\phi_{L})&=\Lambda_{1}^{4}\left[1-\cos\left(\frac{\phi_{S}}{f_{S}}+\frac{\phi_{L}}{f_{L}}\right)\right]\\\
&\hphantom{{}={}}+\Lambda_{2}^{4}\left(1-\cos\frac{\phi_{L}}{f_{L}}\right).\end{split}$
(2)

The canonically normalized axion field variables $\phi_{S}$ and $\phi_{L}$ are
naturally recast as angular variables via the definition
$\theta_{S}\equiv\phi_{S}/f_{S}$ and $\theta_{L}\equiv\phi_{L}/f_{L}$.
Redefining $\Lambda_{1}^{4}\equiv m^{2}f^{2}$ and
$\Lambda_{2}^{4}\equiv\mu^{2}m^{2}\mathcal{F}^{2}f^{2}$, the axion masses
are222 The interaction-basis axions $\phi_{S}$ and $\phi_{L}$ are not exact
mass eigenstates, making this definition ambiguous. For $\mathcal{F}\gg 1$,
the distinction between the two bases is small, and so we often neglect the
distinction in our heuristic discussions. The effect is not, however,
quantitatively negligible for all parts of the parameter space we consider,
and it is always included in our results.  $m_{S}=m$ and $m_{L}=\mu m_{S}$ and
their decay constants are $f_{S}=f$ and $f_{L}=\mathcal{F}f$, respectively. In
terms of these variables, Eq. 2 takes the form

$\displaystyle\begin{split}V(\theta_{L},\theta_{S})&=m^{2}f^{2}\Big{[}\left(1-\cos\left(\theta_{S}+\theta_{L}\right)\right)\\\
&\hphantom{{}={}m^{2}f^{2}\Big{[}}+\mu^{2}\mathcal{F}^{2}\left(1-\cos\theta_{L}\right)\Big{]}.\end{split}$
(3)

We focus on the range $\mathcal{F}>1$ where $f_{S}<f_{L}$, and we refer to
$\phi_{S}$ and $\phi_{L}$ as the “short” and “long” axion respectively in
reference to the size of their decay constants. (The regime with $f_{S}>f_{L}$
does not exhibit nonlinear resonances.) The short and long axions then form a
“friendly pair” when $0.7\lesssim\mu<1$, corresponding to an $\mathcal{O}(1)$
coincidence in their masses. While Eq. 2 might represent a subsector of a much
larger axiverse, the dynamics of the friendly pair of interest are insensitive
to possible couplings to other axions barring additional coincidences in mass.
Namely, only the relative frequency of coupled oscillators determines the
efficiency of energy transfer between them, so the actual instanton scales
$\Lambda_{i}$ and decay constants $f_{i}$ matter only insofar as they
(together) determine the axion masses.

In the early Universe, the misalignment mechanism initializes each axion at an
approximately spatially homogeneous value away from the late-time minimum; a
natural assumption, barring anthropic and other considerations, is that
$\theta_{I}(t_{\text{initial}})=\mathcal{O}(1)$, where the capital index $I$
runs over axion flavors. The axions remain frozen at their misaligned values
until the Hubble rate $H$ drops below their masses. Since the two axion masses
are comparable, the long axion initially has $\mathcal{O}(\mathcal{F}^{2})$
times more energy than the short axion. In the absence of couplings between
the axions, this imbalance would persist to their present-day abundance.

The same conclusion holds for coupled axions as well, so long as the masses of
the axions are well separated. At large field values, however, interactions
can substantially shift the axion oscillation frequency from its ground state
value. Reference [18] showed that coupled axions with a decay constant
hierarchy $\mathcal{F}\gtrsim 3$ and sufficiently close masses
$0.75\lesssim\mu<1$ tend to align their frequencies in a process called
autoresonance, illustrated in Fig. 1.

Figure 1:  Homogeneous dynamics (i.e., ignoring the effect of spatial
fluctuations) of friendly axions with mass ratio $\mu=0.75$ and decay constant
ratio $\mathcal{F}=20$. Here $t$ and $\tau$ are cosmic and conformal time
coordinates, respectively. The top panel depicts the evolution of the energy
density in the short (pink) and long (blue) axions, while the bottom panel
displays the field values $\phi_{I}/f_{I}$. Transparent curves of the same
color denote the corresponding results for uncoupled axions in the top panel.

Specifically, interactions drive the short axion (with the smaller decay
constant) to dynamically adjust its oscillation amplitude to a fixed value in
order to match its frequency to the long axion’s, as evident in the lower
panel of Fig. 1. Consequently, the short axion energy density does not dilute
like cold matter but instead remains fixed (as in the top panel of Fig. 1) by
siphoning energy from the long axion. If the fields remain spatially
homogeneous, this energy transfer runs until backreaction disrupts the precise
phase locking of the two fields. Autoresonance then ends when
$\bar{\rho}_{S}/\bar{\rho}_{L}\simeq 2\mathcal{F}^{2}(1-\mu)^{2}$ for
$\mathcal{F}^{2}\gg(1-\mu^{2})^{-1}$, representing a near-complete transfer of
the available energy density to the short axion. In other words, when
autoresonance runs to completion, the energy density at late times in the dark
sector is virtually entirely in the short axion—an outcome opposite to what
one would expect from free evolution.

The boost to the late-time energy density of the short axion relative to the
scenario of independent axions is of great importance for direct detection
experiments. Laboratory haloscopes probe the couplings of axion DM to SM
states, which are typically higher-dimension operators suppressed by the axion
decay constant $f_{a}$. For example, axions are expected to couple to SM
photons via an interaction of the form:

$\mathcal{L}\supset-\frac{g_{a\gamma\gamma}}{4}\phi
F_{\mu\nu}\tilde{F}^{\mu\nu}$ (4)

where $g_{a\gamma\gamma}=C_{a\gamma\gamma}\alpha_{\text{QED}}/2\pi f_{a}$ is
the axion-photon coupling [13, 32] and $C_{a\gamma\gamma}$ an order-unity
constant. As discussed in Ref. [18], when all axions evolve independently from
$\mathcal{O}(1)$ initial misalignment angles, the final energy density
$\rho_{a,0}$ of each axion is proportional to $f_{a}^{2}$. In this case, the
signal strength $\rho_{a,0}g_{a\gamma\gamma}^{2}$ is roughly independent of
$f_{a}$; as such, at a given mass any axion produced by the standard
misalignment mechanism would be similarly hard to see. In a scenario with
friendship however, the boosted late-time energy density of the short axion is
$\bar{\rho}_{S,0}\propto f_{L}^{2}$ when autoresonance completes, but the
coupling to SM photon is $g_{S\gamma\gamma}\simeq\alpha/2\pi f_{S}$. Thus the
signal strength $\bar{\rho}_{S,0}g_{S\gamma\gamma}^{2}$ of the short axion is
enhanced by $\mathcal{F}^{2}$, making it much more accessible to axion
haloscopes. The effect, however, would be reversed for the long axion: its
energy density is suppressed by $\sim\mathcal{F}^{2}$ compared to standard
misalignment scenarios. In this picture, seeing both friendly axions would
therefore be difficult.

The description of autoresonance given so far assumes the fields remain
approximately spatially homogeneous, but large spatial fluctuations in the
axions can prevent the completion of the energy transfer. The coherent
oscillations of the short axion induce a time-dependent effective mass that
resonantly amplifies fluctuations of the short axion, much like that which
characterizes preheating after inflation [33, 34, 35] (see Refs. [36, 37, 38,
39] for reviews) and large misalignment [16]. Large-amplitude fluctuations of
the short axion can collapse under attractive self-interactions into
oscillons—finite-lifetime, nontopological bound structures with densities of
$\mathcal{O}(m^{2}f^{2})$ and radii of $\mathcal{O}(1/m)$. Such oscillons
explore large field values for the short axion and thus continue to experience
large interactions, but, being nonperturbative objects, are difficult to treat
analytically. Reference [18] presented preliminary evidence that oscillon
nucleation occurs for $\mathcal{F}\gtrsim 6$ and that oscillons quench
autoresonance if they form early enough, setting a limit on the energy density
transfer for $\mathcal{F}\gtrsim 20$. The remainder of this paper investigates
the impact of the nonlinear dynamics of autoresonance and oscillon formation
on the predictions of friendly axion scenarios through the use of $3+1$
dimensional numerical simulations.

##  III Results

We now present numerical solutions for the fully nonlinear, friendly axion
system. We implement numerical simulations of the axions’ classical equations
of motions with pystella [40, 41, 42], discretizing these equations onto a
$3D$, periodic, regularly spaced grid, computing spatial derivatives via
fourth-order centered differencing, and utilizing a fourth-order Runge-Kutta
method for time integration. Further details are provided in App. A. Except
where otherwise stated, all results use grids with $N^{3}=1024^{3}$ points, a
comoving side length $L=1.5\,\pi/m$ and conformal timestep $\Delta\tau=\Delta
x/10=L/10N$. The simulations begin with a numerical solution to the linearized
system of equations starting at a time when the Hubble rate $H\ll m$ (see App.
A for further details). The $3D$ evolution begins when $H=m$, corresponding to
a conformal time $m\tau_{m}=1$ and cosmic time $mt_{m}=1/2$. The scale factor
is normalized relative to $a_{m}\equiv a(t_{m})$.

Of the free parameters in the model, the decay constant ratio $\mathcal{F}$
has the strongest effect on the dynamics. The mass ratio and initial
misalignments mainly determine whether or not autoresonance occurs at all,
whereas the decay constant ratio determines the size of nonlinear backreaction
and even whether fluctuations are sizeably enhanced at all. Therefore, for
most simulations we pick fiducial values $\mu=0.75$,
$\theta_{L}(0,\mathbf{x})=0.8\,\pi$, and $\theta_{S}(0,\mathbf{x})=0$, and run
simulations for varying values of $\mathcal{F}$.333 So long as we choose
$\theta_{L}(0,\mathbf{x})$ large enough that the axions experience
autoresonance, the initial misalignment angles are essentially inconsequential
[18]. On the other hand, the choice of the relatively detuned mass ratio
$\mu=0.75$ is made to reduce the runtime of the simulations, as smaller $\mu$
cause perturbations to grow faster (see Appendix C of Ref. [18]) and shortens
the oscillon lifetime (explained in Sec. III.3 below).

Figure 2:  Projected density contrast for the short axion, Eq. 5, at eight
snapshots illustrating key moments during the course of a simulation with mass
ratio $\mu=0.8$ and decay constant ratio $\mathcal{F}=50$. The onset of strong
nonlinearity is observed between the third and fourth panel, followed by the
collapse of large overdensities into spherically symmetric oscillon
configurations. At this point, autoresonance ceases at the homogeneous level,
and energy stops flowing from the long axion to the short axion. The oscillons
persist, driven by the long axion in a form of localized autoresonance (see
Sec. III.3). Between the sixth, seventh, and eighth panel, the oscillons
dissipate, leaving remnant overdensities that would eventually seed small-
scale galactic substructure.

For $\mathcal{F}\lesssim 6$, spatial perturbations do not grow large enough to
form oscillons and the results of the simulations are described completely by
Ref. [18]. For larger $\mathcal{F}$, fluctuations of $\phi_{S}$ indeed
collapse into oscillons as anticipated by Ref. [18]. We present a broad
overview of the dynamics of oscillon formation in Fig. 2, plotting two
dimensional projections of the energy density in the short axion at various
times over the course of a simulation.444 To be explicit, we display the
energy density projected (averaged) along one axis of the simulation volume,
e.g.,  $\displaystyle\rho_{S}(x,y)$
$\displaystyle=\frac{1}{\left\langle\rho_{S}(x,y,z)\right\rangle}\frac{1}{L}\int_{0}^{L}\mathrm{d}z\,\rho_{S}(x,y,z).$
(5) Such a projected quantity presents more information about the full volume
than a single two dimensional slice but also underestimates the magnitude of
overdensities (since, e.g., any given oscillon occupies only a small fraction
of space along the $z$ axis).  The field begins in a nearly homogeneous state
in the first panel, in which the initial adiabatic fluctuations are too small
to be seen. The second and third panels depict the linear enhancement of
fluctuations by parametric resonance as the amplification of local
overdensities. Fluctuations become nonlinear at a time $mt_{\mathrm{nl}}\sim
100$, resulting in large overdensities that quickly collapse under attractive
self-interactions into the oscillons apparent in the fifth panel. These
oscillons radiate energy and begin to dissipate one by one around $mt\gtrsim
2000$. Eventually, no bound objects remain and nonlinear interactions cease to
dominate the dynamics, although significant density fluctuations remain.

The interplay between the persistence of homogeneous autoresonance and the
onset of nonlinearity has important consequences for the final distribution of
energy between the two axions. We discuss these dynamics in Sec. III.1,
comparing to the results of Ref. [18]. In Sec. III.2 we compute the
gravitational wave production from friendly axions, finding possible
signatures for hyperlight subcomponents in the CMB $B$-mode polarization.
Finally, in Sec. III.3 we demonstrate that the oscillons that form continue to
experience autoresonance long after the spatially averaged fields cease to
resonate, and we discuss the implications for oscillon lifetimes.

###  III.1 Evolution of energy densities

Having established the importance of nonlinear dynamics for a large portion of
parameter space, we now investigate how nonlinear density fluctuations impact
the final distribution of energy between the two axions (and, as a
consequence, their relic abundances today). We first study the evolution of
each axion’s energy density in Fig. 3 for three representative values of
$\mathcal{F}$, comparing the result of simulations to that of a homogeneous
analysis. To avoid ambiguities in the final partition of energy densities we
work in the mass basis

$\displaystyle\nu_{h}$ $\displaystyle\equiv\phi_{S}\cos\eta+\phi_{L}\sin\eta,$
(6a) $\displaystyle\nu_{l}$
$\displaystyle\equiv-\phi_{S}\sin\eta+\phi_{L}\cos\eta,$ (6b)
$\displaystyle\cos 2\eta$
$\displaystyle\equiv\frac{1-\mu^{2}-\mathcal{F}^{-2}}{\sqrt{4\mathcal{F}^{-2}+(1-\mu^{2}-\mathcal{F}^{-2})^{2}}},$
(6c)

where the heavy state $\nu_{h}$ is composed mostly of the short axion, and the
light state $\nu_{l}$ is composed mostly of the long axion in the limit
$\mathcal{F}\gg 1$ (see Appendix A of Ref. [18] for a complete discussion).

Figure 3:  Evolution of energy density in the mass eigenstates (Eq. 6) for
three simulations with mass ratio $\mu=0.75$ and decay constant ratios
$\mathcal{F}=10$, $20$, and $40$ by row. Each panel depicts the simulation
result for the heavy and light states ($\nu_{h}$ and $\nu_{l}$) in pink and
blue, respectively, as well as the corresponding results for a purely
homogeneous calculation in thin black and gray. The shaded region denotes the
time when order-one density fluctuations are present, which we define as the
time when more than 5% of the energy in the short axion resides in
overdensities $\rho_{S}(\mathrm{x})/\bar{\rho}_{S}>10$. Shortly after these
large overdensities form, they either dissipate or coalesce into oscillons, so
the gray regions are decent proxies for the presence of oscillons. Note that
in contrast to Fig. 1 we here plot the energy density of the mass-basis fields
rather than the interaction-basis ones.

Each panel exhibits an initial phase of homogeneous, autoresonant energy
transfer and the onset of nonlinearity that quenches autoresonance, at which
point the energy density departs from the trend of the homogeneous result.
From analytic estimates of growth rate for the fastest growing mode [18],
nonlinearity occurs at approximately555This result accounts for both Hubble
friction and the slight decay of the initial metric perturbation before the
fastest-growing mode starts growing (see Eqs. C17 and C18 in Ref. [18] and the
surrounding discussion, fixing $\delta\omega=\mu-1$).

$\displaystyle mt_{\mathrm{nl}}$ $\displaystyle\approx
17.6\frac{1-0.1\log(1-\mu)}{1-\mu},$ (7)

in good agreement with $mt_{\mathrm{nl}}\approx 80$ observed in Fig. 3. The
ultimate partitioning of energy depends primarily on the precise timing of
nonlinearity and oscillon formation relative to the (would-be) completion of
autoresonance, a point which we detail below. We now describe these two
regimes of $\mathcal{F}$ in detail.

For $6\lesssim\mathcal{F}\lesssim 20$, oscillons nucleate after the short
axion’s energy density first exceeds the long axion’s. At roughly the same
time, autoresonance ends and $\bar{\rho}_{h}$ ceases to be roughly constant,
instead decaying approximately as $a^{-3}$ like nonrelativistic matter.
Contrary to the homogeneous analysis of Ref. [18], however, we observe in this
range that nonperturbative dynamics in fact enable energy transfer from the
short axion back to the long axion, resulting in late-time
(near-)equilibration of the two axion energy densities. This phenomenon is
most evident in the top panel of Fig. 3 ($\mathcal{F}=10$), where the heavy
and light axions’ energy densities asymptote toward a common value.

Interactions between the two axions are strongest where the field values are
largest, suggesting that oscillons play a key role in reversing energy
transfer. Inside an oscillon, the field amplitude oscillates with a period
$\omega<m$ due to its binding energy. Since the long axion’s natural frequency
is $\mu m<m$, an oscillon can provide a locus for more efficient energy
transfer from the short axion back to the long axion.666 In fact, during
autoresonance the short axion is driven at exactly the frequency $\mu m$. When
fluctuations grow nonperturbative the oscillon frequencies will thus remain
close to $\mu m$.  Indeed, for most decay constant ratios
$6\lesssim\mathcal{F}\lesssim 20$, the end of autoresonance and formation of
oscillons is associated with a substantial transfer of energy to the light
axion. For $6\lesssim\mathcal{F}\lesssim 10$, the final stage of energy
transfer to the light axion occurs in discrete jumps that appear to coincide
with the death of individual oscillons. In all cases, we observe that most of
the radiation from the heavy axion into the light axion is into
semirelativistic modes, as one would expect if oscillons are responsible for
equilibration. However, at larger $\mathcal{F}$ equipartition is nearly
achieved by the time oscillons form anyway; the subsequent evolution is more
continuous, obfuscating any association between oscillon death and energy
transfer. While nonlinear effects are evidently crucial, identifying the
specific mechanism for energy flow in general is challenging and would require
software infrastructure and computational resources—in particular, substantial
storage of $3D$ snapshots and analysis thereof—well beyond the reach of the
present work.

The middle panel with $\mathcal{F}=20$ represents the marginal case where
oscillons form at nearly the exact time that the heavy axion’s energy density
first reaches that of the light axion. For larger values $\mathcal{F}\gtrsim
20$, $t_{\mathrm{nl}}$ and oscillon formation occur before the heavy axion
dominates the sector’s energy density, terminating autoresonant energy
transfer to the heavy axion. As shown in the bottom panel of Fig. 3
($\mathcal{F}=40$), the energy density in both axions then decays as
approximately $a^{-3}$. In this case the backreaction effects at play for
smaller $\mathcal{F}$ are too suppressed to enable substantial energy transfer
by the oscillons. The trends for yet larger $\mathcal{F}$ are qualitatively
similar: parametric resonance proceeds at the same rate and oscillons form at
a similar time. The final ratio of energy densities
$\bar{\rho}_{h}/\bar{\rho}_{l}$ thus receives a constant boost due to the
period of autoresonance but still decreases as $1/\mathcal{F}^{2}$.

Having discussed the dynamics that control the distribution of energy between
the two axions, we now summarize the full $\mathcal{F}$-dependence of the
late-time energy fractions,

$\displaystyle\Xi_{I}$
$\displaystyle\equiv\left.\frac{\bar{\rho}_{I}}{\bar{\rho}_{h}+\bar{\rho}_{l}}\right|_{\text{late
time}},$ (8)

where $I=h,l$. The final partitioning changes qualitatively at a critical
decay constant ratio $\mathcal{F}_{\star}$ for which nonlinearities become
important (at $t_{\mathrm{nl}}$) just as the heavy axion’s energy density
first matches the light one’s (via autoresonant energy transfer). From our
simulations we find $\mathcal{F}_{\star}\approx 20$ for $\mu=0.75$; this value
depends on the mass ratio in the same manner as $t_{\mathrm{nl}}$ (c.f. Eq.
7). This timing separates two distinct regimes: one of near-equilibration due
to nonlinear effects at $\mathcal{F}<\mathcal{F}_{\star}$ and a
$1/\mathcal{F}^{2}$ suppression of the heavy-axion abundance via the early end
of autoresonance at larger $\mathcal{F}$. Both regimes are well captured by

$\displaystyle\Xi_{h}$
$\displaystyle\sim\begin{dcases}\frac{1}{2}&6\lesssim\mathcal{F}\lesssim\mathcal{F}_{\star}\\\
\frac{1}{1+1.3(\mathcal{F}/\mathcal{F}_{\star})^{2}}\hphantom{1-}&\mathcal{F}\gtrsim\mathcal{F}_{\star}\end{dcases}$
(9a) $\displaystyle\Xi_{l}$
$\displaystyle\sim\begin{dcases}\frac{1}{2}&6\lesssim\mathcal{F}\lesssim\mathcal{F}_{\star}\\\
1-\frac{1}{1+1.3(\mathcal{F}/\mathcal{F}_{\star})^{2}}&\mathcal{F}\gtrsim\mathcal{F}_{\star}\end{dcases},$
(9b)

including an empirical factor of $1.3$ that best fits the results from
simulations. We display the corresponding quantities computed directly from
simulations in Fig. 4.

Figure 4:  Late-time partition of energy between the friendly axions
($\Xi_{I}$, Eq. 8) as a function of the decay constant ratio $\mathcal{F}$,
evaluated at the end of the simulations (long after energy transfer ceases).
The mass ratio is fixed to $\mu=0.75$. Results for the heavy and light axions
are respectively in pink and blue for the simulations and black and gray for
the homogeneous computation. The dashed, dark pink and blue curves depict the
empirical fit of Eq. 9. The shaded region indicates the range of $\mathcal{F}$
for which oscillons form after the heavy axion’s energy overtakes the light
one’s (in which case nonlinear effects return a non-negligible amount of
energy to the light axion).

For $\mathcal{F}\lesssim 6$, the energy densities indeed match those predicted
by the homogeneous theory, which are computed in full in Ref. [18]. For such
low decay constant ratios, autoresonance is too brief to support parametric
resonance long enough for perturbations to reach order unity. For
$6\lesssim\mathcal{F}\lesssim 20$, the energy density of the light axion,
instead of being entirely depleted, remains within a factor of between 1 and 4
of the heavy axion’s energy density with a precise dependence on $\mathcal{F}$
beyond the sophistication of Eq. 9. At larger $\mathcal{F}\gtrsim 20$, the
$1/\mathcal{F}^{2}$ scaling takes over, parametrically suppressing the heavy
axion’s abundance relative to the light axion’s. Nonetheless, in this range
the heavy axion carries approximately $\mathcal{F}_{\star}^{2}/1.3\sim 310$
times more energy density than it would have had in the absence of
autoresonance.

While the above results fix the mass ratio to $\mu=0.75$, we do not expect our
findings to differ substantially in the parameter space where autoresonance
occurs. However, we do expect the precise thresholds in $\mathcal{F}$ of each
regime to vary mildly with $\mu$. In particular, $\mathcal{F}_{\star}$ varies
in accordance with the change in the timescale of nonlinearity
$t_{\mathrm{nl}}$, which itself varies with $\mu$ via Eq. 7. Limited
simulations corroborate these expectations for $0.70\lesssim\mu\lesssim 0.80$,
the range for which simulation runtimes are tractable, but a more quantitative
assessment of the variation of trends with $\mu$ would require multiple times
more computational resources. For increasingly degenerate axion masses, linear
resonance becomes more relevant and could qualitatively alter our findings for
$0.95\lesssim\mu\leq 1$.

The light axion’s enhanced abundance relative to homogeneous results has
important consequences for direct detection experiments. Although the near-
even partitioning of energy for
$6\lesssim\mathcal{F}\lesssim\mathcal{F}_{\star}$ implies that the heavy axion
is slightly harder to detect than predicted by Ref. [18], it also implies that
the light axion requires only twice the experimental sensitivity as it would
for an uncoupled, misaligned axion of the same mass. The homogeneous
expectation, in contrast, was that the light axion would require
$\mathcal{F}^{2}$ times greater experimental sensitivity to detect. The
nonperturbative equalization of energy density in the sector thus serves to
make the sector as a whole more visible to direct detection experiments.
Observing two axions with similar masses—with the heavier axion substantially
more visible than expectations for single-axion misalignment, and the lighter
one comparably so—is a characteristic signature of friendly dynamics of the
form described here. For the mass range $10^{-10}\;\mathrm{eV}\lesssim
m\lesssim 10^{-3}\;\mathrm{eV}$, many near-future experimental efforts will
probe relevant parameter space (see, e.g., Refs. [43, 44, 45, 46, 47, 48, 49,
50, 51, 52, 53, 54]).

To close, we connect the partitioning of Eq. 9 to present-day abundances by
estimating the net present-day energy density in the sector. The energy
density at horizon crossing is dominated by the light axion, i.e.,
$\bar{\rho}_{\mathrm{tot}}(t_{m})\sim\mu^{2}m^{2}\mathcal{F}^{2}f^{2}\Theta_{L,0}^{2}$,
which subsequently redshifts as $a^{-3}$.777 The $a^{-3}$ redshifting assumes
that the axions are always noninteracting and nonrelativistic. In reality, at
early times the axion interactions during autoresonance cause
$\bar{\rho}_{\mathrm{tot}}$ to redshift slightly slower than $a^{-3}$, and at
later times oscillons radiate mildly relativistic axions such that
$\bar{\rho}_{\mathrm{tot}}$ redshifts slightly faster than $a^{-3}$ (until all
axions become nonrelativistic). Together, these effects amount to only an
$\mathcal{O}(1)$ factor which we neglect for simplicity.  Combined with Eq. 9,
the present-day abundance of the each axion is

$\displaystyle\begin{split}\frac{\Omega_{h,0}}{0.13}&\approx\Xi_{I}\left(\frac{m}{10^{-19}\,\mathrm{eV}}\right)^{1/2}\left(\frac{\mathcal{F}f}{10^{16}\,\mathrm{GeV}}\right)^{2}\left(\mu\Theta_{L,0}\right)^{2},\end{split}$
(10)

with $\Xi_{I}$ set by Eq. 9. Factors accounting for the thermal history of the
SM (i.e., the number of effective relativistic degrees of freedom in the SM
entropy and energy density) change the above result by only an
$\mathcal{O}(1)$ factor over the mass range of interest and are omitted for
simplicity.

###  III.2 Gravitational waves

Rapidly growing fluctuations during and after autoresonance and the resulting
oscillons can both source gravitational waves, again much like the parametric
resonance and oscillon formation that can occur during preheating [55, 56, 57,
58, 59, 60, 61, 62, 63, 64] and single-axion misalignment [16]. In this
section we compute the signal strength generated by friendly axions and
discuss the corresponding constraining power of existing and future
observations. Section A.3 briefly reviews stochastic gravitational wave
backgrounds and the transfer functions required to relate their spectra at
emission to that at the present day.

We begin by estimating the scaling of the peak amplitude of the gravitational
wave spectrum with the short-axion mass $m$ and long axion’s decay constant
$\mathcal{F}f$. Gravitational waves are sourced by the anisotropic part of the
stress tensor (via Eq. 29), whose components scale like the energy density of
the source. The time (relative to $t_{m}$) and wave number (relative to $m$)
of peak emission varies weakly with model parameters (and is entirely
independent of $f$). The spectral abundance of gravitational waves (Eq. 39)
therefore scales with two powers of the fractional energy density of the
source at the time of emission,
$\bar{\rho}_{\mathrm{source}}(t)/\bar{\rho}(t)$. The short axion is the
dominant source of anisotropic stress, which we expect to peak near the time
when the system becomes nonlinear, $t_{\mathrm{nl}}$, when axion gradients are
largest and power is scattered to smaller scales.

To proceed, we follow the model-independent heuristics of Refs. [65, 38].
Approximating the gravitational wave source as a Gaussian peaking at momentum
$k_{\star}$ with width $\sigma$, one may estimate the peak amplitude to be
[65]:

$\displaystyle\Omega_{\mathrm{GW}}(k_{\star})$
$\displaystyle=\frac{27\gamma^{2}\nu^{2}}{\sqrt{\pi}}\frac{k_{\star}}{\sigma}\left(\frac{aH_{p}}{k_{\star}}\right)^{2},$
(11)

at the time the source is maximized. Here $\gamma$ is what fraction of the
Universe’s energy is in the source at the time of the process, $\nu$ measures
how anisotropic the source is, and $H_{p}$ the Hubble parameter at the time of
the process.888 Note that our $\nu^{2}$ corresponds to $\beta w^{2}$ in terms
of the parameters of Ref. [65].  The peak wave number $k_{\star}$ and width
$\sigma$ are straightforward to approximate (or read off of simulation
results), but the anisotropy coefficient $\nu$ is harder to estimate; Ref.
[65] motivates $\nu\sim 10^{-2}$ to $10^{-1}$ for typical processes.
Evaluating Eq. 11 at $t_{\mathrm{nl}}$ and plugging in
$aH=1/\sqrt{2t_{\mathrm{nl}}}$,

$\displaystyle\Omega_{\mathrm{GW}}(k_{\star})$
$\displaystyle\approx\frac{27\nu^{2}}{\sqrt{\pi}}\frac{k_{\star}}{\sigma}\left(\frac{m}{k_{\star}}\right)^{2}\frac{\left(\mu\mathcal{F}\Theta_{L,0}f/M_{\mathrm{pl}}\right)^{4}}{\left[1+1.3\left(\mathcal{F}/\mathcal{F}_{\star}\right)^{2}\right]^{2}},$
(12)

where we have used that the energy density in the short axion at
$t_{\mathrm{nl}}$ is approximately
$1/[1+1.3(\mathcal{F}/\mathcal{F}_{\star})^{2}]$ of the total axion energy
density, and the total axion energy density is given by
$\bar{\rho}_{\mathrm{tot}}(t_{m})(t_{m}/t_{\mathrm{nl}})^{3/2}$. Notice that
the suppression from how far inside the horizon the peak is (the factor of
$[aH_{p}/k_{\star}]^{2}$ in Eq. 11) is exactly compensated by the growth in
time of $\gamma$ (since the homogeneous energy available to source the short
axion redshifts with one fewer power of the scale factor than the SM
radiation). We therefore expect gravitational wave signals from friendly
axions to be only weakly sensitive to the time of nonlinearity
$t_{\mathrm{nl}}$.

By comparing Eq. 12 with the relic abundance (Eq. 10) we may estimate the peak
of the gravitational wave signal as a function of the mass $m$ and relic
abundance of the heavy axion $\Omega_{h}(t_{0})$:

$\displaystyle\begin{split}\Omega_{\mathrm{GW},0}h^{2}&=\frac{27\nu^{2}}{\sqrt{\pi}}\frac{k_{\star}}{\sigma}\left(\frac{m}{k_{\star}}\right)^{2}\left[\Omega_{\mathrm{rad}}(t_{0})h^{2}\right]^{-1/2}\\\
&\hphantom{{}={}}\times\left[\Omega_{h}(t_{0})h^{2}\right]^{2}\left(\frac{m}{H_{100}}\right)^{-1}\end{split}$
(13)

where we have dropped factors of the relativistic degrees of freedom, which
reduce the amplitude by at most a factor of two at early enough times $t_{m}$
such that all SM species are in thermal equilibrium. Here,
$H_{100}=100\,\mathrm{km}/\mathrm{s}/\mathrm{Mpc}=2.13\times
10^{-33}\,\mathrm{eV}$ and $h=H_{0}/H_{100}$. From our simulations we observe
that $k_{\star}/m\approx 9$ and $\sigma\approx k_{\star}/3$. Taking $\nu=1/20$
(for which the estimates agree well with the simulation results) and
considering the regime $\mathcal{F}\lesssim 20$ for which the signal is not
suppressed, the peak amplitude is

$\displaystyle\Omega_{\mathrm{GW},0}(k_{\star})h^{2}$ $\displaystyle\approx
10^{-15}\left(\frac{\Omega_{S}(t_{0})h^{2}}{0.06}\right)^{2}\left(\frac{m}{10^{-21}\,\mathrm{eV}}\right)^{-1},$
(14)

at a present-day frequency of the peak of

$\displaystyle\frac{f_{\mathrm{GW},\star}}{\mathrm{Hz}}$ $\displaystyle\approx
2.8\times
10^{-14}\frac{k_{\star}}{m}\left(\frac{m}{10^{-21}\,\mathrm{eV}}\right)^{1/2}.$
(15)

The amplitude estimate Eq. 14 agrees quite well—within a factor of a few—with
the spectra from simulations when evaluated at $t_{\mathrm{nl}}$. However, the
spectra evaluated at the end of the simulation (after all gravitational wave
production has concluded) are about an order of magnitude larger due to
factors not captured by these simplistic estimates such as the time evolution
of the source after $t_{\mathrm{nl}}$.

From Eqs. 14 and 15 we see that adjusting the mass $m$ to change the peak
frequency by some factor modulates the gravitational wave power spectrum by
two powers of that factor. Signals that could be visible at pulsar timing
arrays (with frequencies of order $10^{-9}\,\mathrm{Hz}$) could therefore not
exceed amplitudes of $10^{-23}$ (about eight orders of magnitude below the
projected sensitivity of the Square Kilometer Array [66, 67]) without
requiring an axion abundance that would overclose the Universe. This is an
unfortunate consequence of the source both being short-lived and redshifting
more slowly than the SM plasma from the time of gravitational wave production
to the present day.999 Despite the dearth of direct gravitational-wave probes
at frequencies between $10^{-15}$ and $10^{-9}$ Hz, two indirect constraints
apply around this range. Precision CMB measurements of the energy density in
radiation in the Universe provide an upper bound on the present-day
gravitational wave spectrum (gravitational waves themselves contributing to
expansion like radiation) [68, 69], currently of order
$\Omega_{\mathrm{GW},0}h^{2}\sim 10^{-6}$ [70, 71]. In addition, recent work
has argued that spectral distortions of the CMB blackbody also probe
gravitational waves in this frequency range [72]. While far-future experiments
could provide tighter constraints than $N_{\mathrm{eff}}$ measurements, these
still are unlikely to be useful probes of friendly axions.  Note that the
stochastic background from single-field inflation, if detectable by future CMB
experiments, is nearly scale invariant and of order $10^{-16}$ [73, 74], far
larger than those possible from friendly axion DM.

On the other hand, existing measurements of the $B$-mode polarization of the
CMB already constrain gravitational waves at frequencies between $10^{-18}$
and $10^{-16}$ Hz, with a most stringent upper limit of
$\Omega_{\mathrm{GW},0}h^{2}\sim 10^{-16}$ at $f_{\mathrm{GW}}\sim
10^{-17}\,\mathrm{Hz}$ [75]. Importantly, the polarization is sourced almost
exclusively at recombination, when the photon visibility function spikes. As a
result, CMB constraints are only relevant for scenarios where gravitational
waves are sourced before this time, i.e., when the Hubble scale is $H\gtrsim
H_{\mathrm{rec}}\sim 3\times 10^{-29}\,\mathrm{eV}$ [70]. In the scenarios we
consider here, the anisotropic stress maximizes after about ten field
oscillations, so the smallest relevant mass is of order
$10^{-27}\,\mathrm{eV}$. Consulting Eq. 15, the peak of the signal at such a
mass corresponds to a frequency of order $2\times 10^{-16}\,\mathrm{Hz}$. Such
hyperlight axion dark matter is well ruled out by fuzzy dark matter
constraints, but could make up some subcomponent of the total dark matter
(depending on the mass) [76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86].

By the preceding argument, CMB measurements mainly probe the infrared tail of
the gravitational wave background from friendly axions. Simulating a large
enough volume to capture these modes while also resolving the nonlinear
dynamics at small scales requires orders of magnitude more computational
resources than available to us. Fortunately, on causal grounds the
gravitational wave spectrum on scales much larger than the relevant dynamical
scales (i.e., $k/a\ll m$) follows a simple power law, independent of the
underlying dynamics [87]. We therefore extrapolate the signals computed in
simulations as decaying with $f_{\mathrm{GW}}^{3}$ at smaller frequencies,
appropriate for infrared, “causal” modes generated inside the horizon and
those generated while superhorizon that reenter during radiation-dominated
era. The lowest-frequency modes in the simulation do nearly follow an
$f_{\mathrm{GW}}^{3}$ scaling, so we expect this approximation to be at worst
conservative. Constraints were similarly derived on the infrared tail of
gravitational waves from a model of early dark energy in Ref. [88]. However,
note that recombination occurs shortly after matter-radiation equality, and
superhorizon causal modes that reenter the horizon during the matter era
instead grow with one power of $f_{\mathrm{GW}}$.101010 See Ref. [87] for a
thorough presentation of the imprint of the expansion history and presence of
free-streaming radiation on causal gravitational waves.  Since the CMB becomes
increasingly less sensitive to $\Omega_{\mathrm{GW}}$ on scales larger than
the horizon at equality, we do not expect the break in the power law to
improve constraints. (On the other hand, it would likely be observable for
causal gravitational wave backgrounds that are indeed observable in the CMB.)
Our simulations themselves also do not account for the transition to matter
domination, instead assuming a radiation-dominated Universe; we expect this to
be entirely sufficient for our estimates here.

We investigate the possibility of CMB constraints in Fig. 5, which displays
the possible signals from friendly axions subcomponents with $\mathcal{F}=20$.

Figure 5:  Gravitational wave backgrounds from friendly axions with fixed
decay constant ratio $\mathcal{F}=20$ and mass ratio $\mu=0.75$. Each curve
depicts the signal for a short axion mass making up a fraction of the total
dark matter $F_{\mathrm{DM}}$ as allowed by CMB and large scale structure data
[83], indicated by the legend. Dashed lines indicate the extrapolation of the
signal computed in the simulations (solid lines) to smaller frequencies as an
$f_{\mathrm{GW}}^{3}$ power law as justified in the text. The transparent gray
curves, with masses for which the aforementioned probes do not provide
constraints, take the axions to make up all of dark matter. Such light dark
matter is ruled out by fuzzy dark matter constraints, but we include these
curves to illustrate the mass dependence of the signals. CMB constraints from
Ref. [75] are superimposed in light green.

To illustrate the constraining power on the present abundance of friendly
axions, we consider masses varying from $10^{-27}$ to $10^{-25}\,\mathrm{eV}$,
each taking the fraction of dark matter in friendly axions to saturate limits
on ultralight axions from CMB and large scale structure data [83]. Because the
CMB is highly sensitive to low-frequency tensor modes (in terms of their
effective energy density) [75], the CMB can place tight constraints on the
fraction of dark matter in hyperlight friendly axions. For
$m=10^{-27}\,\mathrm{eV}$, Planck and BICEP2/Keck [89, 90] allow only a
$0.1\%$ friendly subcomponent.111111 Note that the most recent BICEP/Keck data
release further improved constraints on the tensor-to-scalar ratio by a factor
of $\sim 1.6$, and that CMB-S4 projects to provide upwards of a further factor
of $\sim 30$ [91].  We may obtain a heuristic bound as a function of mass by
extrapolating the $f_{\mathrm{GW}}^{3}$ tail to the frequency corresponding to
the horizon size at matter-radiation equality, $f_{\mathrm{eq}}=1.54\times
10^{-15}\,\mathrm{Hz}$, where Ref. [75] reports an upper limit of
$\Omega_{\mathrm{GW},0}^{\mathrm{bound}}h^{2}=3.2\times 10^{-16}$. Applying
the scalings of Eqs. 14 and 15, the CMB probes

$\displaystyle\frac{\Omega_{h}+\Omega_{l}}{\Omega_{\mathrm{DM}}}\lesssim
10^{-2}\left(\frac{m}{5\times
10^{-27}\,\mathrm{eV}}\right)^{5/4}\sqrt{\frac{\Omega_{\mathrm{GW},0}^{\mathrm{bound}}(f_{\mathrm{eq}})h^{2}}{3.2\times
10^{-16}}}$ (16)

which is only relevant for masses $m\gtrsim 10^{-27}\,\mathrm{eV}$ for which
the signals are produced before recombination. With the bound of Ref. [75],
the limits become irrelevant (i.e., order unity) for masses $m\gtrsim
10^{-25}\,\mathrm{eV}$, well below the present lower limits on the mass of
fuzzy dark matter.

Finally, for illustrative purposes Fig. 5 includes signals for friendly axions
with larger masses $m>10^{-25}\,\mathrm{eV}$ taken to make up all of the dark
matter. While these scenarios themselves are ruled out by lower limits on the
mass of axion dark matter, they demonstrate that gravitational waves from
friendly axions would likely be unobservably small at allowed masses. Even for
$m=10^{-20}\,\mathrm{eV}$ the peak of the signal is only $\sim 10^{-15}$,
consistent with the preceding discussion.

###  III.3 Driven oscillons

Beyond their important role in the nonperturbative dynamics of friendly
axions, oscillons have long been a subject of interest [22, 23, 25, 26, 27,
28]. As nontopological excitations, oscillons generically decay by radiating
semirelativistic modes, making their lifetime an interesting object of study
[29, 30, 31]. In models with friendly axions, the lifetime of an oscillon can
be extended beyond its in-vacuum expectation because of energy transfer from
the long axion [18]. In this section, we show that short axion oscillons are
driven by the long axion via autoresonance, obeying equations analogous to
those for the homogeneous fields.121212 Other related nonlinear wave equations
exhibit oscillons/solitons driven by autoresonance, e.g., Refs. [92, 93, 94,
95, 96].  We then provide analytic results to quantify the oscillon lifetime
enhancement and numerical results to support our analysis. Because of the
importance of nonlinear dynamics, it is more convenient to discuss oscillons
in the interaction basis of the short and long axion, which we adopt
consistently throughout this section.

An oscillon is a quasiperiodic, quasilocalized excitation of the axion field,
and its fundamental frequency $\omega$ is smaller than the rest mass of the
axion $m$. After its initial formation, the binding energy per particle inside
an isolated oscillon,

$\displaystyle E_{\text{bind}}=m-\omega,$ (17)

is a decreasing function of time, and consequently the characteristic size of
the oscillon

$\displaystyle r_{\text{osc}}\approx\frac{\pi}{\sqrt{2mE_{\text{bind}}}},$
(18)

is an increasing function of time. Over time its rate of radiation falls off
approximately exponentially with this growing separation of scales [31]:

$\displaystyle
P_{\text{rad}}/f^{2}\sim\exp\left[-\frac{r_{\text{osc}}}{\lambda_{\text{rad}}}\right]\sim\exp\left[-\frac{3\omega}{\sqrt{2m(m-\omega)}}\right],$
(19)

where $\lambda_{\mathrm{rad}}=3\omega$ is the radiation wavelength due to
$3\to 1$ processes, which dominate the decay for most $\omega$ in potentials
with parity symmetry. Thus as oscillons age and $\omega$ increases toward $m$,
the oscillon radiates more and more slowly. However, the oscillon is not
infinitely long-lived: a geometrical consequence of living in three dimensions
is that an oscillon must die at a frequency $\omega_{\text{crit}}<m$ (see
e.g., Refs. [97, 98, 31]). An oscillon therefore typically spends the majority
of its lifetime at a frequency close to $\omega_{\text{crit}}$.

The mechanism of radiation (which is discussed in detail in Refs. [98, 30,
31]) does not play an important role in our analysis, so we may simply
characterize the decay rate of the oscillon by its instantaneous lifetime:

$\displaystyle T_{\text{inst}}(\omega)$
$\displaystyle=\frac{E_{\text{osc}}(\omega)}{P_{\text{rad}}(\omega)}\equiv\frac{1}{\Gamma_{\text{inst}}(\omega)},$
(20)

where $E_{\text{osc}}$ is the total bound energy in the oscillon and
$P_{\text{rad}}$ is the power radiated by the oscillon, all measured while the
oscillon is at a frequency $\omega$. [The precise values of
$E_{\mathrm{osc}}(\omega)$ and $P_{\mathrm{rad}}(\omega)$ may be calculated
using the software package  Simple Oscillon.]

As shown in earlier sections, friendly axions admit oscillon solutions too,
although only the short axion is likely to form them. One may gain insight
into the dynamics of short oscillons by taking a spherically symmetric ansatz
with a fixed radial profile, $\phi_{S}(t,r)=\Phi_{S}(t)R_{S}(r)$, in the
background of a homogeneous long axion, $\phi_{L}(t,r)=\Phi_{L}(t)$.131313
Note that $\phi_{S}(t,r)$ is in general not separable, since as the oscillon
radiates energy it adjusts its radial profile, but neglecting such changes is
sufficient on timescales much shorter than the oscillon lifetime. Moreover, as
demonstrated below, the long axion supplies exactly enough energy to drive the
short oscillon at a fixed amplitude via autoresonance, thus maintaining its
radial profile.  To be self-consistent, we work in the limit where the short
axion does not backreact onto the long axion, i.e., $f_{S}\ll f_{L}$. This
approximation neglects the effects of radiation, which we include via an
artificial damping term with coefficient $\Gamma_{\mathrm{inst}}$. In the
small $\Phi_{L}/f_{L}$ limit, keeping only linear terms,

$\displaystyle\ddot{\Phi}_{L}+\frac{3}{2t}\dot{\Phi}_{L}$
$\displaystyle=-m^{2}\mu^{2}\Phi_{L},$ (21a)
$\displaystyle\ddot{\Phi}_{S}+\left(\frac{3}{2t}+\Gamma_{\text{inst}}\right)\dot{\Phi}_{S}$
$\displaystyle=-V_{0}^{\prime}(\Phi_{S})-V_{1}^{\prime}(\Phi_{S})\frac{\Phi_{L}}{f_{L}}.$
(21b)

Here $V_{0}$ and $V_{1}$ are effective potentials derived by integrating out
$R_{S}$.141414 To integrate out the radial profile, one substitutes
$\phi_{S}=\Phi_{S}(t)R_{S}(r)$ into the action and integrates over space. The
qualitative features of the resulting effective action for $\Phi_{S}$ are
insensitive to the choice of $R_{S}(r)$ so long as it is monotonically
decreasing and changes over a radial length scale of at least $1/m$.
Ultimately, the precise form of these potentials is unimportant, their salient
feature being that they possess attractive nonlinearities. These equations are
thus precisely in the same form as those for the homogeneous system studied in
Ref. [18], and autoresonance between $\Phi_{S}$ and $\Phi_{L}$ therefore
occurs for a broad range of initial conditions (although the likelihood that
any given patch of space leads to oscillon formation and autoresonance is most
easily assessed with simulations).

Figure 6 demonstrates that autoresonance between the short oscillons and the
long axion occurs in our full $3+1D$ simulations.

Figure 6:  Central amplitude of a short-axion oscillon (left vertical axis,
blue curves) and average amplitude of the long axion (right vertical axis, red
curves), measured in a cosmological simulation with $N^{3}=512^{3}$
gridpoints, box length $L=\pi/m$, mass ratio $\mu=0.75$, and decay constant
ratio $\mathcal{F}=50$. Note the break in the horizontal axis between
$520\lesssim mt\lesssim 850$.

The short-axion oscillon oscillates at the same frequency $\omega=\mu m$ as
the driver $\left\langle\phi_{L}\right\rangle$, conclusively demonstrating
local nonlinear resonance inside the oscillon. The phase offset between the
short and long axion evolves from roughly $\pi/3$ at early times (left side of
Fig. 6) to about $\pi/2$ just before oscillon death (right side), indicating
an increasingly efficient energy transfer from the long axion to the short
oscillon. The energy transfer rate peaks when the phase shift reaches $\pi/2$,
and subsequently decreases as the long axion’s amplitude redshifts. In analogy
to homogeneous autoresonance, the long axion can then no longer support the
short axion against its own radiation.

Using the equations of motion Eq. 21 for the spherically symmetric system, we
can solve for the dynamics of the short oscillon and the long driver to arrive
at the driven oscillon lifetime $t_{\text{death}}$ in terms of its
instantaneous vacuum lifetime $T_{\text{inst}}(\mu)$. The details of this
calculation are described in App. B, and we summarize the results here. First,
the long driver amplitude must be large enough that it supplies sufficient
energy to the oscillon, leading to

$\displaystyle m\mu t_{\text{death}}\approx\left[m\mu
T_{\text{inst}}(\mu)\right]^{4/3}.$ (22)

We compare this analytic scaling to simulations in Fig. 7, verifying the $[\mu
T_{\mathrm{inst}}(\mu)]^{4/3}$ dependence in the range $0.73\lesssim\mu\leq
0.81$. but spherically symmetric simulations verify the scaling of Eq. 22 out
to $\mu=0.89$.

Figure 7:  Driven oscillon lifetime measured in $3+1D$ cosmological
simulations (red dots) versus the analytic prediction of Eq. 22 (blue line) at
a decay constant ratio $\mathcal{F}=1000$ (large enough that backreaction is
never important). The analytic curve only contains one free parameter
corresponding to an order 1 multiplicative constant in Eq. 22 (see App. B).
Note that the driver frequency $\omega_{\mathrm{driver}}$ simply corresponds
to the mass of the long axion, $\mu m$.

At lower frequencies the driven oscillon lifetime is shorter than the vacuum
lifetime: once the oscillon falls off autoresonance, it rapidly dumps its
energy into the long axion field, cutting its life short. Larger values of
$\mu$ require longer ($3+1D$) simulation runtimes than our computational
resources permit.

The oscillon also backreacts on the driver, inducing spatial perturbations.
These fluctuations remain small only until

$\displaystyle m\mu t_{\text{death}}\lesssim\mathcal{F}^{8/3}.$ (23)

This scaling is demonstrated in Fig. 8: the spherically symmetric solutions to
the full coupled nonlinear wave equations (see Eq. 48) exhibit
$\mathcal{F}^{8/3}$ behavior for $\mathcal{F}\lesssim 40$, at which point the
lifetime saturates the driving limit, Eq. 22.

Figure 8:  Driven oscillon lifetime versus decay constant ratio $\mathcal{F}$
in spherically symmetric simulations with $\mu=0.85$ (red) with superposed
$\mathcal{F}^{8/3}$ scaling predicted by Eq. 23 (blue). The wiggles in the red
curve are due to the precise phase relationship between the oscillon and the
driver, a feature that would be washed out in any realistic simulation with
randomly sampled initial phases.

Performing a quantitative analysis of the parametric dependence on
$\mathcal{F}$ in $3+1D$ simulations would require simulations with larger mass
ratio $\mu$, which, as discussed previously, exceed the limits of our
computing resources. However, we do observe qualitative behavior—lifetimes
that increase with $\mathcal{F}$ and saturate at the driving limit in Fig.
7—similar to that of Fig. 8 in the set of simulations presented above.

Finally, the oscillon siphons energy from the long axion, depleting the
latter’s local energy density. Nearby regions of space then resupply this
region with long axion; requiring that this resupply rate exceeds the
depletion rate due to the oscillon leads to the final constraint

$\displaystyle t_{\text{death}}\lesssim\mathcal{F}^{2}T_{\text{inst}}(\mu).$
(24)

One can check that there is no region of $\mathcal{F}-T_{\mathrm{inst}}$
parameter space where Eq. 24 dominates the lifetime. Taken together, the
lifetime is

$\displaystyle m\mu t_{\text{death}}$
$\displaystyle\approx\min\left(\left[m\mu
T_{\text{inst}}(\mu)\right]^{4/3},\mathcal{F}^{8/3}\right).$ (25)

While the lifetimes of these driven oscillons are parametrically enhanced
relative to their in-vacuum expected lifespans, they are still far too short-
lived to be of any cosmological relevance. Nonetheless, these novel
dynamics—interesting in their own right—potentially broaden the class of
scalar field theories that admit oscillons surviving into the present day. We
discuss this possibility and associated challenges in App. B.2.

##  IV Discussion

Nonlinear effects in the early Universe can have a drastic impact on the late-
time distribution of energy in dark sectors. The “friendly” axion system of
Ref. [18] provides a concrete example model, where nonlinearities dominate the
dynamics of both the background and fluctuations and have important
consequences for direct detection experiments. In this paper, we numerically
evolved the full system, showing that large axion fluctuations—in particular
the nucleation of oscillons—work to equilibrate the relic densities of the two
axions for a moderate ratio of the heavy axion’s decay constant to that of the
light one ($6\lesssim\mathcal{F}\lesssim 20$). For smaller decay constant
ratios $\mathcal{F}\lesssim 6$, spatial fluctuations have a negligible effect
on the dynamics, and we recover the results of homogeneous computations from
Ref. [18]. At larger ratios $\mathcal{F}\gg 20$, oscillon nucleation prevents
the heavy axion from ever attaining a substantial abundance.

The novel dynamics in the intermediate-$\mathcal{F}$ regime position friendly
axions to be positively identified as two-component dark matter by direct
detection experiments. The lighter axion’s abundance is reduced by no more
than a factor of about two, in sharp contrast to expectations based on
homogeneous approximations in which its abundance would be parametrically
depleted [18]. The heavier axion’s abundance (and therefore detection
prospects) is still parametrically enhanced (by a factor of
$\approx\mathcal{F}^{2}/2$), but only at a moderate cost to the visibility of
the lighter axion. Many upcoming axion direct detection experiments [43, 44,
45, 50] would potentially be sensitive to both axions in a friendly pair
having masses within the experiment’s sensitivity band. Direct detection of
axion dark matter with a decay constant substantially smaller than that
expected in standard misalignment scenarios should prompt a search for a
second, more weakly coupled axion at a nearby mass.

We also computed the stochastic gravitational wave background produced by
oscillon nucleation in a friendly axion sector. If friendly axions compose all
of the dark matter, the present-day strain is well out of reach of near-future
gravitational wave experiments, but the cosmic microwave background
polarization does constrain (and in the future may probe) hyperlight friendly
pairs making up a subcomponent of dark matter. Density and vector
perturbations are also produced in these scenarios; their effect on the CMB
(and other cosmological observables) is less straightforward to evaluate, but
they may well provide even more stringent constraints than just the (as-yet
unobserved) primordial $B$-mode polarization.

Finally, for $\mathcal{F}\gtrsim 20$, although autoresonance is quenched by
oscillon production (preventing the axions’ energies from equalizing) our
simulations demonstrate that short-axion oscillons produced in the early
Universe are driven by the long axion background, parametrically extending
their lifetimes. For the specific friendly axion potential studied here (Eq.
3), driven oscillons can live about an order of magnitude longer than their
undriven counterparts. Though they are still not long-lived enough to be
astrophysically relevant, even at the lightest possible axion masses, this may
not be the case for other scalar potentials. Similar dynamics in other coupled
axion theories may lead to driven oscillons that could naturally live until
the present day, with numerous possible observational signatures including
gravitational lensing [99], optical lensing [100], and electromagnetic bursts
[101, 102, 103, 104].

Our results are qualitatively insensitive to the amplitude of the initial
primordial curvature perturbations. However, the size does determine the
precise minimum and maximum decay constant ratios for which the two axion
energy densities equalize. For simplicity, we used a scale-invariant initial
power spectrum with magnitude set by the Planck measurements at CMB scales,
but in reality adiabatic fluctuations are red-tilted on large scales and are
much less constrained on smaller scales. If there is less initial power at
small scales, the time to nonlinearity and oscillon formation increases,
allowing friendly pairs with larger decay constant ratios $\mathcal{F}\gg 20$
to achieve equipartition. We also note that whether the initial axion
perturbations are adiabatic as above or seeded directly in the axion field
(i.e., isocurvature perturbations) does not measurably affect any of our
results, as corroborated by simulations with purely isocurvature initial
conditions.

The string axiverse is a rare example of a low-energy signature of quantum
gravity, most of whose novel predictions reside at the grand-unified or string
scales, far outside experimental reach. In general, an axiverse can comprise a
multitude of light, coupled axions; this work provides further evidence of the
outsized impact nonlinear effects and interactions have on the phenomenology
of the axiverse. The friendly model considered here, for which nonperturbative
dynamics revise predictions by multiple orders of magnitude, is only a
prototypical example; further work is necessary to understand the
phenomenology of fully realistic axiverses and the critical role played by
nonlinear dynamics.

###### Acknowledgements.

We thank Masha Baryakhtar, Davide Racco, and the anonymous referee for
thoughtful feedback on this manuscript, Savas Dimopoulos for useful
discussions on “oscillon longevity centers,” and Dmitriy Zhigunov for helpful
conversations about strongly nonlinear fields. D.C. is grateful for the
support of the Stanford Institute for Theoretical Physics (SITP), the National
Science Foundation under Grant No. PHY-2014215, and the Gordon and Betty Moore
Foundation under Grant No. GBMF7946. O.S. is supported by a DARE fellowship
from the Vice Provost for Graduate Education at Stanford University. J.O.T. is
supported by the ARCS Foundation, and is thankful to the Perimeter Institute
for its hospitality during the final stages of completing this manuscript.
Z.J.W. is supported by the Department of Physics and the College of Arts and
Sciences at the University of Washington. This work used the Extreme Science
and Engineering Discovery Environment (XSEDE) [105], which is supported by
National Science Foundation Grant No. ACI-1548562; simulations were run
through allocation No. TG-PHY200037 on the Anvil cluster at Purdue University,
Bridges-2 at the Pittsburgh Supercomputing Center which is supported by NSF
Grant No. ACI-1928147, and Expanse at the San Diego Supercomputer Center. This
research was supported in part by Perimeter Institute for Theoretical Physics.
Research at Perimeter Institute is supported by the Government of Canada
through the Department of Innovation, Science and Economic Development and by
the Province of Ontario through the Ministry of Research, Innovation and
Science. Simulations in this work were implemented with pystella [42], which
is available at github.com/zachjweiner/pystella and makes use of the Python
packages PyOpenCL [106], Loopy [107], mpi4py [108, 109], mpi4py-fft [110], and
NumPy [111]. This work also made use of the packages SciPy [112], matplotlib
[113], SymPy [114], and CMasher [115].

##  Appendix A Equations of motion and numerical implementation

For completeness, we enumerate the evolution equations as implemented in
simulations. We use a conformal-time, “mostly plus” FLRW metric in the
conformal Newtonian gauge with line element

$\displaystyle\begin{split}\mathrm{d}s^{2}&=-a(\tau)^{2}\left[1+2\Phi(\tau,\mathbf{x})\right]\mathrm{d}\tau^{2}\\\
&\hphantom{{}={}}+a(\tau)^{2}\left[\left\\{1-2\Phi(\tau,\mathbf{x})\right\\}\delta_{ij}+h_{ij}(\tau,\mathbf{x})\right]\mathrm{d}x^{i}\mathrm{d}x^{j}.\end{split}$
(26)

Here $\Phi$ is the Newtonian potential and $h_{ij}$ the transverse
($\partial_{i}h_{ij}=0$) and traceless ($\delta^{ij}h_{ij}=0$) tensor
perturbation. We neglect scalar anisotropic stress and vector perturbations.
We define the comoving Hubble parameter $\mathcal{H}\equiv\partial_{\tau}a/a$
(while the standard Hubble parameter is
$H\equiv\partial_{t}a/a=\mathcal{H}/a$).

Well before matter-radiation equality, the axions make a negligible
contribution to the expansion of the Universe; the solution to the Friedmann
equations in a radiation Universe is $a(\tau)/a(\tau_{m})=\tau/\tau_{m}$. We
take $\tau_{m}=m^{-1}$ so that the scale factor is normalized to the time when
$H=m$. Because tensor perturbations (i.e., gravitational waves) from inflation
are very small and their subsequent production by axion production is
suppressed by $(f_{L}/M_{\mathrm{pl}})^{2}$, we neglect their backreaction
onto the axion fields. For the same reason, we neglect the contribution of the
axions to the Newtonian potential.

###  A.1 Equations of motion

The Euler-Lagrange equation for the axions reads

$\displaystyle\begin{split}\partial_{\tau}^{2}\phi_{I}&=-\left[2\mathcal{H}(\tau)-4\partial_{\tau}\Phi\right]\partial_{\tau}\phi_{I}+\left(1+4\Phi\right)\partial_{i}\partial_{i}\phi_{I}\\\
&\hphantom{{}={}}-a(\tau)^{2}\left(1+2\Phi\right)\frac{\partial
V}{\partial\phi_{I}},\end{split}$ (27)

where the potential $V$ is defined in Eq. 2 and $I=S$ or $L$. (We use repeated
Latin indices $i$, $j$, $k$, etc., to denote spatial components that are
contracted with the Kronecker delta regardless of their placement.) In a
Universe dominated by a single fluid with equation of state $w$ and a sound
speed $c_{s}^{2}\equiv\delta P/\delta\rho$, the Einstein equations for $\Phi$
may be rearranged into

$\displaystyle\begin{split}\partial_{\tau}^{2}\Phi&=-(2+3c_{s}^{2})\mathcal{H}\partial_{\tau}\Phi-\mathcal{H}\partial_{\tau}\Phi\\\
&\hphantom{{}={}}-3\left(c_{s}^{2}-w\right)\mathcal{H}^{2}\Phi+c_{s}^{2}\partial_{i}\partial_{i}\Phi,\end{split}$
(28)

where both $w$ and $c_{s}^{2}$ are $1/3$ in the radiation-dominated era.
Finally, the tensor perturbations evolve according to

$\displaystyle\partial_{\tau}^{2}h_{ij}+2\mathcal{H}\partial_{\tau}h_{ij}-\partial_{k}\partial_{k}h_{ij}$
$\displaystyle=\frac{2a(\tau)^{2}}{M_{\mathrm{pl}}^{2}}\Pi^{i}_{\hphantom{i}j}.$
(29)

Gravitational waves are sourced by the transverse and traceless anisotropic
stress tensor $\Pi_{ij}$ whose Fourier modes are given in terms of the full
stress tensor $T_{ij}$ by

$\displaystyle\Pi_{ij}(\tau,\mathbf{k})$
$\displaystyle=\left(P_{il}(\mathbf{k})P_{jm}(\mathbf{k})-\frac{1}{2}P_{ij}(\mathbf{k})P_{lm}(\mathbf{k})\right)T_{lm}(\tau,\mathbf{k})$
(30)

with

$\displaystyle P_{ij}(\mathbf{k})$
$\displaystyle=\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}.$ (31)

Unlike the Newtonian potential, which is dominated by the standard adiabatic
perturbations in the radiation fluid, the gravitational waves are only sourced
by the axions. However, the sourced gravitational waves are also suppressed by
$(f_{L}/M_{\mathrm{pl}})^{2}$, for which reason it is entirely sufficient to
evaluate the stress tensor as if in a homogeneous FLRW Universe:

$\displaystyle\begin{split}T^{i}_{\hphantom{i}j}&=\sum_{I}\partial_{i}\phi_{I}\partial_{j}\phi_{I}\\\
&\hphantom{{}={}}-\delta^{i}_{\hphantom{i}j}\left(\frac{1}{2}\sum_{I}\partial_{\mu}\phi_{I}\partial^{\mu}\phi_{I}+V(\phi_{S},\phi_{L})\right).\end{split}$
(32)

The final term purely contributes to the trace of the stress-energy tensor and
may be dropped when computing the metric tensor perturbations. Because the
backreaction of gravitational waves on the axions is negligible, we may simply
ignore any initial amplitude generated during inflation and consider the
evolution of Eq. 29 from zero initial conditions.

###  A.2 Initial conditions

Imposing that $\Phi$ was frozen outside of the horizon ($k\tau\ll 1$) to its
primordial value generated during inflation sets $\Phi(\tau\ll
1/k,\mathbf{k})=\Phi_{0}(\mathbf{k})$ and $\partial_{\tau}\Phi(\tau\ll
1/k,\mathbf{k})=0$. In this case, we find the solution

$\displaystyle\Phi(\tau,\mathbf{k})$
$\displaystyle=\Phi_{0}(\mathbf{k})\frac{\sin y-y\cos y}{y^{3}},$ (33)

where $y\equiv\sqrt{w}k\tau$. The primordial curvature perturbation is
characterized by a dimensionless power spectrum $\Delta_{\Phi}^{2}(k)$ defined
by

$\displaystyle\frac{k^{3}}{2\pi^{2}}\left\langle\Phi_{0}(\mathbf{k}_{1})\Phi_{0}(\mathbf{k}_{2})\right\rangle$
$\displaystyle=(2\pi)^{3}\delta(\mathbf{k}_{1}+\mathbf{k}_{2})\Delta_{\Phi}^{2}(k),$
(34)

which, for standard slow-roll inflation, is nearly scale invariant with
amplitude $\sim 10^{-9}$ [70]. To avoid specifying the mass scale of the
problem, we neglect the small spectral tilt and simply take a scale-invariant
spectrum.

We begin the simulations when $H=m$ (at $\tau=\tau_{1}$), solving the
linearized system of equations from an early time when all relevant wave
numbers are well outside the horizon (i.e., when $k\ll aH$) to determine
initial conditions. The axions are initialized as Gaussian-random fields with
mean field value and velocity at $\tau_{1}$ set according to the solution to
the homogeneous system. The fluctuations are set in Fourier space to match the
linearized solution $\phi_{I,k}(\tau_{1})$, multiplied by a Gaussian-random
complex number (normalized such that the mean squared modulus is unity). For
each axion and wave vector $\mathbf{k}$ on the grid, we therefore set

$\displaystyle\phi_{I}(\tau_{1},\mathbf{k})$
$\displaystyle=\phi_{I,k}(\tau_{1})\sqrt{-\ln U_{1}(\mathbf{k})}e^{2\pi
iU_{2}(\mathbf{k})}$ (35)
$\displaystyle\partial_{\tau}\phi_{I}(\tau_{1},\mathbf{k})$
$\displaystyle=\partial_{\tau}\phi_{I,k}(\tau_{1})\sqrt{-\ln
U_{1}(\mathbf{k})}e^{2\pi iU_{2}(\mathbf{k})},$ (36)

where $U_{1}(\mathbf{k})$ and $U_{2}(\mathbf{k})$ are random variates between
0 and 1, each of which are the same for both axions. The Newtonian potential
is initialized analogously (and with random numbers for each wave vector
matching those for the axions) using the analytic solution Eq. 33.

We separately implemented simulations that neglect metric perturbations,
instead initializing the axion with a scale-invariant spectrum of isocurvature
fluctuations (with amplitude comparable to that of the adiabatic initial
conditions). At the level of the analysis performed in Sec. III, the results
are effectively unchanged, since the initial size and early evolution of
fluctuations are important only insofar as they affect the time to
nonlinearity. Nonlinear interactions quench autoresonance and distribute power
across length scales, largely washing out any detailed features of the initial
conditions. Since amplification (due to parametric resonance) does not begin
immediately after modes start oscillating, quantitative results retain only a
small, logarithmic dependence on initial conditions.

###  A.3 Gravitational wave backgrounds

Gravitational waves carry an effective energy density which, deep inside the
horizon (i.e., $k\gg\mathcal{H}$), is [116, 117, 75]

$\displaystyle\rho_{\mathrm{GW}}(\tau)=\frac{M_{\mathrm{pl}}^{2}}{4a^{2}}\left\langle\mkern
1.25mu\overline{\mkern-1.25mu\partial_{\tau}h_{ij}\partial_{\tau}h_{ij}\mkern-1.25mu}\mkern
1.25mu\right\rangle,$ (37)

where the bar denotes a time average (i.e., over oscillations). The relic
abundance of gravitational waves today per logarithmic wave number is

$\displaystyle\Omega_{\mathrm{GW}}(\tau,k)=\frac{1}{\bar{\rho}(\tau)}\frac{\mathrm{d}{\rho_{\mathrm{GW}}(\tau)}}{\mathrm{d}{\ln
k}}.$ (38)

Substituting the inverse Fourier transform of $h_{ij}$ into Eq. 37 permits
rewriting Eq. 38 in terms of the dimensionless power spectrum of
$\partial_{\tau}h_{ij}$ (defined in analogy to Eq. 34) as

$\displaystyle\Omega_{\mathrm{GW}}(\tau,k)=\frac{1}{12\mathcal{H}(\tau)^{2}}\sum_{i,j}\Delta^{2}_{\partial_{\tau}h_{ij}}(\tau,k),$
(39)

after plugging in $\bar{\rho}(\tau)=3M_{\mathrm{pl}}^{2}H(\tau)^{2}$. The
spectrum evaluated in the early Universe is related to that at the present day
(at $\tau_{0}$) by the transfer function

$\displaystyle\frac{\Omega_{\mathrm{GW}}(\tau_{0},k)h^{2}}{\Omega_{\mathrm{GW}}(\tau,k)}$
$\displaystyle=\Omega_{\mathrm{rad}}(\tau_{0})h^{2}\frac{g_{\star}(\tau)}{g_{\star}(\tau_{0})}\left(\frac{g_{\star
S}(\tau)}{g_{\star S}(\tau_{0})}\right)^{-4/3},$ (40)

and would be observed at present-day frequencies related to the wave number
$k$ by

$\displaystyle\begin{split}f_{\mathrm{GW}}&=\frac{k/2\pi
a(\tau)}{\sqrt{H(\tau)M_{\mathrm{pl}}}}\left[\Omega_{\mathrm{rad}}(\tau_{0})H_{0}^{2}M_{\mathrm{pl}}^{2}\right]^{1/4}\\\
&\hphantom{{}={}}\times\left(\frac{g_{\star}(\tau)}{g_{\star}(\tau_{0})}\right)^{1/4}\left(\frac{g_{\star
S}(\tau)}{g_{\star S}(\tau_{0})}\right)^{-1/3}.\end{split}$ (41)

Here $g_{\star}$ and $g_{\star S}$ are the numbers of relativistic degrees of
freedom in energy and entropy density, respectively. Note that the present-day
abundance of radiation is $\Omega_{\mathrm{rad}}(\tau_{0})h^{2}\approx
4.2\times 10^{-5}$ [70] and that $H_{0}/h\equiv
100\,\mathrm{km}\,\mathrm{s}^{-1}/\mathrm{Mpc}\approx 3.24\times
10^{-18}\,\mathrm{Hz}$.

###  A.4 Numerical implementation

We discretize the evolution equations, Eqs. 27, 28 and 29, onto a three
dimensional, regularly spaced grid with periodic boundary conditions.
Following Refs. [41, 40, 88], we evolve the gravitational wave equation of
motion (Eq. 29) under the replacement of the transverse-traceless stress
tensor (which is only easily calculated in Fourier space) with the full stress
tensor (Eq. 32). The transverse-traceless projection is instead performed on
the tensor field itself only when outputting gravitational wave spectra,
drastically reducing the required number of fast Fourier Transforms (which in
distributed-memory contexts are a major bottleneck). For similar reasons,
rather than use the analytic solution Eq. 33 for the Newtonian potential
(which requires forward and inverse Fourier transforms at each step), we
simply evolve Eq. 28 in position space.

Spatial derivatives are approximated with fourth-order centered differencing,
while integration in time is implemented with a “low-storage,” fourth-order
Runge-Kutta method [118]. We have verified that the results are insensitive to
the precise choice of box length $L$ and are consistent with simulations with
the same physical volume but $1280^{3}$ and $1536^{3}$ gridpoints (compared to
the $1024^{3}$ used for main results), as well as at the same resolution with
a box length $1.25$ and $1.5$ times larger. We have also checked for a
representative case that results at a fixed resolution and volume are
qualitatively insensitive to the the random seed used in generating stochastic
initial conditions.

##  Appendix B Details of driven oscillons

In Sec. III.3, we demonstrated that short-axion oscillons remain in
autoresonance with the nearly homogeneous long axion, often extending their
lifetime well beyond the in-vacuum expectation. The physics of this localized
autoresonance is mostly captured by the effective equations of motion (Eq. 21)
obtained by integrating out the spatial profile of the oscillon, demonstrating
that the short oscillon is well approximated as a one dimensional driven
nonlinear oscillator. However, Eq. 21 fails to capture some important effects
that ultimately limit the lifetime of the driven oscillon, namely Eqs. 22, 23
and 24. In what follows, we derive the lifetime bounds in Eqs. 22, 23 and 24
in Sec. B.1 and discuss the possibility of cosmologically long-lived oscillons
in more general multiaxion potentials in Sec. B.2.

###  B.1 Lifetime bounds

Most long-lived oscillons are well approximated by a single-harmonic ansatz,
so $\Phi_{S}$ oscillates predominantly at a single frequency. On
autoresonance, this frequency approximately matches the long-axion mass ($\mu
m$) up to transient oscillations around this stable point. Thus on
autoresonance we may approximate

$\displaystyle\Phi_{S}$ $\displaystyle\approx f_{S}\Theta_{S}^{0}\sin(m\mu
t+\delta),$ (42) $\displaystyle\Phi_{L}$ $\displaystyle\approx
f_{L}\Theta_{L}^{0}\left(\frac{t}{t_{0}}\right)^{-3/4}\sin(m\mu t),$ (43)

where $\delta$ is a phase-offset. While $\delta$ does have a time dependence,
as discussed in App. B.3 of [18], it is slow compared to the oscillatory
timescale and thus we can approximate $\delta$ as a constant. The power
transferred from the long to the short axion is

$\displaystyle\frac{P_{L\to S}}{V_{\text{osc}}}$
$\displaystyle=\left(\dot{\Phi}_{S}\partial_{\Phi_{S}}-\dot{\Phi}_{L}\partial_{\Phi_{L}}\right)V_{\text{int}},$
(44)

where for our purposes we may approximate the interaction potential by a mass-
mixing term,

$\displaystyle V_{\text{int}}\sim m^{2}f_{S}\Phi_{S}\frac{\Phi_{L}}{f_{L}}.$
(45)

Thus, the time-averaged power transfer is

$\displaystyle\langle P_{L\to S}\rangle\approx-m\mu
E_{\text{osc}}\frac{\Theta_{L}^{0}}{\Theta_{S}^{0}}\left(\frac{t}{t_{0}}\right)^{-3/4}\sin\delta,$
(46)

where we’ve taken $m^{2}f_{S}^{2}(\Theta_{S}^{0})^{2}V_{\text{osc}}\approx
E_{\text{osc}}$. As shown in Ref. [18], the phase $\delta$ tends toward
$-\pi/2$ toward the end of autoresonance, and thus we take $\delta=-\pi/2$
when calculating the driven oscillon lifetime. The oscillon is supported by
the driver when its radiated power $P_{\mathrm{rad}}$ is smaller than the
maximum power transfer from the long axion $\langle P_{L\to S}\rangle$, from
which we obtain the approximate time of death $t_{\mathrm{death}}$,

$\displaystyle m\mu
t_{\text{death}}\approx\left(m\mu\frac{\Theta_{L}^{0}}{\Theta_{S}^{0}}T_{\text{inst}}(\mu)\right)^{4/3},$
(47)

taking $t_{0}\approx 1/m\mu$. Here
$T_{\mathrm{inst}}(\omega)=E_{\text{osc}}(\omega)/P_{\text{rad}}(\omega)$ is
the instantaneous oscillon lifetime at the frequency $\omega$, which may be
longer or shorter than the in-vacuum oscillon lifetime. Provided that the
driving frequency $m\mu$ corresponds to a slow-decaying part of the oscillon
lifecycle, this result shows that the driven oscillon lifetime is
parametrically enhanced relative to its vacuum lifetime.

So far, we have neglected the backreaction of $\theta_{S}$ onto $\theta_{L}$
in order to approximate $\theta_{L}$ as a homogeneous field. While the details
of backreaction are complicated, it is simple to obtain an estimate for when
$\theta_{S}$ causes ${\cal O}(1)$ fluctuations in $\theta_{L}$ which then
terminate the autoresonance. To make this estimate, we observe in the
equations of motion,

$\displaystyle\square\theta_{L}+m^{2}\mathcal{F}^{-2}\sin(\theta_{S}+\theta_{L})+m^{2}\mu^{2}\sin\theta_{L}$
$\displaystyle=0$ (48a)
$\displaystyle\square\theta_{S}+m^{2}\mu^{2}\sin(\theta_{S}+\theta_{L})$
$\displaystyle=0,$ (48b)

that $\theta_{L}$’s evolution depends on $\theta_{S}$ multiplied by
$\mathcal{F}^{-2}$. Comparing this term to the final term in Eq. 48a, we see
that backreaction becomes important when

$\displaystyle\theta_{L}\sim\mathcal{F}^{-2}\theta_{S}.$ (49)

Assuming that $\theta_{S}$ has a constant $\mathcal{O}(1)$ amplitude, as
inside an oscillon, and that $\theta_{L}$ decays like cold matter as in Eq.
43, we find that spatial fluctuations in $\theta_{L}$ become order one at the
time

$\displaystyle m\mu
t_{\text{death}}\sim\left(\frac{\Theta_{L}^{0}}{\Theta_{S}^{0}}\right)^{4/3}\mathcal{F}^{8/3}.$
(50)

After this time, the long axion no longer serves as a good driver and quickly
dephases with the short axion oscillon, siphoning its energy and causing its
rapid death. We plot this predicted $\mathcal{F}^{8/3}$ scaling versus the
observed lifetime in a spherically symmetric driven oscillon simulation in
Fig. 8. The $\mathcal{F}^{8/3}$ scaling at small $\mathcal{F}$ is eventually
replaced by approximately constant scaling at larger $\mathcal{F}$ when the
lifetime bound Eq. 47 takes over.

As the oscillon siphons energy from the long axion, it depletes the energy
density in its local environment of volume

$\displaystyle V_{\text{dep}}\sim\frac{P_{\text{rad}}t}{\rho_{L}}.$ (51)

The surrounding long axion flows into the depleted volume at a velocity which
we approximate assuming the inner depleted region is diminished by ${\cal
O}(1)$:

$\displaystyle
v=\frac{\mathrm{d}\omega}{\mathrm{d}k}=\frac{k}{m_{L}}\sim\frac{2\pi}{m\mu
R_{\text{dep}}}.$ (52)

where $V_{\text{dep}}=4\pi/3R_{\text{dep}}^{3}$. Using this velocity, we
calculate the power flowing from the environment into the depleted region:

$\displaystyle P_{\text{env}\to\text{dep}}=4\pi R_{\text{dep}}^{2}v\rho_{L}.$
(53)

Solving for the time at which $P_{\text{env}\to\text{dep}}=P_{\text{rad}}$, we
find the following upper bound on the oscillon lifetime

$\displaystyle m\mu
t_{\text{death}}\lesssim\left(\frac{\Theta_{L}^{0}}{\Theta_{S}^{0}}\right)^{2}\mathcal{F}^{2}mT_{\text{inst}}(\mu).$
(54)

###  B.2 General potentials and cosmological longevity

For the simplest two-axion potentials, the lifetime of driven oscillons is
still far shorter than the age of the Universe. Generalizations of Eq. 3 of
the form

$\displaystyle
V(\theta_{S},\theta_{L})=V_{S}(\theta_{S}+\theta_{L})+V_{L}(\theta_{L}),$ (55)

could plausibly be constructed so that driven oscillons survive into the
present day. The “vacuum lifetime” of the short oscillon
$T_{\mathrm{inst}}(\mu)$ is then the instantaneous lifetime of single-field
oscillons in the $V_{S}(\theta_{S})$ potential.

During hierarchical structure formation, a driven short-axion oscillon finds
itself in long-axion halos of increasing size, which could make it
energetically possible for the oscillon to live indefinitely, although there
are some significant uncertainties. A stationary oscillon in a static long
axion halo would eventually deplete its local environment of long axion,
starving itself of the energy it needs to survive. In a realistic galactic
halo, however, the region inside the oscillon is constantly being replenished
via the virial motion of the long axion. To take advantage of this energy
source, the oscillon must also adjust its phase to match that of the driver
over one virial timescale; it is not clear whether particularly long-lived
oscillons can perform this phase alignment without some efficient dissipation
mechanism (such as dynamical friction due to photon radiation). We thus leave
the question of cosmologically long-lived driven oscillons to future work.

## References

  * Peccei and Quinn [1977a] R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn [1977b] R. D. Peccei and H. R. Quinn, Constraints Imposed by CP Conservation in the Presence of Instantons, Phys. Rev. D 16, 1791 (1977b).
  * Weinberg [1978] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40, 279 (1978).
  * Witten [1984] E. Witten, Some Properties of O(32) Superstrings, Phys. Lett. B 149, 351 (1984).
  * Banks and Dine [1996] T. Banks and M. Dine, Couplings and scales in strongly coupled heterotic string theory, Nucl. Phys. B 479, 173 (1996), arXiv:hep-th/9605136 .
  * Svrcek and Witten [2006] P. Svrcek and E. Witten, Axions In String Theory, JHEP 06, 051, arXiv:hep-th/0605206 .
  * Arvanitaki  _et al._ [2010] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String Axiverse, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie [1983] L. F. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120, 137 (1983).
  * Turner [1983] M. S. Turner, Coherent Scalar Field Oscillations in an Expanding Universe, Phys. Rev. D 28, 1243 (1983).
  * Marsh [2016] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Daido  _et al._ [2015] R. Daido, N. Kitajima, and F. Takahashi, Domain Wall Formation from Level Crossing in the Axiverse, Phys. Rev. D 92, 063512 (2015), arXiv:1505.07670 [hep-ph] .
  * Kitajima  _et al._ [2018] N. Kitajima, J. Soda, and Y. Urakawa, Gravitational wave forest from string axiverse, JCAP 10, 008, arXiv:1807.07037 [astro-ph.CO] .
  * Arvanitaki  _et al._ [2020] A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, Large-misalignment mechanism for the formation of compact axion structures: Signatures from the QCD axion to fuzzy dark matter, Phys. Rev. D 101, 083014 (2020), arXiv:1909.11665 [astro-ph.CO] .
  * Co  _et al._ [2020a] R. T. Co, L. J. Hall, and K. Harigaya, Axion Kinetic Misalignment Mechanism, Phys. Rev. Lett. 124, 251802 (2020a), arXiv:1910.14152 [hep-ph] .
  * Cyncynates  _et al._ [2022] D. Cyncynates, T. Giurgica-Tiron, O. Simon, and J. O. Thompson, Resonant nonlinear pairs in the axiverse and their late-time direct and astrophysical signatures, Phys. Rev. D 105, 055005 (2022), arXiv:2109.09755 [hep-ph] .
  * Co  _et al._ [2020b] R. T. Co, L. J. Hall, K. Harigaya, K. A. Olive, and S. Verner, Axion Kinetic Misalignment and Parametric Resonance from Inflation, JCAP 08, 036, arXiv:2004.00629 [hep-ph] .
  * Eröncel  _et al._ [2022] C. Eröncel, R. Sato, G. Servant, and P. Sørensen, ALP Dark Matter from Kinetic Fragmentation: Opening up the Parameter Window, (2022), arXiv:2206.14259 [hep-ph] .
  * Eröncel and Servant [2022] C. Eröncel and G. Servant, ALP Dark Matter Mini-Clusters from Kinetic Fragmentation, (2022), arXiv:2207.10111 [hep-ph] .
  * Kudryavtsev [1975] A. Kudryavtsev,  _Solitonlike solutions for a Higgs scalar field_ , Tech. Rep. (Institute of Theoretical and Experimental Physics, 1975).
  * Makhankov [1978] V. G. Makhankov, Dynamics of Classical Solitons In Nonintegrable Systems, Phys. Rept. 35, 1 (1978).
  * Gleiser [1994] M. Gleiser, Pseudostable bubbles, Phys. Rev. D 49, 2978 (1994), arXiv:hep-ph/9308279 .
  * Kolb and Tkachev [1994] E. W. Kolb and I. I. Tkachev, Nonlinear axion dynamics and formation of cosmological pseudosolitons, Phys. Rev. D 49, 5040 (1994), arXiv:astro-ph/9311037 .
  * Salmi and Hindmarsh [2012] P. Salmi and M. Hindmarsh, Radiation and Relaxation of Oscillons, Phys. Rev. D 85, 085033 (2012), arXiv:1201.1934 [hep-th] .
  * Amin  _et al._ [2012] M. A. Amin, R. Easther, H. Finkel, R. Flauger, and M. P. Hertzberg, Oscillons After Inflation, Phys. Rev. Lett. 108, 241302 (2012), arXiv:1106.3335 [astro-ph.CO] .
  * Kawasaki  _et al._ [2020] M. Kawasaki, W. Nakano, and E. Sonomoto, Oscillon of Ultra-Light Axion-like Particle, JCAP 01, 047, arXiv:1909.10805 [astro-ph.CO] .
  * Olle  _et al._ [2021] J. Olle, O. Pujolas, and F. Rompineve, Recipes for oscillon longevity, JCAP 09, 015, arXiv:2012.13409 [hep-ph] .
  * Zhang  _et al._ [2020] H.-Y. Zhang, M. A. Amin, E. J. Copeland, P. M. Saffin, and K. D. Lozanov, Classical Decay Rates of Oscillons, JCAP 07, 055, arXiv:2004.01202 [hep-th] .
  * Cyncynates and Giurgica-Tiron [2021] D. Cyncynates and T. Giurgica-Tiron, Structure of the oscillon: The dynamics of attractive self-interaction, Phys. Rev. D 103, 116011 (2021), arXiv:2104.02069 [hep-ph] .
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2020, 083C01 (2020).
  * Traschen and Brandenberger [1990] J. H. Traschen and R. H. Brandenberger, Particle Production During Out-of-equilibrium Phase Transitions, Phys. Rev. D 42, 2491 (1990).
  * Kofman  _et al._ [1994] L. Kofman, A. D. Linde, and A. A. Starobinsky, Reheating after inflation, Phys. Rev. Lett. 73, 3195 (1994), arXiv:hep-th/9405187 .
  * Kofman  _et al._ [1997] L. Kofman, A. D. Linde, and A. A. Starobinsky, Towards the theory of reheating after inflation, Phys. Rev. D 56, 3258 (1997), arXiv:hep-ph/9704452 .
  * Bassett  _et al._ [2006] B. A. Bassett, S. Tsujikawa, and D. Wands, Inflation dynamics and reheating, Rev. Mod. Phys. 78, 537 (2006), arXiv:astro-ph/0507632 .
  * Allahverdi  _et al._ [2010] R. Allahverdi, R. Brandenberger, F.-Y. Cyr-Racine, and A. Mazumdar, Reheating in Inflationary Cosmology: Theory and Applications, Ann. Rev. Nucl. Part. Sci. 60, 27 (2010), arXiv:1001.2600 [hep-th] .
  * Amin  _et al._ [2014] M. A. Amin, M. P. Hertzberg, D. I. Kaiser, and J. Karouby, Nonperturbative Dynamics Of Reheating After Inflation: A Review, Int. J. Mod. Phys. D 24, 1530003 (2014), arXiv:1410.3808 [hep-ph] .
  * Lozanov [2019] K. D. Lozanov, Lectures on Reheating after Inflation, (2019), arXiv:1907.04402 [astro-ph.CO] .
  * Adshead  _et al._ [2020a] P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner, Constraining axion inflation with gravitational waves from preheating, Phys. Rev. D 101, 083534 (2020a), arXiv:1909.12842 [astro-ph.CO] .
  * Adshead  _et al._ [2020b] P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner, Constraining Axion Inflation with Gravitational Waves across 29 Decades in Frequency, Phys. Rev. Lett. 124, 171301 (2020b), arXiv:1909.12843 [astro-ph.CO] .
  * Weiner [2021] Z. J. Weiner, Stencil Solvers for PDEs on GPUs: An Example From Cosmology, Computing in Science & Engineering 23, 55 (2021).
  * Brouwer  _et al._ [2022a] L. Brouwer _et al._ , Introducing DMRadio-GUT, a search for GUT-scale QCD axions, (2022a), arXiv:2203.11246 [hep-ex] .
  * Alesini  _et al._ [2017] D. Alesini, D. Babusci, D. Di Gioacchino, C. Gatti, G. Lamanna, and C. Ligi, The KLASH Proposal, (2017), arXiv:1707.06010 [physics.ins-det] .
  * Stern [2016] I. Stern, ADMX Status, PoS ICHEP2016, 198 (2016), arXiv:1612.08296 [physics.ins-det] .
  * BRA [2021] The BRASS experiment, https://www1.physik.uni-hamburg.de/iexp/gruppe-horns/forschung/brass.html (2021), accessed: 2021-08-11.
  * Lasenby [2020] R. Lasenby, Microwave cavity searches for low-frequency axion dark matter, Phys. Rev. D 102, 015008 (2020), arXiv:1912.11056 [hep-ph] .
  * Berlin  _et al._ [2020] A. Berlin, R. T. D’Agnolo, S. A. R. Ellis, C. Nantista, J. Neilson, P. Schuster, S. Tantawi, N. Toro, and K. Zhou, Axion Dark Matter Detection by Superconducting Resonant Frequency Conversion, JHEP 07 (07), 088, arXiv:1912.11048 [hep-ph] .
  * Berlin  _et al._ [2021] A. Berlin, R. T. D’Agnolo, S. A. R. Ellis, and K. Zhou, Heterodyne broadband detection of axion dark matter, Phys. Rev. D 104, L111701 (2021), arXiv:2007.15656 [hep-ph] .
  * Brouwer  _et al._ [2022b] L. Brouwer _et al._ (DMRadio), DMRadio-m3: A Search for the QCD Axion Below $1\,\mu$eV, (2022b), arXiv:2204.13781 [hep-ex] .
  * Beurthey  _et al._ [2020] S. Beurthey _et al._ , MADMAX Status Report, (2020), arXiv:2003.10894 [physics.ins-det] .
  * McAllister  _et al._ [2017] B. T. McAllister, G. Flower, J. Kruger, E. N. Ivanov, M. Goryachev, J. Bourhill, and M. E. Tobar, The ORGAN Experiment: An axion haloscope above 15 GHz, Phys. Dark Univ. 18, 67 (2017), arXiv:1706.00209 [physics.ins-det] .
  * Nagano  _et al._ [2019] K. Nagano, T. Fujita, Y. Michimura, and I. Obata, Axion Dark Matter Search with Interferometric Gravitational Wave Detectors, Phys. Rev. Lett. 123, 111301 (2019), arXiv:1903.02017 [hep-ph] .
  * Liu  _et al._ [2019] H. Liu, B. D. Elwood, M. Evans, and J. Thaler, Searching for Axion Dark Matter with Birefringent Cavities, Phys. Rev. D 100, 023548 (2019), arXiv:1809.01656 [hep-ph] .
  * Khlebnikov and Tkachev [1997] S. Y. Khlebnikov and I. I. Tkachev, Relic gravitational waves produced after preheating, Phys. Rev. D 56, 653 (1997), arXiv:hep-ph/9701423 .
  * Easther and Lim [2006] R. Easther and E. A. Lim, Stochastic gravitational wave production after inflation, JCAP 04, 010, arXiv:astro-ph/0601617 .
  * Easther  _et al._ [2007] R. Easther, J. T. Giblin, Jr., and E. A. Lim, Gravitational Wave Production At The End Of Inflation, Phys. Rev. Lett. 99, 221301 (2007), arXiv:astro-ph/0612294 .
  * Easther  _et al._ [2008] R. Easther, J. T. Giblin, and E. A. Lim, Gravitational Waves From the End of Inflation: Computational Strategies, Phys. Rev. D 77, 103519 (2008), arXiv:0712.2991 [astro-ph] .
  * Garcia-Bellido and Figueroa [2007] J. Garcia-Bellido and D. G. Figueroa, A stochastic background of gravitational waves from hybrid preheating, Phys. Rev. Lett. 98, 061302 (2007), arXiv:astro-ph/0701014 .
  * Dufaux  _et al._ [2007] J. F. Dufaux, A. Bergman, G. N. Felder, L. Kofman, and J.-P. Uzan, Theory and Numerics of Gravitational Waves from Preheating after Inflation, Phys. Rev. D 76, 123517 (2007), arXiv:0707.0875 [astro-ph] .
  * Zhou  _et al._ [2013] S.-Y. Zhou, E. J. Copeland, R. Easther, H. Finkel, Z.-G. Mou, and P. M. Saffin, Gravitational Waves from Oscillon Preheating, JHEP 10, 026, arXiv:1304.6094 [astro-ph.CO] .
  * Lozanov and Amin [2019] K. D. Lozanov and M. A. Amin, Gravitational perturbations from oscillons and transients after inflation, Phys. Rev. D 99, 123504 (2019), arXiv:1902.06736 [astro-ph.CO] .
  * Amin  _et al._ [2018] M. A. Amin, J. Braden, E. J. Copeland, J. T. Giblin, C. Solorio, Z. J. Weiner, and S.-Y. Zhou, Gravitational waves from asymmetric oscillon dynamics?, Phys. Rev. D 98, 024040 (2018), arXiv:1803.08047 [astro-ph.CO] .
  * Hiramatsu  _et al._ [2021] T. Hiramatsu, E. I. Sfakianakis, and M. Yamaguchi, Gravitational wave spectra from oscillon formation after inflation, JHEP 03, 021, arXiv:2011.12201 [hep-ph] .
  * Giblin and Thrane [2014] J. T. Giblin and E. Thrane, Estimates of maximum energy density of cosmological gravitational-wave backgrounds, Phys. Rev. D 90, 107502 (2014), arXiv:1410.4779 [gr-qc] .
  * Janssen  _et al._ [2015] G. Janssen _et al._ , Gravitational wave astronomy with the SKA, PoS AASKA14, 037 (2015), arXiv:1501.00127 [astro-ph.IM] .
  * Schmitz [2021] K. Schmitz, New Sensitivity Curves for Gravitational-Wave Signals from Cosmological Phase Transitions, JHEP 01, 097, arXiv:2002.04615 [hep-ph] .
  * Maggiore [2000] M. Maggiore, Gravitational wave experiments and early universe cosmology, Phys. Rept. 331, 283 (2000), arXiv:gr-qc/9909001 .
  * Dvorkin  _et al._ [2022] C. Dvorkin _et al._ , The Physics of Light Relics, in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.07943 [hep-ph] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Pagano  _et al._ [2016] L. Pagano, L. Salvati, and A. Melchiorri, New constraints on primordial gravitational waves from Planck 2015, Phys. Lett. B 760, 823 (2016), arXiv:1508.02393 [astro-ph.CO] .
  * Kite  _et al._ [2021] T. Kite, A. Ravenni, S. P. Patil, and J. Chluba, Bridging the gap: spectral distortions meet gravitational waves, Mon. Not. Roy. Astron. Soc. 505, 4396 (2021), arXiv:2010.00040 [astro-ph.CO] .
  * Smith  _et al._ [2006] T. L. Smith, M. Kamionkowski, and A. Cooray, Direct detection of the inflationary gravitational wave background, Phys. Rev. D 73, 023504 (2006), arXiv:astro-ph/0506422 .
  * Caprini and Figueroa [2018] C. Caprini and D. G. Figueroa, Cosmological Backgrounds of Gravitational Waves, Class. Quant. Grav. 35, 163001 (2018), arXiv:1801.04268 [astro-ph.CO] .
  * Clarke  _et al._ [2020] T. J. Clarke, E. J. Copeland, and A. Moss, Constraints on primordial gravitational waves from the Cosmic Microwave Background, JCAP 10, 002, arXiv:2004.11396 [astro-ph.CO] .
  * Iršič  _et al._ [2017] V. Iršič, M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker, First constraints on fuzzy dark matter from Lyman-$\alpha$ forest data and hydrodynamical simulations, Phys. Rev. Lett. 119, 031302 (2017), arXiv:1703.04683 [astro-ph.CO] .
  * Armengaud  _et al._ [2017] E. Armengaud, N. Palanque-Delabrouille, C. Yèche, D. J. E. Marsh, and J. Baur, Constraining the mass of light bosonic dark matter using SDSS Lyman-$\alpha$ forest, Mon. Not. Roy. Astron. Soc. 471, 4606 (2017), arXiv:1703.09126 [astro-ph.CO] .
  * Dentler  _et al._ [2022] M. Dentler, D. J. E. Marsh, R. Hložek, A. Laguë, K. K. Rogers, and D. Grin, Fuzzy dark matter and the Dark Energy Survey Year 1 data, Mon. Not. Roy. Astron. Soc. 515, 5646 (2022), arXiv:2111.01199 [astro-ph.CO] .
  * Bozek  _et al._ [2015] B. Bozek, D. J. E. Marsh, J. Silk, and R. F. G. Wyse, Galaxy UV-luminosity function and reionization constraints on axion dark matter, Mon. Not. Roy. Astron. Soc. 450, 209 (2015), arXiv:1409.3544 [astro-ph.CO] .
  * Hlozek  _et al._ [2015] R. Hlozek, D. Grin, D. J. E. Marsh, and P. G. Ferreira, A search for ultralight axions using precision cosmological data, Phys. Rev. D 91, 103512 (2015), arXiv:1410.2896 [astro-ph.CO] .
  * Kobayashi  _et al._ [2017] T. Kobayashi, R. Murgia, A. De Simone, V. Iršič, and M. Viel, Lyman-$\alpha$ constraints on ultralight scalar dark matter: Implications for the early and late universe, Phys. Rev. D 96, 123514 (2017), arXiv:1708.00015 [astro-ph.CO] .
  * Baxter  _et al._ [2019] E. J. Baxter _et al._ (DES), Dark Energy Survey Year 1 Results: Methodology and Projections for Joint Analysis of Galaxy Clustering, Galaxy Lensing, and CMB Lensing Two-point Functions, Phys. Rev. D 99, 023508 (2019), arXiv:1802.05257 [astro-ph.CO] .
  * Laguë  _et al._ [2022] A. Laguë, J. R. Bond, R. Hložek, K. K. Rogers, D. J. E. Marsh, and D. Grin, Constraining ultralight axions with galaxy surveys, JCAP 01 (01), 049, arXiv:2104.07802 [astro-ph.CO] .
  * Hlozek  _et al._ [2018] R. Hlozek, D. J. E. Marsh, and D. Grin, Using the Full Power of the Cosmic Microwave Background to Probe Axion Dark Matter, Mon. Not. Roy. Astron. Soc. 476, 3063 (2018), arXiv:1708.05681 [astro-ph.CO] .
  * Dalal and Kravtsov [2022] N. Dalal and A. Kravtsov, Excluding fuzzy dark matter with sizes and stellar kinematics of ultrafaint dwarf galaxies, Phys. Rev. D 106, 063517 (2022), arXiv:2203.05750 [astro-ph.CO] .
  * Flitter and Kovetz [2022] J. Flitter and E. D. Kovetz, Closing the window on fuzzy dark matter with the 21-cm signal, Phys. Rev. D 106, 063504 (2022), arXiv:2207.05083 [astro-ph.CO] .
  * Hook  _et al._ [2021] A. Hook, G. Marques-Tavares, and D. Racco, Causal gravitational waves as a probe of free streaming particles and the expansion of the Universe, JHEP 02, 117, arXiv:2010.03568 [hep-ph] .
  * Weiner  _et al._ [2021] Z. J. Weiner, P. Adshead, and J. T. Giblin, Constraining early dark energy with gravitational waves before recombination, Phys. Rev. D 103, L021301 (2021), arXiv:2008.01732 [astro-ph.CO] .
  * Akrami  _et al._ [2020] Y. Akrami _et al._ (Planck), Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Ade  _et al._ [2018] P. A. R. Ade _et al._ (BICEP2, Keck Array), BICEP2 / Keck Array x: Constraints on Primordial Gravitational Waves using Planck, WMAP, and New BICEP2/Keck Observations through the 2015 Season, Phys. Rev. Lett. 121, 221301 (2018), arXiv:1810.05216 [astro-ph.CO] .
  * Abazajian  _et al._ [2022] K. Abazajian _et al._ (CMB-S4), CMB-S4: Forecasting Constraints on Primordial Gravitational Waves, Astrophys. J. 926, 54 (2022), arXiv:2008.12619 [astro-ph.CO] .
  * Friedland and Shagalov [2005] L. Friedland and A. Shagalov, Excitation of multiphase waves of the nonlinear schrödinger equation by capture into resonances, Physical Review E 71, 036206 (2005).
  * Friedland and Shagalov [2006] L. Friedland and A. Shagalov, Emergence and control of breather and plasma oscillations by synchronizing perturbations, Physical Review E 73, 066612 (2006).
  * Maslov  _et al._ [2007] E. M. Maslov, L. A. Kalyakin, and A. G. Shagalov, Breather resonant phase locking by an external perturbation, Theoretical and Mathematical Physics 152, 1173 (2007).
  * Batalov and Shagalov [2011] S. Batalov and A. Shagalov, Control of solitons, Physical Review E 84, 016603 (2011).
  * Karachalios  _et al._ [2019] N. I. Karachalios, P. Kyriazopoulos, and K. Vetas, Excitation of peregrine-type waveforms from vanishing initial conditions in the presence of periodic forcing, Zeitschrift für Naturforschung A 74, 371 (2019).
  * Fodor  _et al._ [2009] G. Fodor, P. Forgacs, Z. Horvath, and M. Mezei, Radiation of scalar oscillons in 2 and 3 dimensions, Phys. Lett. B 674, 319 (2009), arXiv:0903.0953 [hep-th] .
  * Fodor [2019] G. Fodor,  _A review on radiation of oscillons and oscillatons_ , Ph.D. thesis, Wigner RCP, Budapest (2019), arXiv:1911.03340 [hep-th] .
  * Van Tilburg  _et al._ [2018] K. Van Tilburg, A.-M. Taki, and N. Weiner, Halometry from Astrometry, JCAP 07, 041, arXiv:1804.01991 [astro-ph.CO] .
  * Prabhu [2020] A. Prabhu, Optical Lensing by Axion Stars: Observational Prospects with Radio Astrometry, (2020), arXiv:2006.10231 [astro-ph.CO] .
  * Prabhu and Rapidis [2020] A. Prabhu and N. M. Rapidis, Resonant Conversion of Dark Matter Oscillons in Pulsar Magnetospheres, JCAP 10, 054, arXiv:2005.03700 [astro-ph.CO] .
  * Buckley  _et al._ [2021] J. H. Buckley, P. S. B. Dev, F. Ferrer, and F. P. Huang, Fast radio bursts from axion stars moving through pulsar magnetospheres, Phys. Rev. D 103, 043015 (2021), arXiv:2004.06486 [astro-ph.HE] .
  * Amin and Mou [2021] M. A. Amin and Z.-G. Mou, Electromagnetic Bursts from Mergers of Oscillons in Axion-like Fields, JCAP 02, 024, arXiv:2009.11337 [astro-ph.CO] .
  * Amin  _et al._ [2021] M. A. Amin, A. J. Long, Z.-G. Mou, and P. Saffin, Dipole radiation and beyond from axion stars in electromagnetic fields, JHEP 06, 182, arXiv:2103.12082 [hep-ph] .
  * Towns  _et al._ [2014] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peterson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr, Xsede: Accelerating scientific discovery, Computing in Science & Engineering 16, 62 (2014).
  * Klöckner  _et al._ [2012] A. Klöckner, N. Pinto, Y. Lee, B. Catanzaro, P. Ivanov, and A. Fasih, PyCUDA and PyOpenCL: A Scripting-Based Approach to GPU Run-Time Code Generation, Parallel Computing 38, 157 (2012).
  * Klöckner [2014] A. Klöckner, Loo.py: transformation-based code generation for GPUs and CPUs, in  _Proceedings of ARRAY ‘14: ACM SIGPLAN Workshop on Libraries, Languages, and Compilers for Array Programming_ (Association for Computing Machinery, Edinburgh, Scotland., 2014).
  * Dalcín  _et al._ [2008] L. Dalcín, R. Paz, M. Storti, and J. D’Elía, MPI for Python: Performance improvements and MPI-2 extensions, Journal of Parallel and Distributed Computing 68, 655 (2008).
  * Dalcín  _et al._ [2005] L. Dalcín, R. Paz, and M. Storti, MPI for Python, Journal of Parallel and Distributed Computing 65, 1108 (2005).
  * Dalcin, Lisandro and Mortensen, Mikael and Keyes, David E [2019] Dalcin, Lisandro and Mortensen, Mikael and Keyes, David E, Fast parallel multidimensional FFT using advanced MPI, Journal of Parallel and Distributed Computing 10.1016/j.jpdc.2019.02.006 (2019).
  * Harris  _et al._ [2020] C. R. Harris _et al._ , Array programming with NumPy, Nature 585, 357 (2020), arXiv:2006.10256 [cs.MS] .
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python, Nature Meth. 17, 261 (2020), arXiv:1907.10121 [cs.MS] .
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2D Graphics Environment, Comput. Sci. Eng. 9, 90 (2007).
  * Meurer  _et al._ [2017] A. Meurer _et al._ , SymPy: symbolic computing in Python, PeerJ Comput. Sci. 3, e103 (2017).
  * van der Velden [2020] E. van der Velden, CMasher: Scientific colormaps for making accessible, informative and ’cmashing’ plots, The Journal of Open Source Software 5, 2004 (2020), arXiv:2003.01069 [eess.IV] .
  * Abramo  _et al._ [1997] L. R. W. Abramo, R. H. Brandenberger, and V. F. Mukhanov, The Energy - momentum tensor for cosmological perturbations, Phys. Rev. D 56, 3248 (1997), arXiv:gr-qc/9704037 .
  * Brandenberger and Takahashi [2018] R. Brandenberger and T. Takahashi, Back-Reaction of Gravitational Waves Revisited, JCAP 07, 040, arXiv:1805.02424 [astro-ph.CO] .
  * Carpenter and Kennedy [1994] M. Carpenter and C. Kennedy, Fourth-order 2N-storage Runge-Kutta schemes, NASA reports TM, 109112 (1994).