Nonperturbative structure in coupled axion sectors and implications for direct detection

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 H 0 ) oscillate near the bottom of their potential and may be treated as free, davidcyn@uw.edu osimon@stanford.edu jedidiah@stanford.edu zweiner@uw.edu 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: 1 The canonically normalized axion field variables φ S and φ L are naturally recast as angular variables via the definition θ S ≡ φ S /f S and θ L ≡ φ L /f L . Redefining Λ 4 1 ≡ m 2 f 2 and Λ 4 2 ≡ µ 2 m 2 F 2 f 2 , the axion masses are 2 m S = m and m L = µm S and their decay constants are f S = f and f L = Ff , respectively. In terms of these variables, Eq. 1 takes the form V (θ L , θ S ) = m 2 f 2 (1 − cos (θ S + θ L )) + µ 2 F 2 (1 − cos θ L ) .
(2) 1 Throughout, we work in units where = c = 1. We also define the reduced Planck mass M pl = 1/ √ 8πG ≈ 2.44 × 10 18 GeV. 2 The interaction-basis axions φ S and φ L are not exact mass eigenstates, making this definition ambiguous. For F 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.
We focus on the range F > 1 where f S < f L , and we refer to φ S and φ 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 µ < 1, corresponding to an O(1) coincidence in their masses. While Eq. 1 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 Λ 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 θ I (t initial ) = 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 O(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 F 3 and sufficiently close masses 0.75 µ < 1 tend to align their frequencies in a process called autoresonance, illustrated in Fig. 1. 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 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 exam- Homogeneous dynamics (i.e., ignoring the effect of spatial fluctuations) of friendly axions with mass ratio µ = 0.75 and decay constant ratio F = 20. Here t and τ 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 φI /fI . Transparent curves of the same color denote the corresponding results for uncoupled axions in the top panel.
ple, axions are expected to couple to SM photons via an interaction of the form: where g aγγ = C aγγ α QED /2πf a is the axion-photon coupling [13,32] and C aγγ an order-unity constant. As discussed in Ref. [18], when all axions evolve independently from O(1) initial misalignment angles, the final energy density ρ a,0 of each axion is proportional to f 2 a . In this case, the signal strength ρ a,0 g 2 aγγ 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 ρ S,0 ∝ f 2 L when autoresonance completes, but the coupling to SM photon is g Sγγ α/2πf S . Thus the signal strengthρ S,0 g 2 Sγγ of the short axion is enhanced by 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 ∼ 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 O(m 2 f 2 ) and radii of 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 F 6 and that oscillons quench autoresonance if they form early enough, setting a limit on the energy density transfer for F 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 π/m and conformal timestep ∆τ = ∆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 m (see App. A for further details). The 3D evolution begins when H = m, corresponding to a conformal time mτ m = 1 and cosmic time mt m = 1/2. The scale factor is normalized relative to a m ≡ a(t m ).
Of the free parameters in the model, the decay constant ratio 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 µ = 0.75, θ L (0, x) = 0.8 π, and θ S (0, x) = 0, and run simulations for varying values of F. For F 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 F, fluctuations of φ 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. 4 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 parametof the simulations, as smaller µ cause perturbations to grow faster (see Appendix C of Ref. [18]) and shortens the oscillon lifetime (explained in Sec. III C below). 4 To be explicit, we display the energy density projected (averaged) along one axis of the simulation volume, e.g., 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).
ric resonance as the amplification of local overdensities.
Fluctuations become nonlinear at a time mt nl ∼ 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 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 A, comparing to the results of Ref. [18]. In Sec. III B 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 C 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.

A. 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 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 where the heavy state ν h is composed mostly of the short axion, and the light state ν l is composed mostly of the long axion in the limit F 1 (see Appendix A of Ref. [18] for a complete discussion). 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 approximately 5 in good agreement with mt nl ≈ 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 F in detail. For 6 F 20, oscillons nucleate after the short axion's energy density first exceeds the long axion's. At roughly the same time, autoresonance ends andρ 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 (F = 10), where the heavy and light axions' energy densities asymptote toward a common value. Each panel depicts the simulation result for the heavy and light states (ν h and ν 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 ρS(x)/ρ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.
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 ω < m due to its binding energy. Since the long axion's natural frequency is µm < m, an oscillon can provide a locus for more efficient energy transfer from the short axion back to the long axion. 6 Indeed, for most decay constant ratios 6 F 20, the end of autoresonance and formation of oscillons is associated with a substantial transfer of energy to the light axion. For 6 F 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 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 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 F 20, t 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 (F = 40), the energy density in both axions then decays as approximately a −3 . In this case the backreaction effects at play for smaller F are too suppressed to enable substantial energy transfer by the oscillons. The trends for yet larger F are qualitatively similar: parametric resonance proceeds at the same rate and oscillons form at a similar time. The final ratio of energy densitiesρ h /ρ l thus receives a constant boost due to the period of autoresonance but still decreases as 1/F 2 .
Having discussed the dynamics that control the distribution of energy between the two axions, we now summarize the full F-dependence of the late-time energy fractions, where I = h, l. The final partitioning changes qualitatively at a critical decay constant ratio F for which nonlinearities become important (at t nl ) just as the heavy axion's energy density first matches the light one's (via autoresonant energy transfer). From our simulations we find F ≈ 20 for µ = 0.75; this value depends on the mass ratio in the same manner as t nl (c.f. Eq. 6). This timing separates two distinct regimes: one of near-equilibration due to nonlinear effects at F < F and a 1/F 2 suppression of the heavy-axion abundance via the early end of autoresonance at larger F. Both regimes are well captured by 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. For F 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 F 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 F beyond the sophistication of Eq. 8. At larger F 20, the 1/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 F 2 /1.3 ∼ 310 times more energy density than it would have had in the absence of autoresonance. While the above results fix the mass ratio to µ = 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 F of each regime to vary mildly with µ. In particular, F varies in accordance with the change in the timescale of nonlinearity t nl , which itself varies with µ via Eq. 6. Limited simulations corroborate these expectations for 0.70 µ 0.80, the range for which simulation runtimes are tractable, but a more quantitative assessment of the variation of trends with µ 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 µ ≤ 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 F F 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 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 eV m 10 −3 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. 8 to presentday 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.,ρ tot (t m ) ∼ µ 2 m 2 F 2 f 2 Θ 2 L,0 , which subsequently redshifts as a −3 . 7 Combined with Eq. 8, the present-day abundance of the each axion is 7 The a −3 redshifting assumes that the axions are always noninteracting and nonrelativistic. In reality, at early times the axion interactions during autoresonance causeρtot to redshift slightly slower than a −3 , and at later times oscillons radiate mildly relativistic axions such thatρtot redshifts slightly faster than a −3 (until all axions become nonrelativistic). Together, these effects amount to only an O(1) factor which we neglect for simplicity.
with Ξ I set by Eq. 8. 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 O(1) factor over the mass range of interest and are omitted for simplicity.

B. 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. Appendix 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 shortaxion mass m and long axion's decay constant Ff . Gravitational waves are sourced by the anisotropic part of the stress tensor (via Eq. A4), 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. A14) therefore scales with two powers of the fractional energy density of the source at the time of emission, ρ source (t)/ρ(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 nl , when axion gradients are largest and power is scattered to smaller scales.
To proceed, we follow the model-independent heuristics of Refs. [38,65]. Approximating the gravitational wave source as a Gaussian peaking at momentum k with width σ, one may estimate the peak amplitude to be [65]: at the time the source is maximized.
Here γ is what fraction of the Universe's energy is in the source at the time of the process, ν measures how anisotropic the source is, and H p the Hubble parameter at the time of the process. 8 The peak wave number k and width σ are straightforward to approximate (or read off of simulation results), but the anisotropy coefficient ν is harder to estimate; Ref. [65] motivates ν ∼ 10 −2 to 10 −1 for typical processes. Evaluating Eq. 10 at t nl and plugging in aH = 1/ √ 2t nl , where we have used that the energy density in the short axion at t nl is approximately 1/[1 + 1.3(F/F ) 2 ] of the total axion energy density, and the total axion energy density is given byρ tot (t m )(t m /t nl ) 3/2 . Notice that the suppression from how far inside the horizon the peak is (the factor of [aH p /k ] 2 in Eq. 10) is exactly compensated by the growth in time of γ (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 nl . By comparing Eq. 11 with the relic abundance (Eq. 9) we may estimate the peak of the gravitational wave signal as a function of the mass m and relic abundance of the heavy axion Ω h (t 0 ): 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 km/s/Mpc = 2.13 × 10 −33 eV and h = H 0 /H 100 . From our simulations we observe that k /m ≈ 9 and σ ≈ k /3. Taking ν = 1/20 (for which the estimates agree well with the simulation results) and considering the regime F 20 for which the signal is not suppressed, the peak amplitude is at a present-day frequency of the peak of The amplitude estimate Eq. 13 agrees quite well-within a factor of a few-with the spectra from simulations when evaluated at t 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 nl . From Eqs. 13 and 14 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 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 shortlived and redshifting more slowly than the SM plasma from the time of gravitational wave production to the present day. 9 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 Bmode polarization of the CMB already constrain gravitational waves at frequencies between 10 −18 and 10 −16 Hz, with a most stringent upper limit of Ω GW,0 h 2 ∼ 10 −16 at f GW ∼ 10 −17 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 H rec ∼ 3 × 10 −29 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 eV. Consulting Eq. 14, the peak of the signal at such a mass corresponds to a frequency of order 2 × 10 −16 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 m) follows a simple power law, independent of the underlying dynamics [87]. We therefore extrapolate the signals computed in simulations as decaying with f 3 GW 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 Each curve depicts the signal for a short axion mass making up a fraction of the total dark matter FDM 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 3 GW 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. the simulation do nearly follow an f 3 GW 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 GW . 10 Since the CMB becomes increasingly less sensitive to Ω 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 F = 20. To illustrate the con-10 See Ref. [87] for a thorough presentation of the imprint of the expansion history and presence of free-streaming radiation on causal gravitational waves.
straining power on the present abundance of friendly axions, we consider masses varying from 10 −27 to 10 −25 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 eV, Planck and BICEP2/Keck [89,90] allow only a 0.1% friendly subcomponent. 11 We may obtain a heuristic bound as a function of mass by extrapolating the f 3 GW tail to the frequency corresponding to the horizon size at matter-radiation equality, f eq = 1.54 × 10 −15 Hz, where Ref. [75] reports an upper limit of Ω bound GW,0 h 2 = 3.2 × 10 −16 . Applying the scalings of Eqs. 13 and 14, the CMB probes which is only relevant for masses m 10 −27 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 10 −25 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 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 eV the peak of the signal is only ∼ 10 −15 , consistent with the preceding discussion.

C. 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. 12 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 ω is smaller than the rest mass of the axion m. After its initial formation, the binding energy per particle inside an isolated oscillon, is a decreasing function of time, and consequently the characteristic size of the oscillon is an increasing function of time. Over time its rate of radiation falls off approximately exponentially with this growing separation of scales [31]: where λ rad = 3ω is the radiation wavelength due to 3 → 1 processes, which dominate the decay for most ω in potentials with parity symmetry. Thus as oscillons age and ω 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 ω crit < m (see e.g., Refs. [31,97,98]). An oscillon therefore typically spends the majority of its lifetime at a frequency close to ω crit . The mechanism of radiation (which is discussed in detail in Refs. [30,31,98]) does not play an important role in our analysis, so we may simply characterize the decay rate of the oscillon by its instantaneous lifetime: where E osc is the total bound energy in the oscillon and P rad is the power radiated by the oscillon, all measured while the oscillon is at a frequency ω. [The precise values of E osc (ω) and P rad (ω) 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, φ S (t, r) = Φ S (t)R S (r), in the background of a homogeneous long axion, φ L (t, r) = Φ L (t). 13 To be self-consistent, we work in the limit where 13 Note that φ S (t, r) is in general not separable, since as the oscillon the short axion does not backreact onto the long axion, i.e., f S f L . This approximation neglects the effects of radiation, which we include via an artificial damping term with coefficient Γ inst . In the small Φ L /f L limit, keeping only linear terms, Here V 0 and V 1 are effective potentials derived by integrating out R S . 14 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 Φ S and Φ 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. The short-axion oscillon oscillates at the same frequency ω = µm as the driver φ L , conclusively 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. 14 To integrate out the radial profile, one substitutes φ S = Φ S (t)R S (r) into the action and integrates over space. The qualitative features of the resulting effective action for Φ 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. demonstrating local nonlinear resonance inside the oscillon. The phase offset between the short and long axion evolves from roughly π/3 at early times (left side of Fig. 6) to about π/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 π/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. 20 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 death in terms of its instantaneous vacuum lifetime T inst (µ). 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 We compare this analytic scaling to simulations in Fig. 7, verifying the [µT inst (µ)] 4/3 dependence in the range 0.73 µ ≤ 0.81. but spherically symmetric simulations verify the scaling of Eq. 21 out to µ = 0.89. 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 µ 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 This scaling is demonstrated in Fig. 8: the spherically symmetric solutions to the full coupled nonlinear wave equations (see Eq. B7) exhibit F 8/3 behavior for F 40, at which point the lifetime saturates the driving limit, Eq. 21. Performing a quantitative analysis of the parametric dependence on F in 3 + 1D simulations would require simulations with larger mass ratio µ, which, as discussed previously, exceed the limits of our computing resources. However, we do observe qualitative behavior-lifetimes that increase with 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 One can check that there is no region of F − T inst parameter space where Eq. 23 dominates the lifetime. Taken together, the lifetime is 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 dynamicsinteresting 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 F 20). For smaller decay constant ratios F 6, spatial fluctuations have a negligible effect on the dynamics, and we recover the results of homogeneous computations from Ref. [18]. At larger ratios F 20, oscillon nucleation prevents the heavy axion from ever attaining a substantial abundance.
The novel dynamics in the intermediate-F regime position friendly axions to be positively identified as twocomponent 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 ≈ 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 nearfuture 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 F 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. 2), 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 F 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.
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  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 Here Φ is the Newtonian potential and h ij the transverse (∂ i h ij = 0) and traceless (δ ij h ij = 0) tensor perturbation. We neglect scalar anisotropic stress and vector perturbations. We define the comoving Hubble parameter H ≡ ∂ τ a/a (while the standard Hubble parameter is 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(τ )/a(τ m ) = τ /τ m . We take τ 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 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.

Equations of motion
The Euler-Lagrange equation for the axions reads where the potential V is defined in Eq. 1 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 2 s ≡ δP/δρ, the Einstein equations for Φ may be rearranged into where both w and c 2 s are 1/3 in the radiation-dominated era. Finally, the tensor perturbations evolve according to Gravitational waves are sourced by the transverse and traceless anisotropic stress tensor Π ij whose Fourier modes are given in terms of the full stress tensor T ij by 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 pl ) 2 , for which reason it is entirely sufficient to evaluate the stress tensor as if in a homogeneous FLRW Universe: The final term purely contributes to the trace of the stressenergy 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. A4 from zero initial conditions.

Initial conditions
Imposing that Φ was frozen outside of the horizon (kτ 1) to its primordial value generated during inflation sets Φ(τ 1/k, k) = Φ 0 (k) and ∂ τ Φ(τ 1/k, k) = 0. In this case, we find the solution where y ≡ √ wkτ . The primordial curvature perturbation is characterized by a dimensionless power spectrum ∆ 2 Φ (k) defined by which, for standard slow-roll inflation, is nearly scale invariant with amplitude ∼ 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 τ = τ 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 aH) to determine initial conditions. The axions are initialized as Gaussian-random fields with mean field value and velocity at τ 1 set according to the solution to the homogeneous system. The fluctuations are set in Fourier space to match the linearized solution φ I,k (τ 1 ), multiplied by a Gaussian-random complex number (normalized such that the mean squared modulus is unity). For each axion and wave vector k on the grid, we therefore set where U 1 (k) and U 2 (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. A8. 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.

Gravitational wave backgrounds
Gravitational waves carry an effective energy density which, deep inside the horizon (i.e., k H), is [75,116,117] where the bar denotes a time average (i.e., over oscillations). The relic abundance of gravitational waves today per logarithmic wave number is Substituting the inverse Fourier transform of h ij into Eq. A12 permits rewriting Eq. A13 in terms of the dimensionless power spectrum of ∂ τ h ij (defined in analogy to Eq. A9) as after plugging inρ(τ ) = 3M 2 pl H(τ ) 2 . The spectrum evaluated in the early Universe is related to that at the present day (at τ 0 ) by the transfer function and would be observed at present-day frequencies related to the wave number k by Here g and g S are the numbers of relativistic degrees of freedom in energy and entropy density, respectively. Note that the present-day abundance of radiation is Ω rad (τ 0 )h 2 ≈ 4.2 × 10 −5 [70] and that H 0 /h ≡ 100 km s −1 /Mpc ≈ 3.24 × 10 −18 Hz.

Numerical implementation
We discretize the evolution equations, Eqs. A2 to A4, onto a three dimensional, regularly spaced grid with periodic boundary conditions. Following Refs. [40,41,88], we evolve the gravitational wave equation of motion (Eq. A4) under the replacement of the transverse-traceless stress tensor (which is only easily calculated in Fourier space) with the full stress tensor (Eq. A7). The transversetraceless 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. A8 for the Newtonian potential (which requires forward and inverse Fourier transforms at each step), we simply evolve Eq. A3 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 C, 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. 20) 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. 20 fails to capture some important effects that ultimately limit the lifetime of the driven oscillon, namely Eqs. 21 to 23. In what follows, we derive the lifetime bounds in Eqs. 21 to 23 in App. B 1 and discuss the possibility of cosmologically long-lived oscillons in more general multiaxion potentials in App. B 2.

Lifetime bounds
Most long-lived oscillons are well approximated by a single-harmonic ansatz, so Φ S oscillates predominantly at a single frequency. On autoresonance, this frequency approximately matches the long-axion mass (µm) up to transient oscillations around this stable point. Thus on autoresonance we may approximate where δ is a phase-offset. While δ 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 δ as a constant. The power transferred from the long to the short axion is where for our purposes we may approximate the interaction potential by a mass-mixing term, Thus, the time-averaged power transfer is where we've taken m 2 f 2 S (Θ 0 S ) 2 V osc ≈ E osc . As shown in Ref. [18], the phase δ tends toward −π/2 toward the end of autoresonance, and thus we take δ = −π/2 when calculating the driven oscillon lifetime. The oscillon is supported by the driver when its radiated power P rad is smaller than the maximum power transfer from the long axion P L→S , from which we obtain the approximate time of death t death , taking t 0 ≈ 1/mµ. Here T inst (ω) = E osc (ω)/P rad (ω) is the instantaneous oscillon lifetime at the frequency ω, which may be longer or shorter than the in-vacuum oscillon lifetime. Provided that the driving frequency mµ 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 θ S onto θ L in order to approximate θ L as a homogeneous field. While the details of backreaction are complicated, it is simple to obtain an estimate for when θ S causes O(1) fluctuations in θ L which then terminate the autoresonance. To make this estimate, we observe in the equations of motion, θ L + m 2 F −2 sin(θ S + θ L ) + m 2 µ 2 sin θ L = 0 (B7a) θ S + m 2 µ 2 sin(θ S + θ L ) = 0, (B7b) that θ L 's evolution depends on θ S multiplied by F −2 . Comparing this term to the final term in Eq. B7a, we see that backreaction becomes important when Assuming that θ S has a constant O(1) amplitude, as inside an oscillon, and that θ L decays like cold matter as in Eq. B2, we find that spatial fluctuations in θ L become order one at the time 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 F 8/3 scaling versus the observed lifetime in a spherically symmetric driven oscillon simulation in Fig. 8. The F 8/3 scaling at small F is eventually replaced by approximately constant scaling at larger F when the lifetime bound Eq. B6 takes over.
As the oscillon siphons energy from the long axion, it depletes the energy density in its local environment of volume The surrounding long axion flows into the depleted volume at a velocity which we approximate assuming the inner depleted region is diminished by O (1): where V dep = 4π/3R 3 dep . Using this velocity, we calculate the power flowing from the environment into the depleted region: Solving for the time at which P env→dep = P rad , we find the following upper bound on the oscillon lifetime

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. 2 of the form could plausibly be constructed so that driven oscillons survive into the present day. The "vacuum lifetime" of the short oscillon T inst (µ) is then the instantaneous lifetime of single-field oscillons in the V S (θ S ) potential. During hierarchical structure formation, a driven shortaxion 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.