The Structure of the Oscillon: The Dynamics of Attractive Self-Interaction

Real scalar fields with attractive self-interaction may form self-bound states, called oscillons. These dense objects are ubiquitous in leading theories of dark matter and inflation; of particular interest are long-lived oscillons which survive past $14$ Gyr, offering dramatic astrophysical signatures into the present day. We introduce a new formalism for computing the properties of oscillons with improved accuracy, which we apply to study the internal structure of oscillons and to identify the physical mechanisms responsible for oscillon longevity. In particular, we show how imposing realistic boundary conditions naturally selects a near-minimally radiating solution, and how oscillon longevity arises from its geometry. Further, we introduce a natural vocabulary for the issue of oscillon stability, which we use to predict new features in oscillon evolution. This framework allows for new efficient algorithms, which we use to address questions of whether and to what extent long-lived oscillons are fine-tuned. Finally, we construct a family of potentials supporting ultra-long-lived oscillons, with lifetimes in excess of $10^{17}$ years.

Axions are real scalar fields predicted to exist in many extensions of the Standard Model. One of the bestmotivated is the QCD axion, which emerges as the pseudo-Nambu-Goldstone boson of a broken U (1)-axial symmetry, known as Peccei-Quinn (PQ) symmetry [1]. The PQ breaking scale f a , known as the axion decay constant, suppresses the axion's self-interactions and its coupling to the Standard Model (SM). To avoid impacting stellar cooling rates, axion-SM interactions must be highly suppressed, forcing f a to be in the deep UV, f a 10 10 GeV [2][3][4]. As the universe cools below the QCD scale Λ QCD ≈ 200 MeV, strong dynamics generate a periodic potential for the axion, whose VEV cancels the strong sector's CP-violating phase, thus resolving the strong-CP problem. The separation between the QCD scale and the PQ scale forces the axion's mass m a ∼ Λ 2 QCD /f a to be smaller than 1 meV, potentially by many orders of magnitude [5,6].
Furthermore, axionic degrees of freedom emerge in great numbers from realistic string compactifications, collectively known as the Axiverse. Like the QCD axion, these axion-like particles (ALPs) are generally described by two parameters: their mass m and the decay constant f . Generic ALPs are also expected to have naturally small masses, which are exponentially suppressed by the string instanton action. The precise form of the ALP potential depends on the specifics of the UV theory it descends from, leaving its low-energy dynamics effectively unconstrained [7].
Axions (both the QCD axion and ALPs) come equipped with natural production mechanisms, such as the vacuum misalignment mechanism, making them wellmotivated dark matter candidates [8][9][10][11][12][13][14][15][16]. Of particular phenomenological interest are ultralight axions, whose To interpret this figure, we recognize that each color corresponds to approximately a single lifetime. Therefore, thin regions contain the most significant features, while broad regions, such as the value of the potential near φ/f = ±π, are the least significant for determining the lifetime. As the central part of the potential approaches a free theory, the oscillon must grow in spatial extent because of weak selfinteraction, leading to decoupling of the large bound oscillon from the short wavelength radiation (see section III A). On the other hand, some self-interaction is necessary to delay energetic death, which is why the purple potentials are much longer-lived than the red ones (see sections III C and IV).
As the densest object in this family of bound axionic structures, oscillons promise dramatic astrophysical signatures, and have therefore been the subject of intense scrutiny [64][65][66][67][68][69][70]. Oscillons have a finite lifespan, and such phenomena crucially rely on oscillons that are cosmologically long-lived. Since dark matter axions are constrained by Lyman-α forest measurements to be at least 10 −21 eV in mass, their oscillation period is at most 0.1 years [17][18][19][20][22][23][24]. Therefore, oscillons that survive 14 Gyr until the present day must be stable for at least 10 11 oscillations. Simulating an oscillon this long-lived is at the upper limit of current computational capabilities [57,58], and thus indirect methods are required to study longer-lived oscillons.
Significant progress has been made towards understanding the structure and evolution of oscillons in the last two decades, building on improved computational resources and theoretical understanding [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63]. Of central theoretical importance are artificial, exactly-periodic so-lutions of the equations of motion, which have been used to approximate the oscillon's instantaneous profile and radiation rate. In rare instances, in which the oscillon is known to be infinitely long-lived, this approximation is exact, and the solution is called a breather : a finite energy periodic solution of the equations of motion. The most famous such example is the 1+1 dimensional sine-Gordon breather, which is stabilized by an infinite set of conserved quantities [71]. Breathers are not known to exist in 3+1 dimensions. Relaxing the breather's finite energy constraint, we find the periodic solutions known as quasibreathers. These constructions have an infinite amount of energy residing in their standing-wave tails. These radiative tails can be understood as an approximation of the oscillon's classical radiation amplitude, which can be used to estimate the oscillon's lifetime.
In this paper, we further develop the quasibreather technique into a framework for understanding the classical properties of oscillons, unifying several observations made in the literature, and addressing key conceptual questions about the harmonic structure and stability of oscillons. By imposing realistic boundary conditions, we introduce the physical quasibreather (PQB) as the member of the quasibreather family closest to a radiating oscillon, and arrive at an improved method for calculating oscillon properties, such as lifetime, radial profile, linear stability, and frequency content. In the limit of long lifetimes, our method becomes especially efficient, since semi-perturbative techniques may be employed to rapidly compute oscillon radiation. We apply our new methods to systematically study oscillon lifetimes in periodic axion potentials, allowing us to probe the genericity of long-lived oscillons. Moreover, we apply our framework to expand on existing studies of long-lived oscillons in monodromy potentials.
We summarize our study of oscillon lifetimes in periodic potentials with parity in the form of longevity landscapes, such as the one depicted in figure 1a. There, we scan the coefficients V n of the axion potential V (φ) defined as Here, the field φ is the axion field, m its mass, and f its decay constant. The particular slice through the space of coefficients in figure 1a is defined by the choice to treat V 1 and V 2 as free parameters, while fixing V 3 = −1, and forcing V 4 to satisfy the mass constraint, with all other V n set to zero. Our numerical techniques based on the PQB formalism have allowed us to perform this parameter sweep in 96 CPU-hours, parallelized down to a few hours of wall-clock time. We see that the landscape is broken down into "islands of longevity," where neighboring potentials sustain oscillons that are similarly longlived. While most of this space supports oscillons in the range 10 2 − 10 4 oscillations, these few tunable parame-ters in the potential are enough to allow for oscillons that may live up to 10 14 cycles.
The distinct islands in figure 1a correspond to the action of two mechanisms that suppress oscillon radiation, which we identify as totally destructive self-interference and geometric decoupling. Together, these two effects comprise the form-factor of the oscillon coupling to radiation, but we separate them because of their distinct imprints on the oscillon life-cycle, as depicted in figure  2. Further, the cliffs in figure 1a represent destructive interference peaks entering unreachable frequencies beyond the point of energetic death, where the oscillon is forced to dissipate because of energy conservation. Here, we briefly review these three effects.
Destructive interference: The bound bulk of the oscillon is a nearly coherent object, oscillating at frequency ω. Through the interaction terms of integer order φ n+1 , the oscillon bulk behaves as a nearly coherent source of radiation at multiples of the fundamental frequency nω. Similar to a diffraction experiment, certain geometries lead to totally destructive interference, exponentially confining certain radiation channels at exceptional frequencies. When the dominant radiation channel destructively interferes, the radiated power experiences a sudden 'dip.' Geometric decoupling: The size of the oscillon is inversely proportional to the binding energy per particle √ m 2 − ω 2 , which blows up as ω approaches the rest mass m (see figure 3). In this limit, the oscillon grows much larger than the wavelengths of radiation 2π/nω, causing a separation of scales. As this separation grows, the smooth oscillon bulk decouples from radiation, which manifests as an exponential decrease in radiated power towards the end of the oscillon's lifetime.
Energetic death: As the oscillon radiates away its energy, the binding energy per particle decreases, reducing the oscillon's central amplitude and increasing its radius. In three or more spatial dimensions, weak self-interactions result in a volume growing faster than can be accommodated by the decreasing central amplitude. Therefore, at frequencies ω approaching the mass m, there is a point past which an external energy source is necessary for the oscillon to remain bound. At this point, the oscillon is forced to undergo a rapid process of dissipation, which we call energetic death.
These mechanisms explain the structure of the longevity landscape observed in figure 1a. An island of longevity starts when a point of destructive interference (a 'dip') emerges from low frequencies (green contour). As the dip migrates toward higher frequencies, its effect is enhanced by geometric decoupling, causing lifetime to increase until the dip moves beyond the frequency of energetic death, resulting in a longevity 'cliff' (blue contour).
In order to obtain these results, we have applied our PQB formalism to estimate the evolution of extremely long-lived spherically symmetric oscillons in isolation. In doing so, we have made the implicit assumption that the physical oscillon has relaxed into a state near the PQB. In order to check that this assumption is valid, we have performed a detailed linear stability analysis of the PQB to spherical and non-spherical perturbations, and we have presented evidence that unstable modes remain small enough that our procedure stays predictive (appendix C).
This paper is structured as follows. Section II introduces the main object of study, the physical quasibreather. The oscillon is identified as living in the basin of attraction of the PQB, which naturally captures notions of oscillon stability. Section III uses the PQB formalism to understand the mechanisms of longevity briefly discussed above, and derives the minimum radiation condition. Section IV applies the mechanisms of oscillon longevity and death to construct a family of potentials supporting ultra-long-lived oscillons. Section V applies our techniques to study the genericity of long-lived oscillons, and introduces local and global measures of finetuning. Section VI applies our formalism to well-known potentials in the literature, re-deriving and expanding on previous results. Finally, the appendices provide a detailed technical overview of our formalism, and contain an exhaustive treatment of linear stability, as well as our numerical workflow. Appendix A provides the mathematical basis of the PQB. Appendix B details the numerical procedure for obtaining the PQB. Appendix C details our linear and nonlinear stability analysis of the PQB. Appendix D provides technical formulae relevant for computing the PQB and its linear stability. Appendix E details our explicit numerical simulations.

II. THE PHYSICAL QUASIBREATHER
The nonlinear wave equation we study in this paper is of the generic form Here, f is the scale of self-interaction, known as the axion decay constant. The overdot represents time differentiation, ∇ 2 is the usual flat-space Laplacian, and V represents differentiation of the potential V (φ/f ) with respect to the field φ. An oscillon is a finite-energy solution of (2) that is quasibound by self-interactions. In 3+1 dimensions, which is the focus of our study, all known oscillons have a finite lifetime because they radiate classical scalar waves. To understand whether a potential hosts cosmologically relevant oscillons, one needs a robust computational formalism for obtaining these classical radiation rates. Here, we introduce the physical quasibreather formalism for computing the oscillon radiation and lifetime, while leaving the more technical details to appendix A.
FIG. 2. This plot illustrates the mechanisms of oscillon longevity and death described in section III. Here, we plot the power carried out of the oscillon in the dominant radiating harmonics as a function of the oscillon frequency ω. The fundamental frequency ω increases with time, and therefore may be interpreted as a time coordinate (see figure 3). For simplicity, we consider a scalar potential with parity symmetry, leading to radiation at odd multiples of ω due to n → 1 processes. Towards higher frequencies, the size of the oscillon 2π/ √ m 2 − ω 2 is much larger than the radiation wavelength 2π/(nω), leading to the geometric decoupling of radiation. As the oscillon becomes more diffuse, its volume grows faster than its amplitude shrinks, forcing an early energetic death. At exceptional frequencies, certain radiative harmonics vanish as a consequence of destructive self-interference.
A physical potential V represents interactions between an integer number of particles, and therefore possesses a well defined Taylor series. Consequently, a field oscillating at fundamental frequency ω will only couple to integer multiples of ω. Thus, one may look for quasibreather solutions: spherically symmetric, exactly periodic solutions of the equation of motion (2) of the form where δ n are constant phases, with δ 1 = 0 by the choice of a time coordinate. The harmonic profiles S n (r, ω) divide into bound modes n < m/ω and radiative modes n > m/ω. Solutions of this form were first introduced in [72] and have since been used throughout the oscillon literature to obtain approximate oscillon solutions (see [54] for a complete review). Although (3) is a periodic solution of the equations of motion, it is not an infinitely long-lived oscillon; the far-field tails of the radiative harmonics S n>m/ω decay like r −1 , and therefore contribute infinite energy. These unphysical, infinite energy radiative tails have been problematic when interpreting quasibreathers as approximate oscillons. Furthermore, finding a quasibreather of a specific frequency is underdetermined: there are as many different quasibreathers of frequency FIG. 4. The radial profile of the physical quasibreather (PQB) (solid) and its orthogonal deformation (OD) (dashed) for the sine-Gordon (SG) oscillon at ω = 0.92m, plotted against radius in units of the mass m −1 . In the limit where the radiation tails are small, this serves as an instantaneous approximation of the internal structure of the oscillon. The first quasibreather harmonic S1 is exponentially bound, defining the oscillon bulk. The third harmonic S3 is the dominant radiation mode, followed by the fifth, seventh, and so on. The spatial and temporal phase of the OD are 90 degrees out of phase with the PQB in the radiative region, representing outgoing radiation.
ω as there are radiative degrees of freedom, representing the choice of central amplitudes S n>m/ω (r = 0, ω). One proposal to resolve this ambiguity is to pick the quasibreather with the minimum radiation amplitude, in an attempt to minimize the influence of the unphysical radiation tails (see e.g. [54]). Here, we introduce a different criterion for choosing the quasibreather closest to a physical oscillon. Instead of demanding that the radiative tails FIG. 5. The PQB trajectory of the harmonic amplitudes S1 and S3 (red) is plotted on top of the level sets of the effective potential. The set of all initial conditions corresponding to quasibreathers is outlined in dotted blue. The particular example plotted here is of the sine-Gordon equation for ω = 0.5m. are minimized, we will require that the quasibreather is perturbatively close to a radiating solution of (2).
To this end, we introduce the orthogonal deformation (OD) whose temporal phase is 90 degrees offset from that of the quasibreather (3). Note that the sum over n only includes the frequencies corresponding to modes with radiative tails, nω > m. When added to the standing wave quasibreather (3), the orthogonal deformation allows for travelling modes (see figure 4). We then define the family of physical quasibreathers (PQB), θ PQB , parametrized by ω, as those quasibreathers which may be orthogonally deformed θ PQB → θ PQB + θ OD to satisfy purely outgoing boundary conditions at leading order in θ OD (i.e. θ PQB + θ OD must satisfy the Sommerfeld radiation condition [73]). Note, we will use subscripts to refer either to a general quasibreather θ QB or to a physical quasibreather θ PQB with an OD partner that together satisfy the Sommerfeld radiation condition. The radiative boundary conditions are enforced at spatial infinity, where the wave equation (2) is well approximated by the Klein-Gordon equation. In this region, the OD and the radiative tails of the PQB are of the same amplitude because they represent purely outgoing radiation. Because θ PQB is a solution of the equations of motion, the perturbation θ OD must backreact at second order O(θ 2 OD ), and it must obey a homogeneous linear equation to the same order [74]. Therefore, a PQB with small radiative tails must have an OD that is small everywhere, compared to the PQB central amplitude. The infinite lifetime limit is the limit of no radiation, and in this case, the PQB approaches a finite energy oscillon. Therefore, the PQB will be the central object in our study of long-lived oscillons.
To summarize, the following three objects are pointwise close to one another: the finite energy oscillon, the PQB, and the orthogonally deformed PQB. This proximity forms the basis of an expansion of the oscillon, which we fully develop in appendix A, where the oscillon is understood to be a stable perturbation of the orthogonally deformed PQB. As such, oscillon properties (instantaneous frequency, stability, radiation, etc.) may be understood as originating from the nearest PQB. Furthermore, we develop the Floquet analysis of linear perturbations to the deformed PQB in appendix C. Although linear stability turns out to be a sufficient criterion for the existence of an oscillon, it is not a necessary condition, since stable orbits can (and do) emerge at higher orders. In other words, linear instability does not imply the dissolution of the oscillon, since nonlinearities control the size of the linearly unstable perturbations. This effect has important phenomenological consequences for the nature of the oscillon evolution (for examples, see figures 11,13,18a,18b). Specifically, slow quasiperiodic oscillations around the PQB profile emerge in linearly unstable regions, with amplitude that depends strongly on initial conditions. Below, section II A provides a minimal technical review of our framework, which will be useful in understanding the qualitative features of oscillon evolution in section III. Afterwards, section II B outlines the steps in the numerical workflow of computing the PQB and OD, as well as the associated oscillon properties such as lifetime.

A. The mode equations
At each stage of its life-cycle, the oscillon may be viewed as close to a particular physical quasibreather. This description becomes increasingly precise in the infinite lifetime limit, where radiation goes to zero and the oscillon evolves slowly. Because the oscillon spends a long time in the vicinity of a particular physical quasibreather, the notion of the instantaneous frequency ω becomes well defined. Physically, ω then behaves like an adiabatic parameter, although formally it serves as an index to label which physical quasibreather the oscillon is closest to at a given time. The fact that the oscillon does remain close to the physical quasibreather family is a consequence of its attractive properties, which we make precise in appendix C.
We are now in position to introduce the mode equations, which describe the spatial profile of the physical quasibreather at a given frequency ω. In the interest of a pedagogical introduction, we will consider the particularly simple case of a single bound harmonic S 1 for a potential with parity V (θ) = V (−θ), and we will keep only the first radiative harmonic S 3 .
As outlined above, the potential V is Taylor expandable, and therefore factorizes into a sequence of integer harmonics of the fundamental frequency ω. By restricting to V (θ) = V (−θ), only the odd harmonics are coupled to one another, allowing for the following decomposition where the dots refer to terms proportional to higher frequencies nω, and terms that contain the small harmonics S n , n ≥ 5. Inserting the quasibreather and the orthogonal deformation into the equations of motion, we arrive at the orthogonally deformed mode equations To fully specify the solution to this system, we must provide 6 boundary conditions: regularity at the origin regularity at spatial infinity and radiative boundary conditions [73] To understand these equations, it is helpful to visualize the evolution of S 1 and S 3 as the coordinates of a point particle rolling down a hill, where r is now the time coordinate, and the initial stationary particle is placed so that it arrives at the saddle located at the origin when r → ∞ (see figure 5). Out of the continuum of quasibreather initial conditions S 1 (0), S 3 (0) satisfying this constraint, the orthogonal deformation selects only one, corresponding to the PQB.

B. Calculation workflow
Here we review the workflow of estimating the oscillon lifetime in the physical quasibreather framework, leaving a more detailed presentation to the appendices.
1. The harmonics S n of the PQB may be thought of as existing in two categories. The perturbative harmonics are those S n whose amplitude is everywhere small enough that self-interaction can be safely neglected. Those S n for which this is not true are called non-perturbative. Typically, only a few non-perturbative harmonics are needed to achieve numerical convergence. The physical intuition for whether a harmonic may be treated perturbatively or not is whether it contributes significantly to the binding energy compared to the flux radiated per cycle. In other words, a good rule of thumb for whether a harmonic is perturbative is whether its central amplitude is significantly larger than the leading orthogonal deformation at the origin.
2. The non-perturbative harmonics (which must include S 1 ) are calculated using a shooting technique, in which the S n 's are propagated from the origin to an outer boundary at r = r out . At this point, the Sommerfeld radiation condition (10) is used to calculate the OD, c n (r out ) and c n (r out ). From these final conditions, the c n are propagated back to the origin in the background of the non-perturbative S n . One then checks whether the backwards propagated c n 's satisfy regularity at the origin. We perform a search over initial conditions S n (0) until regularity is satisfied for all c n 's.

3.
Having computed the non-perturbative harmonics, an arbitrary number of perturbative harmonics may be computed to linear order by solving a sparse matrix equation. In other words, once the hard work of computing the non-perturbative harmonics is done, one may compute the full spectrum of the oscillon to arbitrary harmonic order with little computational cost. One may then re-shoot the non-perturbative harmonics in the background of the perturbative harmonics to account for linear back-reaction, repeating until converged.
4. The result of these calculations is a semi-nonperturbative expression for the physical quasibreather S n and its orthogonal deformation c n . The radiation power in each harmonic is easily computed as P n = 2πr 2 (nω) (nω) 2 − 1 S 2 n + c 2 n evaluated at the outer boundary. The sum n P n is the PQB approximation to the total power P radiated by the oscillon.

5.
Having calculated the outgoing power P as a function of the PQB frequency ω, we may approximate the lifetime of the oscillon near the physical quasibreather trajectory as T = dω(dE B /dω)/P , where E B is the bound energy in the oscillon, defined as the difference between the PQB and OD energy (see appendix A 4).
We provide a public implementation of this protocol for the case of a single non-perturbative harmonic in potentials with parity -a fast and easy-to-use tool to obtain ballpark estimates of oscillon properties at larger frequencies [75].

III. THE OSCILLON LIFE-CYCLE
Here we review and expand upon previous literature results [36,37,[45][46][47][48][49][50][51][52][53][54][55][56][57][58] in order to identify the main mechanisms responsible for oscillon longevity and death. We point out two distinct effects contributing to oscillon longevity: geometric decoupling and destructive interference, both of which may be thought of together as the form-factor of the oscillon coupling to radiation. It is important to separate form-factor into these two effects because they intervene at different times, and have different consequences for oscillon evolution. Often, an oscillon's lifetime is dominated by one mechanism or the other, while the longest lived oscillons take advantage of both simultaneously. Separately, as the oscillon ages and grows more diffuse, it will inevitably undergo an energetic death, beyond which its energy would be forced to unphysically increase. These three effects are all pointed out in figure 2, which depicts the typical radiation history of an oscillon. Below, we provide a semi-quantitative overview of these three effects.

A. Geometric decoupling
Recall that the oscillon is a smooth, nearly coherent object, coupling to integer multiples n of its fundamental frequency ω through many-to-one interactions at leading order φ n+1 . As the oscillon radiates binding energy throughout its life, its fundamental frequency increases towards m (see figure 3), and its typical size 2π/ √ m 2 − ω 2 blows up, where √ m 2 − ω 2 is the binding energy per particle. Therefore, a natural separation of scales occurs between the length scale of radiation 2π/nω and the size of the oscillon, leading to an exponential suppression of the oscillon's coupling to radiative modes nω, n ≥ 2. According to a standard Riemann-Lebesgue suppression argument, the ratio of the nω harmonic amplitude to the fundamental harmonic central value scales as γ n , with where G is an order 1 geometrical factor, used here as a stand-in for the exact shape of the oscillon. The fact that the geometrical factor G is in the exponent shows that even modest changes in the oscillon's shape can dramatically change its lifespan, emphasizing the importance of accurately resolving the oscillon geometry. Moreover, because the factor ω/ √ m 2 − ω 2 becomes larger as ω approaches m, the differences between potentials will be exaggerated in this limit, while low-frequency oscillons will typically be similar to one another (see figures 11 FIG. 6. A physical model for an oscillon radiating into the third harmonic. The black line represents the background oscillon source φ 3 , while the red lines represent the amplitude of the radiated field. The spherical symmetry of the oscillon imposes boundary conditions at the origin which behave like an optical mirror: an inward propagating spherical wave is reflected at the origin, propagating back outward with the opposite phase. The result is that the oscillon radiation may experience two kinds of self-interference: interference from the physical extent of the source, analogous to diffraction of a laser beam through a finite-width slit, and interference due to the spherical symmetry of the oscillon, represented by the mirror. At certain oscillon frequencies, these two effects conspire to destructively interfere, trapping a nominally free harmonic. and 12 for an example). As a consequence of this growing separation of scales, oscillons whose frequency ω approaches the mass m radiate at increasingly suppressed rates, so that the last phase of the oscillon's life is often the longest. We refer to this general trend as geometric decoupling.

B. Destructive interference and the minimum radiation condition
Throughout the oscillon lifetime, radiative harmonics are subject to self-interference, which is totally destructive at exceptional frequencies. At these points, destructive interference completely confines specific harmonics, and subverts the expected radiation hierarchy implied by geometric decoupling. When the leading harmonic is confined, the overall radiation amplitude shrinks by another global factor of γ. For many especially long-lived oscillons, a period near harmonic confinement dominates the total lifetime. In principle, it is possible to imagine engineering ultra-long-lived oscillons by aligning the destructive interference of multiple harmonics, leading to additional suppression by γ , where is the number of aligned exceptional points. In practice, these constructions are necessarily fine tuned, since each resonance must be aligned to order γ −1 .

Interferometric analogue
The basic physics of oscillon radiation is captured by the physical model in figure 6, which describes an interference experiment reminiscent of the classic Lloyd's mirror. In this simple one-dimensional setup, a coherent, finite-sized, optical source at r > 0, representing the oscillon's coupling to the radiative harmonic, is placed in front of a mirror at r = 0, representing the spherical symmetry of the oscillon. Each point in the source experiences interference both from its reflection, and from its neighbors. Let the spatial location and magnitude of the source be described byJ (r). The direct radiation reaching the observer is therefore On the other hand, the reflected light paths sum up to an amplitude: where, crucially, a half-wavelength path difference is picked up upon reflection at the mirror. This is equivalent to enforcing the usual regularity conditions at the origin in a spherically symmetric field solution. Finally, the observer adds up these contributions coherently, which explicitly leads to an amplitude equal to the sine-transform of the source: In the following section, we derive a similar result from the mode equations of the PQB, and quantify corrections to this simplified picture.

The physical quasibreather picture
In the previous section, we introduced a simple interpretation of the oscillon radiation in terms of the interference of a coherent source with its own reflection. Here we study the mode equations (7), in which the first radiative harmonic S 3 and the orthogonal deformation c 3 are treated as a perturbation of the fundamental S 1 . Under this perturbative assumption, the mode equations for the radiative harmonicS 3 ≡ rS 3 and its orthogonal deformationc 3 ≡ rc 3 further simplify to the frictionless linear systemS Here, k S and k c represent the r-dependent wavenumbers of the third harmonic S 3 and c 3 in the background of the fundamental harmonic S 1 , and rJ 3 (r) corresponds to the 3 → 1 processes generating the radiation. Note, the wavenumber is different for the third harmonic S 3 and its orthogonal deformation c 3 , a distinction explicitly derived in the appendix result (D9). There, we find that the difference between k S and k c appears at sixth order in a Bessel expansion of the background, and therefore is typically small, making the approximation k S = k c quantitatively good in most circumstances. We can solve the linear system analytically in terms of two linearly independent solutions y S 1,2 (r) and y c 1,2 (r) of the homogeneous equations. In this case, the expression for the Green's function is simple and the full solution becomes a sum of homogeneous (defined by initial conditions) and inhomogeneous contributions, of the form: where be the sine-like solution (nonzero derivative at r = 0) and let y S,c 2 be the cosine-like solution (zero derivative at r = 0). Regularity at the origin requires that only sine-like initial conditions are allowed, constraining the cosine-like terms to be zero b H 2 = a H 2 = 0. In the far-field region, all solutions y S,c 1,2 are simple combinations of sines and cosines of frequency k 3 = (3ω) 2 − m 2 . However, orthogonality between y 1 and y 2 is generally not maintained into the far-field. Without loss of generality, we can introduce phase-shifts to express these misalignments, with the understanding that when these phase-shifts are zero, we regain the simple constant-wavenumber Helmholtz solutions. These phase-shifts can in principle be computed in the WKB approximation. Furthermore, we define the orthogonal components y s = sin(k 3 r) and y c = cos(k 3 r) against which we can project the shifted solutions, leading to Substituting, we collect the orthogonal contributions to the radiative tails as where a I 1 and a I 2 are fixed, representing the total inhomogeneous contribution from the oscillon background Radiative boundary conditions (10) match the coefficients of y s and y c between S 3 and c 3 , which uniquely determines the homogeneous degrees of freedom, Consequently, the solution simplifies to In other words, the amplitude of the radiation is always proportional to the inhomogeneous contribution a I 2 . At exceptional frequencies, this contribution is exactly zero and the harmonic experiences totally destructive interference; this is visible in the power versus frequency plots as a sudden drop (see figures 11 and 12 for example). Therefore, in this linear model the condition for totally destructive interference is In the case of a flat wave-number, i.e. Helmholtz system, this is precisely the sine-transform of the source, as predicted by the simple interferometric model. Because totally destructive interference is equivalent to a single constraint on one free parameter ω, we conclude this effect is generic, and not the result of some fine tuning.
To reach this result, we have effectively solved for the physical quasibreather, defined by the choice of a H 1 in (28), at the level of the third harmonic and in a linear approximation. In previous literature (e.g. [54]), a different quasibreather was highlighted as relevant in approximating the oscillon, namely the minimum-radiation quasibreather. This corresponds to a different choice of homogeneous parameters; in this case, the construction of c 3 is irrelevant and the value of a H 1 is chosen such that S 3 is minimized at the level of (25), specifically by picking: We see that this differs from the physical quasibreather answer (28) by an additional a I 2 cos ϕ S 2 tan ϕ c 1 , which is zero in the case when ϕ c 1 = 0, i.e. when the wavenumbers k S (r) and k c (r) are identical functions of r. While typically small, differences between k S (r) and k c (r) appear at higher-orders in the background, and are not guaranteed to be perturbative -as derived below in appendix D 1. Therefore, the minimum-radiation quasibreather and the physical quasibreather are generally close but distinct, and are identical only at the exceptional 'dip' frequency where both predict zero radiative tails.

C. Energetic death
As explained in section III A, the spatial extent of the oscillon increases as it radiates away its binding energy. On the other hand, the balance between self-and binding-energy demands that the oscillon's central amplitude decreases. Depending on the number of spatial dimensions, one effect or the other dominates the oscillon's total energy as ω approaches m. In particular, in three or more spatial dimensions, the volume turns out to grow faster than the central amplitude shrinks. The oscillon's parent PQB also obeys the same scaling relation, and at some point the bound energy in the PQB will necessarily begin to increase. To keep up, the oscillon would need a source of energy; in its absence, the oscillon is forced off the PQB trajectory, in a process we call energetic death.
To make these ideas precise, we can invoke the mode equations (7), in the limit of small central amplitude S 1 (0). Note, because the oscillon's volume is large, it is geometrically decoupled from radiation according to the argument in section III A, and therefore it is safe to neglect backreaction from the radiative harmonics. Keeping only the leading quartic nonlinearity in the potential, S 1 is described by Here, d is the number of spatial dimensions. To extract the scaling of S 1 (0), we match the binding energy of the oscillon to its self-energy, leading to Therefore, the scaling of the central amplitude is independent of dimensions, namely On the other hand, since the spatial extent of the oscillon scales like scale 1/ √ m 2 − ω 2 (as seen in section III A), its volume must increase according to Combining these two scalings results in the oscillon's total energy which grows as ω approaches m for spatial dimension d ≥ 3. In other words, the expectation that the oscillon energy decreases as a function of ω is only true up to a specific frequency strictly less than m. Beyond this point, the oscillon energy is forced to increase as a result of weak self-interaction. Such an increase is unphysical, and the value of ω at which the PQB's energy is minimized sets the moment of death. For an earlier argument along these lines, see [47,54].
For an explicit comparison, take the d = 1 sine-Gordon oscillon, which has a simple analytic form In the ω → m limit, the energy of the sine-Gordon oscillon is exactly 16 √ m 2 − ω 2 , which matches our predicted scaling.
All the examples of oscillons studied in sections IV, V and VI live in three spatial dimensions, and therefore exhibit an energetic death. In other words, for each oscillon there is a specific frequency strictly less than m beyond which the scalar field may no longer exist close to a PQB. After this point, our formalism no longer applies, and the oscillon is considered "dead." Afterwards, gravity may take over leading to the formation of much more diffuse configurations such as axion stars [35,36]. In our numerical simulations, this moment of death is distinctly visible as a "loop," representing the rapid conversion of the oscillon into radiation through 3 → 1 processes (see figures 11, 13, and 14).

IV. A PRESCRIPTION FOR OSCILLON LONGEVITY
Here we provide a procedure for generating potentials that support cosmologically long-lived oscillons. In section III, we explained how the longest-lived oscillons exhibit a combination of geometric decoupling and destructive interference. Geometric decoupling refers to the suppression of radiation when the oscillon size is much larger than the radiation wavelengths, which is especially pronounced at large frequencies ω close to m. For a large oscillon, the interferometric 'fringe pattern' also occurs more rapidly, leading to more instances of destructive interference which further suppresses radiation. Thus, we may find long-lived oscillons by searching for potentials that support large oscillons at frequencies ω close to m. An apparent obstacle to this goal is due to energetic death (see section III C), which limits the frequencies for which the oscillon can have decreasing energy as a function of ω. In the following, we identify a feature in the scalar potential that can stave off energetic death and produce large oscillons.
In section II, we introduced the mode equations (7) obeyed by the radial profiles of the PQB harmonics S n (r), and the sense in which these harmonics may be thought of as the coordinates of a point particle, whose initial condition is tuned so that (S 1 , S 3 , . . . ) = 0 at r = ∞. Here we aim to study the longest-lived oscillons, whose radiation is necessarily small. Moreover, we will focus on large frequencies ω ≈ m, for which higher harmonics are further suppressed by a natural separation of scales ω √ m 2 − ω 2 . Therefore, we will drop the higher harmonics n ≥ 3 in this section's analysis, and work with a simplified 1-dimensional point-particle picture, representing the radial profile of the fundamental mode, S 1 (r).
We now introduce the equations that govern S 1 from first principles, using an effective action technique, equivalent to the PQB formalism for a single bound harmonic. The Lagrangian describing the real scalar φ is Because both V and φ are proportional to f 2 (as in (1)), f is an overall factor in the action, and therefore does not contribute to the dynamics. Hence, for the rest of this section, we will work in units of f = 1. Since we are looking for quasiperiodic, spherically symmetric solutions dominated by the fundamental mode, we substitute φ = S 1 (r) sin ωt and integrate out time, leading to the effective action for S 1 : By integrating out time, we arrive at an action for a point particle S 1 (r), where r acts like a time coordinate, and the resulting effective potential is Finally, the equation of motion for S 1 arising from this effective action carries a 2/r friction term from the spherical Jacobian where V eff represents the derivative of V eff with respect to S 1 . A solution of these equations which describes an oscillon profile needs to respect regularity conditions at r = 0 and r = ∞, corresponding to S (0) = 0 and S(∞) = 0. All solutions which respect regularity at r = ∞ must exponentially decay, since V eff behaves like a quadratic hilltop − 1 2 (m 2 − ω 2 )S 2 1 for small S 1 . From the perspective of the point particle, this means that initial conditions are tuned such that S 1 has just enough energy to climb up the hilltop at 0.
In order to engineer large oscillons, we need the particle S 1 to stay at small velocities so that the oscillon interior spreads out. Initializing on a hillside of V eff is detrimental to this goal, since the slope of V eff controls the speed of S 1 , typically leading to a small oscillon core. On the FIG. 7. Effective potential V eff (S1) for a long-lived oscillon, at three nearby frequencies. The example is obtained using the frustrated quadratic method defined in (44) with m 2 f = 0.9m 2 and b = 2, computed using three Fourier coefficients V1,2,3 with V3 forced to satisfy the mass constraint in (1). We see that as ω passes through the frustrated mass m f , new solutions to the equations of motion (43) emerge, specifically when the local maximum of the effective potential increases to positive values. The balls are placed at the values S1(0) which initialize physical oscillon solutions at the respective frequencies ω. The inset figure shows the trajectories of the smallest-amplitude solutions of (43) for each of the three potentials plotted.
FIG. 8. Lifetime versus frustration for oscillons in frustrated quadratic potentials, computed using three and four Fourier coefficients (see equation (44)). The lifetimes are integrated over the interval ω ∈ [0.8, 0.999] in the one-non-perturbative harmonic PQB formalism. We speculate that introducing more Fourier coefficients leads to longer-lived oscillons, since the frustration mass can be closer to m before self-interactions become repulsive, leading to enhanced geometric decoupling. The line of best fit for three coefficients (dashed purple) is log 10 (mt) = 28(m f /m) 2 − 11, and the best fit with four coefficients (solid blue) is log 10 other hand, releasing S 1 close to a hilltop allows S 1 to remain at low velocities for a time inversely proportional to the initial displacement of S 1 from the hilltop. Therefore, we need to connect hilltop-initialized solutions (i.e. low central velocity) to physical solutions (i.e. which arrive at S 1 (∞) = 0). To compensate for the energy lost by 2/r friction, a physical solution must be initialized with positive potential energy (where we've normalized V eff (0) = 0). Effective potentials which satisfy these two conditions will have non-trivial local maxima, whose hilltop is higher than 0 (see figure 7).
Here we reverse engineer a class of scalar potentials V (φ) which generate effective potentials V eff (S 1 ) with the aforementioned hilltops. This construction makes use of the fact that the mass term in the effective potential V eff (S 1 ) acquires an ω 2 offset compared to the mass term in the scalar potential V (φ). Based on this observation, we introduce the family of frustrated quadratics, whose Fourier coefficients (in the basis expansion (1)) are chosen as the solution to the following optimization problem: where 0 < b < π, and 0 < m f < m is the frustrated mass.
In other words, we are forcing the potential to have mass m at small φ, and a different, smaller mass m f at larger φ. For frequencies ω close to m f , the effective potential V eff will consist of a series of hills and valleys inside the , whose amplitude is controlled by how tightly the objective (44) is optimized. As local hilltops in such potentials rise above the zeropotential line, new oscillon solutions emerge. When the hilltop is precisely at the zero-line, this new solution is wholly unphysical, carrying infinite energy. As ω increases, this hilltop is pushed upwards, and S 1 (0) starts with more potential energy that needs to be dissipated through friction. As a consequence, S 1 (0) starts further from the hilltop so that it begins rolling earlier while the 2/r friction is still active, corresponding an oscillon with a smaller radius and less energy (see the inset of figure 7). Even though this branch may appear at very large ω close to m, this effect guarantees there is some finite range of frequencies over which the energy of these solutions decreases, meaning a physical oscillon can be supported.
In figure 8, we plot the lifetime of oscillons in frustrated quadratics as a function of the frustration m 2 f /m 2 . The frustration mass m f controls the frequency at which new hilltop solutions emerge. As m f increases towards m, the appearance of these new branches occurs at larger frequencies, taking advantage of enhanced geometric decoupling and leading to longer lifetime. Increasing the number of Fourier coefficients in the potential reduces the height of the hilltops in the effective potential, allowing them to emerge at larger frequencies. Further, higher frequencies in the potential pushes the hilltops closer to S 1 = 0, allowing for lower-energy oscillons. We speculate that a fixed number of Fourier coefficients in the potential implies an upper bound on oscillon lifetimes, although we leave this question to future work.

V. IS LONGEVITY FINE-TUNED?
There are many known examples of potentials which support very long-lived oscillons, including those identified in section IV. However, the precise form of these potentials remains largely unconstrained by a UV theory, and therefore it is not clear how to assess whether their longevity is a result of fine-tuning, since the distribution from which the potential coefficients are sampled strongly influences the lifetime. Therefore, we introduce two notions of tuning that attempt to quantify the difficulty of constructing a theory with a long-lived oscillon: 1. Global Tuning asks what fraction of parameter space hosts long-lived oscillons. A typical object of study is the probability distribution of lifetimes, sampled with minimal priors over the potential coefficients in some natural basis.
2. Local Tuning asks whether a given long-lived oscillon is sensitive to variation in its potential parameters. The typical objects of study are the local gradient and curvature of the lifetime with respect to the parameter space at the point in question.
In sections V A and V B, we address the genericity of long-lived oscillons in periodic potentials with parity. The advantage of studying periodic potentials is that they are naturally expanded in the Fourier basis. Without any theoretical priors, a natural scale for the Fourier coefficients is m 2 f 2 , and variations will be of the same size.
One may also be interested in studying the genericity of long-lived oscillons in monodromy potentials [76][77][78].
Since an oscillon has a finite amplitude, one may restrict the aperiodic monodromy potentials to a compact interval, which is fully described by a Fourier expansion. However, any realistic model of axion monodromy is asymptotically a power law, meaning the high frequency modes of the potential are perfectly correlated. To sample the full space of monodromy potentials, one must sample from a distribution that imposes this correlation. In the absence of a reliable way to select coefficients from this distribution, we leave this question to future work. Instead, in section VI A, we scan the one-parameter family of monodromy potentials studied in [50,57,58].

A. Global tuning
In this section, we study the distribution of oscillon lifetimes as a function of the potential coefficients in the Here V3 is constrained such that the mass is fixed to m, with all other V n≥4 = 0. The red region indicates parts of the parameter space where φ = 0 is not a global minimum of the potential, and has significantly shorter lifetimes. The stars indicate potentials for which we have compared our formalism with multiple nonperturbative harmonics to direct numerical simulation (see Figure 11). The peninsula of longevity corresponds to the emergence of a frequency at which the third harmonic experiences totally destructive interference at 'dips.' The yellow banding corresponds to the migration of dips to higher frequencies, where geometric decoupling suppresses the fifth harmonic, increasing the impact of the dip. At the upper right of these bands, the dip migrates to frequencies higher than that of energetic death, creating a longevity cliff.
Fourier basis. In particular, we will consider periodic potentials with parity where the sum of the coefficients V n constrains the mass of φ to be m. In Figure 9, we plot the lifetime versus the free variation of the first two coefficients n = 1, 2 with the third constrained so that the sum in (45) is satisfied, with all other V n set to zero. The mass constraint in (45) naturally sets the typical scale of V n to 1. Therefore, we restrict our study to V n in the range [−1, 1], inspired in part by the fact that the QCD axion potential has order 1 coefficients in this basis (see equation (51)). Figures 9 and 10 provide illustrative examples of some important qualitative features of the distribution of oscillon lifetimes. First, we observe islands of longevity, seen in figure 9 as localized regions of exceptionally long lifetimes. In figure 10, this feature is manifested as plateaus   FIG. 10. The distribution of oscillon lifetimes for 1 (yellow), 2 (gray), and 3 (red) degrees of freedom in a periodic potential. We uniformly sample the nDOF-dimensional cube V1, . . . Vn DOF ∈ [−1, 1] restricting the potential such that φ = 0 is a global minimum, and Vn DOF +1 is fixed such that the mass is m, with the remaining Vn set to 0. Lifetimes are computed in the interval ω/m ∈ [0.8, 0.995] in the single-nonperturbative-harmonic approximation. The geometric suppression of the radiative modes means that these frequencies likely dominate the oscillon lifetime, and that the perturbative radiation approximation is typically good. We see that each new degree of freedom is observed to introduce a new island of longevity (island 1 log 10 mT ∈ [0, 4], island 2 log 10 mT ∈ [4,9], island 3 log 10 mT ∈ [9,14]).
FIG. 11. The power radiated by the oscillons in the potentials denoted by stars in figure 9. The dark curves are data from explicit numerical simulations (see appendix E), while the lighter curves are computed in the PQB formalism. The PQB predictions become dotted in regions of linear instability, as computed using the methods described in appendix C. Notice that at low frequencies, the oscillon power curves are of similar magnitude, diverging at larger frequencies due to geometric decoupling, as explained in section III A. The loops at the end of the simulations correspond to the oscillon rapidly converting into 3ω modes past the point of energetic death, causing the measured frequency at the origin to briefly grow larger than 1. in the distribution of lifetimes. We observe that each successive degree of freedom introduces a new longer-lived island of longevity which we observe to be exponentially more long-lived than the last. The probability of landing on one of these islands decreases with a scaling expected to be exponential in the number of degrees of freedom.
With these observations in mind, we introduce a notion of global tuning based on the cumulative probability of finding an oscillon at least as long-lived. Therefore, smaller values of this probability mean more extreme outliers, and thus higher degrees of global tuning. For example, according to our PQB simulations, summarized by figure 10, an oscillon of lifetime log 10 mT = 3 is tuned to one part in 2 (or 50%). Oscillons of lifetime log 10 mT = 7 are tuned to one part in 8 (or 10%), and oscillons of lifetime log 10 mT = 12 are tuned to one part in 400 (or 0.2%). Finally, the longest lived potential we observe in our search lives roughly log 10 mT = 14, and its tuning is roughly one part in 3000, or 0.03%, although longer-lived oscillons may still reside further up the distributional tail.

B. Local tuning
As opposed to global tuning, which deals with the statistics over large volumes of parameter space, local tuning attempts to quantify the sensitivity of an oscillon's lifetime to variations in its potential coefficients. If we understand the lifetime mT as a function of the potential coefficients V ≡ {V 1 , . . . }, we can naturally introduce a local approximation of mT as a function of its gradient and curvature, writing V = V 0 +δ V, and mT ( V 0 ) = mT 0 , In terms of this local approximation, we may quantify the sensitivity of the lifetime to local variations in V, as the minimum relative displacement of the potential coefficients ||δ V||/|| V 0 || necessary to change the lifetime by an order of magnitude log 10 mT 0 ± 1. In other words, our local tuning metric is the solution to the following constrained optimization problem, We denote the minimal value ν ≡ ||δ V||/|| V 0 ||. Consider the potential V = (1, 1/2, −1), for which an oscillon lives approximately log 10 mT = 14. Using the above measure of tuning, and a grid based approximation to the gradient and Hessian, we calculate ν ≈ 0.03. In other words, a 3% variation in the potential parameters corresponds to an order of magnitude change in the lifetime of the oscillon. This is substantially less tuned than one would expect from our global metric, in which this potential is 0.03% tuned. This is a reflection of the structure of the lifetime landscape, which contains islands of stability seen in figures 1a, 9 and 10.

VI. ILLUSTRATIVE EXAMPLES
Here we apply our framework to a series of potentials that have been studied extensively in previous literature, with the aim to reproduce and expand upon known results. The main goal is to show how our methods can accommodate a wide variety of potentials: both with or without parity, and with or without periodicity. We compare the results of our PQB framework to explicit numerical simulation. When simulating the equations of motion explicitly (as in appendix E), the wall-clock time is at least proportional to the oscillon lifetime, which becomes computationally prohibitive for lifetimes beyond 10 10 /m. Our framework can bypass this scaling since time has been explicitly integrated out, allowing us to predict the existence of very long-lived oscillons, well in excess of the lifetimes we can simulate explicitly.
The results of this section are presented in the form of "power vs frequency" curves, which represent the instantaneous flux radiated by the oscillon at a particular fundamental frequency ω. The oscillon fundamental frequency ω monotonically increases with time, and can therefore be thought of as a time coordinate. For a detailed review of how to interpret these plots and their features, see figure 2.

A. Axion monodromy
Monodromy refers to the non-trivial winding of an axionic degree of freedom, which effectively endows the axion with an aperiodic potential at low energies [7,76,77]. In general, the resulting monodromy potentials share the property that they asymptote to a power law at large field values. A common family of potentials which interpolate between the asymptotic power-law and the smallfield mass m is where p scans the asymptotic power-law. The potential (48) has been widely studied, and has been shown to support very long-lived oscillons, in excess of 10 10 cycles [50,52,57,58].
In figure 12, we summarize our study of the oscillon life-cycles as we scan p from −8 to −1. In general, we find good agreement with the results of [58]: the power versus frequency curves and lifetimes broadly match the predictions of [58] in the cases we have mutually studied p = −8, −5, −4, −1 although there are minor differences.  [58]. Here, the power P and binding energy EB are computed as in section II. As p scans from −8 to −1, the third harmonic dip migrates to larger frequencies where the fifth harmonic is further suppressed by geometric decoupling, leading to increased lifetime. To obtain the PQB results, we start with a two-non-perturbative harmonic approximation and used three non-perturbative harmonics to obtain better accuracy near the dips. At frequencies below the dip frequency, we see a small shift in the PQB formalism vs [58], which may arise from the need to use more non-perturbative harmonics at lower frequencies or because the Fourier series representation of (48) converges slowly.
As p increases from −8 to −1, the lifetime of the corresponding oscillon increases dramatically, from 10 6 to 10 10 cycles. This is due to the simultaneous action of the two longevity mechanisms identified in section III. Specifically, the third harmonic experiences totally destructive interference at an exceptional frequency that is larger with increasing p. Therefore, as p grows, the third harmonic dip moves deeper into the frequencies where geometric decoupling dominates, which further suppresses the fifth harmonic.
A natural conjecture is that the oscillons of (48) are unusually long-lived because of the asymptotic power law in the potential. However, our results in figure 12 indicate that longevity is dominated by large frequencies, where field amplitudes are too small for the asymptotic behavior to take over. In particular, for p = −1, the field amplitude at the origin is roughly φ(0) ≈ 1.5f , far too small to be sensitive to the flatness of the potential at large φ/f . Therefore, we conjecture that it is not the asymptotic form, but the details of the connection between 1 2 m 2 φ 2 and φ p that determine the oscillon lifetime.
In the oscillon literature, many examples of extremely long-lived oscillons are obtained with monodromy potentials. Thus, a natural question is whether all monodromy potentials share a common feature leading to longevity, or whether simple examples such as (48) happen to live in a tuned island of longevity. In the language we introduced in section V A, in order to quantify the link between monodromy and longevity we would need to know the probability distribution from which monodromy potentials are chosen. In the absence of this non-trivial construction, we are left with a case-by-case analysis of particular potentials, which, in this probabilistic view, may suffer from sampling bias.
B. φ 4 theory φ 4 theory is the quintessential example of spontaneously broken parity symmetry. It is well known to host moderately long-lived oscillons, which have been studied in previous work [37,44,51,52,55]. Shifting to the broken vacuum and fixing the mass of φ = θ/f to be m, we arrive at the following parity-violating potential In order to properly compute the physical quasibreather in this potential, the first three harmonics C 0 , S 1 , C 2 must be treated non-perturbatively. As is evident from the numerical simulation (figure 13), the φ 4 physical quasibreathers are linearly unstable over the entire range of ω for which the oscillon is long lived. The instability that occurs at linear order in the PQB background is, however, stabilized by self-interaction, leading to quasiperiodic oscillations. These nonlinear oscillations are visible as dense curly-Q's in the Power vs Frequency plot ( figure  13).
FIG. 14. Comparison of the explicitly simulated QCD axion oscillon (black) to the PQB formalism (red) with three non-perturbative harmonics. The radiated power is so large that the orthogonal deformation is non-perturbative, leading to disagreement within a factor of a few, although the shape of the physical quasibreather curve still captures the qualitative features of the simulated oscillon. Namely, it shows that there is a dip around ω = 0.6m where the fifth harmonic is dominant. This region, in which the third harmonic is confined and non-perturbatively large, constitutes most of the oscillon's lifetime.
Our explicit numerical simulation yields an approximate lifetime of 6000/m, which is close to the PQB prediction of 5900/m. This confirms the earlier results in [37].

C. The QCD axion
The QCD axion is the best studied example of a scalar field described by a periodic potential, and could allow for oscillons with cosmological observables. At low temperatures, the QCD axion potential is dominated by strong dynamics, giving rise to the potential [14,36], A simple Taylor expansion about θ = 0 reveals that (50) has a smaller quartic term than the simple cosine potential, which leads us to expect the oscillons in (50) to be shorter-lived than the sine-Gordon oscillon. Indeed, this expectation is confirmed by the physical quasibreather framework and explicit numerical simulation (see figure  14 and appendix E).
In order to compute the lifetime in the physical quasibreather formalism, we calculate the potential's Fourier coefficients Our formalism can accommodate many non-zero Fourier coefficients, although only the first three (and a fourth to normalize the mass) are necessary to converge within 1% of the true potential; adding more terms has a negligible impact on the results of explicit simulation and on the result of our PQB formalism. The result of this analysis is that the QCD axion oscillon is relatively short-lived compared to the average oscillon, living about 400/m. However, it is interesting to observe that this oscillon spends most of its life with a confined third harmonic, undergoing very large centralamplitude oscillations of order 15f . Although the short lifespan of the QCD axion oscillon means that it will only leave its cosmological imprint shortly after formation, the large amplitude and violent deaths of these oscillons may have observational implications.

VII. CONCLUSION
Real scalar fields play a central role in many theories of early universe cosmology and dark matter. Many of these theories predict attractive self-interactions that allow the scalars to form quasistable oscillons. Understanding oscillon lifetime is necessary for determining whether oscillons only play a role in early universe cosmology, or whether they may also survive until the present day and lead to dramatic astrophysical signatures.
In this work, we have expanded the quasibreather approximation into a formalism for computing the properties of oscillons that naturally incorporates realistic boundary conditions. We defined the physical quasibreather by finding initial conditions of the nonlinear wave equation that simultaneously obey radiative boundary conditions and specify a quasibreather solution. As the closest quasibreather to a physical oscillon, the PQB provides a raw approximation for the oscillon profile (see e.g. figure 4) which is increasingly accurate in the limit of long lifetimes. Furthermore, the PQB represents the solutions to which the oscillon is instantaneously and locally attracted to during its evolution. When understood as a stable perturbation to a PQB, the oscillon borrows its properties from its PQB partner, including its radial profile and radiation rate [79]. When the PQB becomes linearly unstable, nonlinear quasiperiodic oscillations emerge, whose size is controlled by higher-order terms, as depicted in figures 13, 18a, and 18b. In other words, linear instability often does not result in the death of the oscillon. Further, we have demonstrated that the PQB and the minimum radiation quasibreather differ at 6th order in the background, explaining the success of the 'minimum radiation quasibreather ansatz. ' Since the PQB offers an accurate description of the oscillon structure, we have used it to understand the oscillon's form-factor and the resulting mechanisms which control longevity. Specifically, as the oscillon radiates its energy away, its central amplitude decreases, causing self-interactions to become weaker; as a result, the oscillon becomes much larger than the radiation wavelengths, suppressing the radiated power (see figure 2). At these high frequencies (ω → m), the large oscillon core naturally leads to the rapid self-interference of radiation. When the self-interference is destructive, the total radiation is suppressed by another power of the formfactor. While both these effects, geometric decoupling and destructive interference, are generic features of oscillon evolution, the longest-lived oscillons are a consequence of these two effects occurring simultaneously at large oscillon frequencies close to m (as in figure 12). Finally, we have understood the physics of oscillon death as a further consequence of weak self-interaction: past a certain critical frequency the energy of the PQB is forced to increase, and the oscillon cannot sustain its proximity to its PQB partner. Using our understanding of the mechanisms responsible for oscillon longevity and death, we have constructed the family of frustrated quadratic potentials which support extremely long-lived oscillons, living more than 10 18 cycles (see section IV).
There are several computational advantages provided by our methodology. First, the oscillon evolution is computed in a time-independent way, separating the physical lifetime of the oscillon from the computational wall-clock time it takes to evolve numerically. Second, our formalism naturally incorporates non-perturbative harmonics, and potentials without even parity. Third, all perturbative harmonics may be efficiently computed by taking advantage of sparse linear algebra. Fourth, we have introduced the Fourier basis as the natural basis for expanding potentials when studying oscillons. In this basis, the form of the mode equations is especially simple, allowing us to write down analytical expansions for the mode potentials that converge everywhere. Fifth, our formalism provides a natural language for studying the stability of oscillons. Finally, the speed of our numerical techniques has enabled us to study extremely long-lived oscillons, and has yielded the first prediction of cosmologically long lived periodic potentials (see figure 8).
Using our efficient numerical techniques, we scan over degrees of freedom in axionic potentials (see figures 1a and 9), allowing us to probe the genericity of long-lived oscillons (see figure 10). Important outcomes of this parameter scan include the identification of features in the lifetime landscape with the mechanisms of longevity (see section III), and the realization that extremal lifetimes may scale at least exponentially in the number of potential degrees of freedom. At the same time, long lifetimes are not particularly fine-tuned, since as few as three degrees of freedom are enough to generate oscillons that survive until last-scattering (4 × 10 6 cycles) with only 15% global tuning (as defined in section V A), and oscillons that live until today (10 11 cycles) with 0.5% global tuning.

ACKNOWLEDGMENTS
We thank Savas Dimopoulos, Asimina Arvanitaki, Sebastian Baum, Robert Lasenby, Horng-Sheng Chia, Anirudh Prabhu, and Dmitriy Zhigunov for stimulating conversations and helpful comments on the manuscript. D.C. is grateful for the support of the National Science Foundation under grant number PHY-1720397, and for the support of the SITP. T.G.T. is grateful for the support of the Department of Energy under grant number DE-SC0020266. Some of the computing for this project was performed on the Sherlock and Farmshare clusters. We would like to thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that contributed to these research results.

Appendix A: The physical quasibreather formalism
In this appendix, we outline the precise definition of the physical quasibreather (PQB) and its orthogonal deformation (OD) introduced in section II. We will see that the orthogonally deformed PQB is an instantaneous solution of the equations of motion that satisfies outgoing boundary conditions. We then find oscillonic solutions of the equations of motion that are perturbations of the orthogonally deformed PQB. By studying the evolution and stability of these perturbations, we arrive at a sense in which the orthogonally deformed PQB can be an attractor, which we apply to study oscillon stability in appendix C.

Quasibreathers
Physical potentials may be interpreted in terms of nparticle interactions, and therefore possess Taylor expansions around their vacua. Consequently, a periodic field configuration with fundamental frequency ω will only couple to modes oscillating with integer multiples of this fundamental frequency. In other words, physical nonlinear wave equations possess periodic orbits, which may be interpreted as a Fourier series in time. This is in contrast to unphysical potentials that may not be interpreted in terms of integer-number particle interactions, which can at best possess quasiperiodic orbits.
For the remainder of the appendix, we move into dimensionless units with m = f = 1. The non-linear wave equation for the field θ in a potential V is then 0 =θ − ∇ 2 θ + V (θ) . (A1) As we have argued above, V must possess a Taylor series, and therefore θ may be expanded as a series of integer harmonics where δ n is a phase, and we have taken spherical symmetry for simplicity. Without loss of generality, we may take δ 1 = 0. We say that a solution of the form (A2) is generated by the frequency ω if S 1 is non-zero, and the only non-zero higher harmonics S n are those that couple to S 1 , consistent with closure of the equations of motion. We then define the quasibreather as the solution generated by ω.
Using this definition, we may compute the generic form of the quasibreather. Consider the generic potential From the symmetries of sine and cosine, we observe that (sin ωt) n = k a k sin n k ωt, n k ∈ N odd , n is odd, k b k cos m k ωt, m k ∈ N even , n is even, (cos ωt) n = k c k cos n k ωt, n k ∈ N odd , n is odd, k d k cos m k ωt, m k ∈ N even , n is even.
The case of parity V (θ) = V (−θ) offers a pleasant simplification, decopling the even harmonics and the odd harmonics from one another. Thus, potentials with parity have quasibreathers of the form and a periodic solution of the form θ = n∈Neven S n (r, ω) sin nωt , although it is not a quasibreather because it is not generated by ω. Quasibreathers in potentials without parity possess expansions where N even contains 0. Thus, we have identified the form of the quasibreathers of the non-linear wave equation when V represents a physical interaction.
Inserting the form (A5) into (A1), one arrives at the set of mode equations where d is number of spatial dimensions, and V n is defined through the equation (A7) Equation (A6) is a system of second order ordinary differential equations, and therefore each degree of freedom S n , C n must be constrained by two boundary conditions.
In order to discuss boundary conditions, we define the number n 0 as the least integer such that n 0 ω > 1, so that bound harmonics have n < n 0 and radiative harmonics have n ≥ n 0 . Regularity at the origin places a non-trivial constraint on all harmonics, that all S n and C n must have zero first derivative at the origin. However, regularity at spatial infinity is only a constraint on the bound modes, n < n 0 ; all radiative harmonics decay geometrically as they propagate to spatial infinity. Thus, for a quasibreather, the radiative harmonics are only constrained by regularity at the origin, and the space of possible quasibreathers has dimension equal to the number of radiative modes. In other words, one has the freedom to pick the amplitude of the radiative modes at the origin, and the result will still be a quasibreather. The authors of [54] have alleviated this ambiguity by picking a specific quasibreather out of this manifold: the minimum radiation quasibreather, whose radiative tails are the smallest. Instead, we pick the physical quasibreather (PQB), defined below, which is perturbatively close to a radiating solution.

The deformed mode equations
A localized field configuration with a finite lifetime necessarily radiates its energy to spatial infinity, and therefore satisfies radiative boundary conditions at spatial infinity. In this section, we introduce the concept of the physical quasibreather (PQB), which is, in a precise sense, the quasibreather closest to a physical configuration satisfying radiative boundary conditions. First, we define the orthogonal deformation (OD) of the quasibreather (A5), which consists of adding 90 • outof-phase components s n and c n to the radiative harmonics in order to satisfy radiative boundary conditions θ = n∈N odd S n (t, r, ω) sin nωt + n∈Neven C n (t, r, ω) cos nωt n∈N ≥n 0 even s n (t, r, ω) sin nωt + n∈N ≥n 0 odd c n (t, r, ω) cos nωt .
(A8) Notice that here we've introduced a time dependence to the modes, which accounts for the fact that a radiating solution cannot have a stationary profile. This formulation will be useful for studying initial conditions of interest, namely, those which specify a quasibreather and orthogonal deformation that together satisfy outgoing boundary conditions. Although (A8) is a vast overparametrization of a single field, we recognize that a so-lution of the deformed mode equations (A9) is also a solution of the full equation of motion (A1). Equation (A9) is obtained from the equations of motion (A1) by substituting (A8) and collecting the terms proportional to sin nωt or cos nωt for a given n, setting them independently to zero. Intuitively, when the time dependence of the harmonic functions S, C, s, c is slow, they have the usual interpretation as the profiles of quasistationary modes. Here the functions V X = V X (S n , C n , s n , c n ), are the mode potentials, in which we have suppressed functional dependence for brevity. The mode potentials are defined by the equation where θ is written in the form of equation (A8), and V Sn , V Cn , V sn , V cn are pure functions of S n , C n , s n , c n . During the evolution of equation (A9), ω is treated as a constant. This is not in contradiction to the usual understanding that the fundamental frequency of the oscillon increases with time. For the purpose of the mode equations (A9), ω is understood as a choice of a fixed parameter, independent of the time variation of the modes S, C, s, c themselves. For certain initial conditions, and for certain choices of ω, there will be periods of time over which S, C, s, c vary slowly, and it is during these periods that ω may be interpreted as the instantaneous frequency of the oscillon.
In other words, there is no a priori reason to choose a particular ω for a particular field configuration, and one may only think of ω as an instantaneous frequency in the context of certain initial conditions. Thus, the following paragraphs are dedicated to specifying initial conditions which allow ω to be interpreted as the instantaneous frequency of an oscillon, where the oscillon is perturbatively close to a quasibreather. The smaller the orthogonal deformation, the better this interpretation is, and the longer it holds. In this sense, ω may be conceptualized as an adiabatic parameter, although one should not confuse it with an externally controlled parameterin our framework, it is a constant that parametrizes the decomposition (A8) of solutions to (A1).
We now specify the following consistent set of initial and boundary conditions, in which we treat s n and c n as linear perturbations.
Here, we take V X = V X (S n (0, r), C n (0, r), 0, 0), and we define δV X ≡ n≥n0 s n (0, r)∂ sn V X + c n (0, r)∂ cn V X (note the absence of a constant term in δV X is a consequence of (a), below). A complete and consistent set of initial and boundary conditions associated with (A9), that exactly specify a quasibreather and orthogonal deformation at t = 0 is (a) Initial Quasibreather: (c) Maximally stationary: Our initial condition (a) selects S n and C n which specify a quasibreather. This quasibreather is one which may be orthogonally deformed to satisfy radiative boundary conditions, and it is this quasibreather which we call the PQB. Because we have broken the time translation symmetry of the quasibreather by satisfying radiative boundary conditions, the modes S, C, s, c are endowed with an irreducible time dependence. The maximally stationary condition (c) shows that this time dependence is proportional to the pointwise small deformations s and c. Since s and c obey a homogeneous system of equations (b), their amplitude everywhere must uniformly go to zero as their amplitude at r = ∞ goes to zero. From (e), we see that c n ∝ S n and s n ∝ C n at spatial infinity. Thus we conclude that if S n and C n possess small radiative tails, then s n and c n become pointwise small everywhere, and the time variation of the modes uniformly approaches zero. This is the limit in which the oscillon is long-lived, and approaches the quasibreather; in this same limit, the interval over which this approximation is FIG. 15. The asymptotic attractor (red) is approached as the inhomogeneous solution goes to zero. The homogeneous terms, representing the initial conditions at t = −t0 cannot converge exactly to zero by the time the inhomogeneous solution passes through zero, and therefore the perturbation never exactly reaches the asymptotic attractor.
valid, during which ω may be thought of as an instantaneous frequency, becomes longer.

The asymptotic attractor
The initial conditions (a-e) specify a solution to the nonlinear wave equation that, at t = 0 is exactly an orthogonally deformed physical quasibreather (to linear order in the deformation). More practically, we want to understand the evolution of the oscillon in the neighborhood of this deformed physical quasibreather, before and after this particular point. To this end, we introduce the perturbation δθ(t, r), which simultaneously absorbs the time dependence of the modes in (a-e), and deviations from the orthogonally deformed physical quasibreather. Specifically, a field configuration θ describing a physical oscillon can be expanded as Crucially, θ PQB is exactly a quasibreather solution, and θ OD is exactly periodic in time, as opposed to δθ, which characterizes the secular evolution of the oscillon in the vicinity of the physical quasibreather at ω. Inserting . Each ellipse centered on a deformed PQB represents the domain of frequencies and field values over which that specific quasibreather is an asymptotic attractor. As the oscillon trajectory enters an attractive region, it moves closer to the attractive deformed PQB. Consequently, it is also drawn into the attractive vicinity of the neighboring PQBs. Therefore, the oscillon is forced to approach the red trajectory as the radii of attraction get larger and larger towards the bottom of the dip. After traversing the dip, the deformed PQB radii of attraction begin to shrink, and the oscillon trajectory begins to diverge from the deformed PQB trajectory. In this latter half of the evolution, we see how the deformed PQB trajectory does not act as a standard attractor, but can still be described as an asymptotic attractor. To see this, notice how the oscillon instantaneously moves closer to the quasibreather when entering each new attractive bubble.
This is a sourced equation, representing the fact that the physical quasibreather with orthogonal deformation does not conserve energy on its own. As a linear equation, δθ may be decomposed into a sum of homogeneous terms, which obey the homogeneous equation and one particular solution δθ P , that obeys the sourced equation (A14), which we take to be identically zero at t = 0. In the absence of homogeneous terms, it is this particular solution δθ P which satisfies the initial conditions (a-e). Therefore, the homogeneous terms represent perturbations around those initial conditions. If the homogeneous solutions of (A14) are stable, then we say that θ PQB + θ OD is an asymptotic attractor.
The usefulness of the construction δθ is that it contains all information about the linear stability of the oscillon [80]. Just like in a standard damped oscillator, linear stability represents an exponential approach to the inhomogeneous solution. In other words, it is enough to study the stability of the homogeneous equation (A15) with the tools of Floquet theory. The full picture of how the oneparameter family of deformed physical quasibreathers, parametrized by the frequency ω, acts like an attractor may be understood in the following picture. Before t = 0, the particular solution δθ P is approaching 0 (see Figure  15). If the homogeneous terms are stable, then the field θ is approaching the deformed physical quasibreather at frequency ω. However, past t = 0, δθ P begins to grow again, causing the field to diverge from this temporary quasibreather partner. This story repeats by choosing the next physical quasibreather to expand around at a nearby frequency ω + dω, such that the attractive region of this new quasibreather has some overlap with the repulsive region of the previous quasibreather at ω (see figure 16). The term "asymptotic attractor" is chosen because of its likeness to the concept of asymptotic series, in which increasing the order of an expansion increases its precision until, at some point, it begins to diverge.

Energetic instability
The physical quasibreather background around which we expand the perturbation δθ is one among a continuum of quasibreathers, parametrized by their fundamental frequency ω. Thus, when we talk about a perturbation δθ, we introduce the notation δθ ω in order to talk about "the perturbation relative to (the deformed physical quasibreather of frequency) ω," where we may omit the parenthetical when it is unambiguous to do so.
In the previous section, we introduced the concept of asymptotic attraction, in which an oscillon may be viewed as approaching a physical quasibreather for a finite period of time. For each quasibreather, there is an epoch of attraction, during which the particular solution δθ P is shrinking towards zero, and an epoch of repulsion, during which δθ P is growing away from zero. Neighboring quasibreathers at ω and ω + dω have particular solutions that cross zero at different absolute times t = t ω and t = t ω+dω respectively; whether t ω < t ω+dω determines whether ω + dω is attractive for some time after ω becomes repulsive. Because the particular solution δθ P encodes the energy flow out of the oscillon, the relative timing of the zero crossings of δθ ω P and δθ ω+dω P may be viewed as a reflection of energy conservation, defining an arrow of time. That is, the oscillon at ω + dω is energetically accessible from ω if t ω < t ω+dω .
This time ordering implies the existence of a relative energy function, whose local monotonicity encodes whether ω + dω is accessible from ω. In other words, the physical quasibreather at ω + dω is energetically accessible from ω if there is a time when δθ ω is a negative energy perturbation relative to ω and δθ ω+dω is a positive energy perturbation relative to ω + dω. However, the energy of the total field configuration is ill-behaved because of the divergent radiative tails. Strictly speaking, because the quasibreather at ω has a different amplitude radiative tail than the quasibreather at ω + dω, δθ cannot be a finite energy perturbation of both quasibreathers. However, the tails are decoupled and do not influence the dynamics of the oscillon bulk. Therefore, our measure of the perturbation energy must be agnostic to the radiation tails.
One might be inclined to count only the energy inside some finite box containing the oscillon bulk. However, such a measure still grows polynomially with the size of the box. One may also try to subtract the radiative tails by removing the 1/r (in d = 3) asymptotic, although again, this depends on an explicit cutoff between the bulk and the tails. Our framework provides a natural resolution to this ambiguity. Specifically, the orthogonal deformation θ OD provides a measure of the radiative tail of θ PQB valid everywhere. It is the energy associated with this orthogonal deformation that we subtract, leading to our definition of the bound energy in the PQB Note, because θ OD and θ PQB are out of phase, this difference will oscillate around an average value, reflecting the uncertainty principle. This definition has the virtue of converging to the deformed physical quasibreather energy when the oscillon is infinitely long-lived, i.e. when all harmonics are confined.
Having provided an unambiguous measure of the bound energy of the physical quasibreather, we may now address the question of when the perturbation δθ may flow between nearby quasibreathers. Because the frequency of the quasibreather, ω, is a decreasing function of binding energy, under the assumption that binding energy decreases as total energy decreases, a family of physical quasibreathers is energetically stable where Inside a region of asymptotic attraction, the radiation power P of the physical quasibreathers is a good approximation for the radiation power of the oscillon. This leads to a standard approximation for the oscillon lifetime, under the assumption that the perturbation δθ may completely relax (the adiabatic assumption) As the oscillon becomes increasingly long-lived, and thus approaches the physical quasibreather, this prediction becomes increasingly precise.

Appendix B: Quasibreather numerical methods
In the previous section, we arrive at the physical quasibreather as the main object of study which may be used to derive the properties of oscillons in a physical potential. In this section, we develop the numerical tools which enable the efficient calculations of physical quasibreathers and their orthogonal deformations.

Linear radiation
Let us begin by supposing that S n and C n are known for n < n pert , and that the remaining S n and C n are perturbatively small everywhere, so that they obey linear equations. Define the perturbation vector and the deformation vector respectively the diagonal matrix of frequencies the source vector J (S 1 , C 2 , . . . ) and the mass matrices V C (S 1 , C 2 , . . . ), V s (S 1 , C 2 , . . . ), which are functions of the non-perturbative harmonics. Finally, we define the Sommerfeld operator S, which together with the Dirichlet-Neumann 1D flat Laplacian acting on each diagonal block ∇ 2 [81], contains the Sommerfeld radiation condition (e), provided in appendix D 1. In this notation, the equations of motion for the perturbation and the deformation can be written as a sparse linear system This form, in which C and s only couple through the boundary term S, is guaranteed because C on its own solves the equations of motion, and hence any backreaction from a perturbation s must come at second order. The explicit forms of S, V C , V s and J C , J s are provided for several cases of interest in the appendix D 1. Note, J s is proportional to the orthogonal deformation of the non-perturbative modes, and is therefore zero when all radiative modes are perturbative. The fact that we may write the equations for the perturbative modes as a well-determined system of equations is a reflection of the fact that the radiative boundary conditions and regularity conditions completely (and uniquely, for the linear modes) specify the physical quasibreather.

Nonlinear harmonics
The perturbative method in the previous section amounts to solving a sparse linear system, a process that is computationally efficient. Thus, given the knowledge of the non-perturbative background harmonics, we can compute the contribution of arbitrarily many additional harmonics at almost no computational cost. Now we must do the dirty work of computing the nonlinear harmonics. Computing n pert − 1 non-perturbative harmonics will amount to shooting a point particle in n pert − 1 dimensions, and tuning its initial condition so that it lands on the saddle-top at the origin.
The physical quasibreather only feels the orthogonal deformation at second order, and therefore we may use the following procedure to compute the deformed PQB to leading order: 1. Choose S n (0) and C n (0) from the n pert − 1 dimensional space of initial conditions.
2. Shoot the harmonics S n , C n from r = 0 to some large finite radius r out by evolving the mode equations (a).
3. Use the radiative boundary conditions (e) to convert the S n to c n and C n to s n at r = r out .
4. Shoot the c n , s n back to the origin in the background resulting from step 2.

5.
Check for regularity at the origin for the c n and s n , and regularity at r out for the bound harmonics. If not regular, adjust S n (0) and C n (0) and repeat from step 2. Since regularity is equivalent to minimizing the first derivative, this can be implemented by a variant of a binary search procedure, e.g. golden section search.
6. Compute any number of perturbative harmonics using the procedure from the previous section in the background of the non-perturbative harmonics.
7. To account for linear backreaction of the perturbative harmonics, re-shoot the non-perturbative harmonics in the background of the perturbative harmonics. This last step is repeated until convergence.
In practice, it is helpful to break down the n pert − 1 dimensional search into n pert − 1 linear searches that are performed hierarchically. The process of nonlinear shooting is sped up by precomputing the potential functions and using table-lookup. This kind of optimization is especially important when dealing with periodic potentials where repeatedly computing Bessel functions is costly.

Branching of the fundamental mode
In section II, we reduce the problem of finding the radial profile of the oscillon to a classical-mechanical 'shoot- The plot shows the effective potential VS 1 (S1) as a function of S1 for positive values of S1; since the potential is parity-symmetric, the S1 < 0 region is the mirror opposite with respect to the S1 = 0 axis. We have adjusted the vertical axis to better illustrate the qualitative features of the potential. Different regions are colored according to the sign of S1(∞) when launched from that location. A shooting solution is represented by a point on the boundary between a black and magenta region. Whereas initially there was only one zero-node shooting solution (marked by the circle), the new potential adds two more zero-node solutions, marked by the stars. Intuitively, the higher the starting point, the further the particle will travel, causing successive solutions to have an increasing number of nodes. However, the combination of 2/r friction and nonlinearities in the potential breaks this intuition. Depending on the potential's convexity at the initial point, the oscillon may lose a widely variable amount of energy to friction. Therefore, it is at these regions of varying curvature that we expect these new solutions to emerge.
ing' problem. In its simplest case of one non-perturbative harmonic S 1 , the problem further simplifies to the rolling of a massive ball in a double-welled potential V S1 (S 1 ) in the presence of 2/r friction. A shooting solution is one which starts at rest at an initial displacement S 1 (0) and ends at S 1 = 0 at r = ∞.
In linear equations, such as the radial hydrogen atom problem, there is exactly one solution for each integer number of nodes (i.e. zero-crossings) of the radial profile S 1 . In our PQB mode equations, strong nonlinearities break this intuition, as depicted in figure 17. Specifically, a small change in the oscillon frequency ω can change the number of solutions with zero nodes by an increment of two, introducing new branches of oscillon solutions (or eliminating them) when the potential possesses non-trivial convexity. If this new branch consists of quasibreathers with lower binding energy, then the original branch may jump to the low-energy branch after the original branch experiences energetic death. In the reverse scenario, oscillons may form on the high-energy branch but the low-energy branch is energetically forbidden from reaching the high-energy branch.
In the oscillons that we study in e.g. figure 1a, many of the longer-lived potentials contain a high-energy branch of very large, low-amplitude oscillons which only exists in a small range of frequencies close to m. One such example is shown in figure 17. All the examples studied in section IV are the result of purposefully introducing these branches at a specific frequency ω ≈ m f .

Appendix C: Floquet analysis
In appendix A we introduced the notion of asymptotic attraction to describe physical oscillons as perturbations of PQBs. From this expansion, we have reduced the problem of oscillon stability to the study of the linear stability of equation (A15). Standard Floquet theory tells us that the result of this analysis can have two outcomes: the equation is linearly unstable, or it possesses oscillatory states exclusively (modulo boundary effects). In other words, the existence of a stable decaying mode implies the existence of a growing mode, and stability must emerge at higher order in perturbation theory, if at all. Here, we address the linear stability of perturbations δθ, and later argue that nonlinear terms stabilize linearly oscillatory modes.

Linear stability analysis
Let us begin by reproducing equation (A15) for ease of reference: the linearized homogeneous equation for the perturbation δθ in the background of θ PQB is, Recall that θ PQB is a periodic solution of the equations of motion, and therefore can induce parametric resonances. Substituting in the form of the quasibreather (A5), we find which, under parity symmetry of V , further simplifies to where V m is defined by (C2). We will leave specific formulae for V m to appendix D 3. Since (A15) is linear, we may Fourier transform t → Ω, and decompose δθ in spherical harmonics. Because the quasibreather background is periodic, it induces couplings between frequencies separated by integer multiples of the fundamental frequency ω. Therefore, let us restrict our analysis to the values of the Fourier transform δθ(Ω, r) at the discrete tower of harmonics defined as Ω n ≡ Ω 0 + nω, n ∈ Z, where the base frequency Ω 0 can be assumed to lie in the interval (0, ω).
Therefore, the Fourier components on this tower, denoted δθ n (Ω 0 , r) ≡ δθ(Ω n , r), will respect a matrix-differential equation: where is the angular momentum number. This is apparently a quadratic eigenvalue problem in the fundamental frequency Ω 0 [82], although as we will see, it becomes an irrational eigenvalue problem upon imposing transparent boundary conditions [83,84]. The eigenvalue solutions Ω 0 characterize the stability or instability of the system: real eigenvalues corresponding to oscillatory motion; if solutions pick up an imaginary part, the mode will be exponentially growing (if Im(Ω 0 ) < 0) or exponentially decaying (if Im(Ω 0 ) > 0). In the absence of transparent boundary conditions, the solutions come in pairs of complex conjugates; in this closed-box scenario, the existence of a stable (i.e. decaying) mode implies the existence of an unstable mode. An instructive example We may gain some insight into the eigenvalues Ω 0 by studying the simpler case of perturbations inside a box for a potential with parity. The matrix differential equation simplifies such that only the sum over even terms in (C4) survives. At leading order we only include the first harmonic's n = ±1 terms; the reason is that large-n harmonics both decouple from the fundamental and become unbound. This allows us to keep only the V 0 and V 2 terms in the equation. Moreover, we make the assumption that Ω 0 is small compared to ω, representing solutions to the perturbation equations with a separation between the fast and slow timescales; this is relevant when we focus on the boundary between periodicity and instability, where Ω 0 will be small in magnitude. This assumption will be supported by the result of the analysis.
The result is the following 2 × 2 matrix-differential equation where we have suppressed the argument of V 0 , V 2 for brevity and σ i are the Pauli matrices with entries of magnitude 1. We may phrase this as a typical eigenvalue problem in Ω 0 by multiplying through with σ z /2ω, leading to where H and A are Hermitian and anti-Hermitian ma- 18a. The Lyapunov characteristic exponent (the eigenvalue Ω0 of (C4) with maximum imaginary part) for the sine-Gordon deformed physical quasibreather (with an error of ±0.005). The perturbation δθ becomes linearly unstable at ω ≈ 0.88. The nearest asymptotically attractive quasibreather is always finitely far away from the oscillon. When ω > 0.88, the linearly unstable mode is therefore always excited, leading to growing quasiperiodic oscillations on top of the deformed quasibreather background (see figure 18b). Note, throughout this band of linear instability, the mass energy dV 1 4 m 2 S 2 1 is monotonically decreasing, in contradiction with [49]. On the plot, we denote the energetic death at ω ≈ 0.94, where the oscillon is forced off the quasibreather trajectory by energy conservation.

trices defined by
In other words, the Ω 0 eigenvalues are the roots of the characteristic polynomial with all-real coefficients defined by the matrix with all-real entries H + A. Consider the case A = 0; the eigenvalues Ω 0 are the eigenvalues of H, which is composed of two mirrored copies of the real spectrum of the single-block operator 1 2ω (ω 2 − ∇ 2 + V 0 ). The addition of A only introduces couplings between these two sectors; since it is also antisymmetric, these couplings are equal and opposite in sign. If we start from a spectrum of H with no overlap between its two sectors, the addition of A  FIG. 18b. The power radiated by a simulated sine-Gordon oscillon versus the central fundamental frequency. On this plot, we've indicated the onset of linear instability ω ≈ 0.88 calculated using our eigenvalue code described in appendix C, and the instance of energetic death ω ≈ 0.94 described in appendix A 4. This figure represents the a consequence of linear instability: growing quasiperiodic oscillations. The specific magnitude of this effect depends on initial conditions and environmental perturbations (see figure 11 for an example where oscillations are suppressed). Whether or not the unstable mode can become large enough to destroy the oscillon, the perturbation itself has a radiation component, which may significantly modify the lifetime. In this particular case, the unstable mode's frequency ω ± ReΩ0 approaches the oscillon frequency ω towards the end of life, leading to growing beats (see figure 18a). The loop of death at the end of the evolution occurs because the central oscillon rapidly becomes a mix of first and third harmonic, causing the central frequency to be larger than 1.
will bring the two mirrored 'ground states' together from Ω 0 = ±E ground to the value of Ω 0 = 0. From the perspective of the characteristic polynomial, this corresponds to the two roots becoming degenerate before turning imaginary. In other words, complex eigenvalues must appear by first passing through an inter-block degeneracy. Therefore, the meeting of the two ground states defines the boundary between periodicity (i.e. an all-real spectrum) and instability (i.e. complex spectrum). If the spectrum of H is bounded below by 0, then the meeting of the two states will produce purely imaginary eigenvalues. This result should be compared to [49]. In general, the symmetries of (C6) that led to this result are only approximate, and therefore we should expect the first nearly-stable eigenvalue to be close to zero in general.
However, because the oscillon lives in an open box, we must ensure that (C4) is endowed with transparent boundary conditions. Such radiative boundary conditions depend on the momentum of the outgoing mode (Ω 0 + nω) 2 − 1.
Eigenvalue problems with radiative boundary conditions have been studied in the non-relativistic limit in [83]. Crucially, the calculations FIG. 19. A visualization of how linear instability emerges in the simplified model of appendix C. The boundary of stability is described by eigenvalues meeting at zero. The plot describes the solutions to the eigenvalue equation (C6) in the case of a simple Gaussian background, in which the fundamental oscillon mode is taken to be S1(r) = A exp −r 2 /2σ 2 . The plot background describes stability as a function of the two Gaussian parameters, the oscillon amplitude A and width σ; for oscillons of sufficient width and amplitude, there are eigenvalues Ω0 with negative imaginary part, and thus the oscillon is unstable. We show the eigenvalues nearest to zero for three points in this parameter space: stable (green), borderline unstable (yellow), and unstable (red). The real eigenvalues closest to the origin become degenerate at zero on the boundary of stability; they further split into purely imaginary conjugates in the instability region.
of [83] depend on the existence of a uniformizing variable u(Ω 0 ) in which the outgoing momentum becomes a rational function of u. As far as we are aware, no such uniformization procedure is known for the relativistic case with two channels, or more generally for any case with more than two channels.
Using series approximations and a uniformizing variable u(Ω 0 ), we show in appendix D 3 that it is possible to express the boundary condition as a polynomial for |Ω 0 | < 1/2 for ω > 1/2 in the case of parity or ω > 3/4 without parity. Using standard linearization techniques, we may reduce this polynomial eigenvalue problem to a generalized eigenvalue problem, for which numerical software is plentiful. This is the method applied to the stability analysis in Figure 18a. Outside the red region, the nonlinear oscillations are centered on y = 0, representing that the oscillations stay bounded independent of phase. However, for |α| < 0.5, only oscillations of a particular phase remain bounded, indicating that y = 0 has become hyperbolic (see left inset). Inset in the plot are two examples of the slow oscillation trajectories. For |α| < 0.5, the red stable trajectories have amplitude larger than 0 and are restricted to a finite interval of phase. This generally nonlinear phenomenon represents a special region of stability within the otherwise unstable phase of the Mathieu parameter space. For |α| > 0.5, the red stable oscillations are restricted to a finite amplitude, but are allowed to have any phase. In both cases, large enough amplitude perturbations grow without bound, represented by the black trajectories.

Nonlinear stabilization
In the previous section, we laid out our numerical method for computing the linear stability of the homogeneous perturbation δθ in the background of the deformed physical quasibreather. Modulo technicalities at the boundary, we found that all linear perturbations either are oscillatory, or come in pairs of exponentially growing and decaying modes. The result is that stability, understood to mean that all homogeneous perturbations shrink, cannot be fully explained at the level of linear Floquet analysis.
Thus, stability must originate at higher order in perturbation theory, if it exists at all. In this section, we identify radiation as the mechanism of stabilization accessible to small oscillatory perturbations; specifically, modes which are periodic in the linear stability analysis will couple to radiative modes at higher orders, providing a channel for dissipation. Therefore, we conclude that a sufficient condition for full nonlinear stability is that all modes are oscillatory at the level of linear perturbation theory. Furthermore, we will see that linear instability does not imply nonlinear instability.
We will explore the effect of adding a nonlinear term to a Floquet-type problem by studying the toy example of the Mathieu equation with a quadratic nonlinear term.
To simplify our analysis, we begin by studying potentials with parity, so that the leading order oscillating contribution to V (θ PQB ) is proportional to cos 2ωt. The leading nonlinear term is then proportional to sin ωt. Thus, we will study the nonlinear generalizations of the Mathieu equation of the form 0 =ÿ + (1 + (α + cos 2t + y sin t))y .
In this toy problem, y represents the perturbation δθ to a physical quasibreather whose potential conserves parity. The fact that the linear term is proportional to cos 2t and the quadratic term y 2 is proportional to sin t is a consequence of the symmetry of the potential, which guarantees that polynomials in y of certain parity have the corresponding oscillatory terms. A standard two-timing analysis, along the lines of [85], with 1 demonstrates that the Mathieu instability bifurcation at |α| = 1/2 is unchanged by the nonlinearity around y = 0. However, one difference is the appearance of regions of stability inside the linearly unstable region |α| < 1/2, although large enough y always implies instability, regardless of α [86]. When y is the smallest scale in the problem, we recover the usual Mathieu equation behavior (see figure 20). In summary, linear periodicity is unchanged for small enough y, although linearly unstable modes may become periodic.
Thus, we should expect that the oscillatory modes of the linear equation (C4) remain oscillatory upon introduction of a nonlinear term as long as they are of small enough amplitude. Moreover, the nonlinear terms may convert an otherwise linearly unstable mode into an oscillatory one. Further, the nonlinear interactions of linearly oscillatory modes will necessarily produce radiation, carrying away energy, causing their amplitude to shrink. Thus, sufficiently small linearly oscillatory modes are stabilized by radiation.

Angular perturbations
In appendix D 3, we develop a calculation scheme to solve for the perturbation δθ as a function of t and r. In order to obtain the perturbation equations for δθ, we performed a spherical harmonic decomposition, resulting in a set of decoupled equations for each mode of angular momentum number . These equations differ by the coefficient of the angular momentum effective potential Because this potential is positive, it acts as a repulsive centrifugal term, reducing the perturbation density at the origin. Hence, we expect that perturbations with more angular momentum are typically more linearly stable, since less of the perturbation lies inside the oscillon bulk, although for low angular momentum, the conclusion is case-dependent. An intuition for this comes from applying the stability phases of the standard Mathieu equation (see figure 21).
A similar 1/r 2 term appears in the effective potential for the perturbation upon removing the (d − 1)/r friction term through a change of variables δθ → r (d−1)/2 δθ. This introduces the effective potential This term differs from the angular momentum term in two important ways. First, it can be of either sign: for d = 1, 3 it vanishes, for d = 2, it is repulsive, and for d ≥ 4 it is attractive. Second, it also influences the quasibreather background itself, whereas the angular momentum terms only influence the non-spherical perturbations. Because this potential influences both the background and the perturbation, its effect on stability depends on the specifics of the nonlinear potential.

Perturbative harmonic formulae
Once we have computed the oscillon's nonperturbative modes S n<npert , C n<npert and their orthogonal deformations c n<npert , s n<npert , we may compute the perturbative modes S n≥npert , C n≥npert and their orthogonal deformations c n≥npert , s n≥npert using the procedure outlined in section B 1. Here, we provide explicit formulae for the Sommerfeld matrix S and the block Laplacian operator ∇ 2 . Upon discretizing space such that r = [0, dr, . . . , (N − 1) dr], the matrix S is comprised of all zeros, except for the lower right entry in each diagonal block. To describe this object, we introduce the following notation. The matrix S has four indices: two upper indices labelling the block, and two lower indices labeling the location within that block. In this notation, the entries in the matrix S in (B3) may be written In the same notation, we may describe the Dirichlet-Neumann block Laplacian operator

Potentials with parity
As described in section V, the Fourier basis is a natural basis to describe any scalar potential, since it is not plagued by the same radius-of-convergence issues of, say, the Taylor basis. Here we provide the harmonic factorization of a general scalar potential with parity (45). Taking the first derivative of (45) with respect to θ, we arrive at the following expression for the self-interaction terms in the non-linear wave equation This expression is specific to the case of a potential with 2π periodicity. To accommodate potentials without periodicity, simply replace θ → θ/θ max where θ ∈ [−πθ max , πθ max ] and V m = θ 2 max . In order to keep the following expressions from getting any more unruly, we will present the explicit formulae for 2π-periodic potentials, since the reader may easily convert these expressions to accommodate general periodicity.
By virtue of the Jacobi-Anger expansion [87] we may compute the harmonic expansion of the potential, evaluated as a function of the PQB harmonics where k = (k 1 , k 2 , . . . ). One may write this more compactly in terms of generalized Bessel functions [87]. From this formula, we obtain the expressions for V Sn in (A9) V Sn (S 1 , . . . , ) sin(nωt) .
In general, we may evaluate the full non-perturbative formulae for the mode-potential derivatives V S N and V c N , n(k s n + k c n ) . (D7) Note, the δs in this equation are Kronecker δs, but we use a parenthetical argument to keep the expression readable. From these expressions, we may derive useful formulae for important cases of interest. Here, we present two examples for illustration, and because the reader may find them particularly useful in generating oscillon profiles of their own. First, in the case that the fundamental mode S 1 dominates and all other modes are perturbative, we have the following source term with J s = 0 and the following mass matrices (J 3−3 (mS 1 ) − J 3+3 (mS 1 )) (J 3−5 (mS 1 ) − J 3+5 (mS 1 )) · · · (J 5−3 (mS 1 ) − J 5+3 (mS 1 )) (J 5−5 (mS 1 ) − J 5+3 (mS 1 )) . . . . . .
with − corresponding to V C and + corresponding to V s . The formulae when there are more non-perturbative harmonics follow the same pattern.

Formulae for linear stability analysis
Here we provide the mathematical details to accompany appendix C. We restrict ourselves to ω > 1/2 in the case of parity and ω > 3/4 otherwise. This restriction is to ensure that the following series approximation converges on the disc |Ω 0 | < 1/2. If one is certain that the unstable eigenvalues occur in a smaller disc, the restrictions on ω may be weakened significantly, and indeed, this is often the case.
The outgoing boundary conditions depend on the momentum of the outgoing mode, which is an irrational function of Ω 0 . In order to convert the irrational eigenvalue problem (C4) into a polynomial eigenvalue problem that may be solved with standard techniques, we need to approximate the momentum (Ω 0 ± nω) 2 − 1 by a poly-nomial. One's first intuition might be that the Taylor series of the momentum expanded about Ω 0 = 0 would be a good approximation. This intuition is good for the higher harmonics, since Ω 0 is often much smaller than nω. However, this series only converges inside the disc |Ω 0 | < 1 − ω for n = 1, which is not sufficient to compute the Lyapunov exponent of the linear perturbation δθ. A more sophisticated approximation is necessary in order to capture the behavior of the momentum as a function of Ω 0 on a disc that remains finite size as ω → 1.
We then Taylor expand around x = 0, yielding the following series that converges on the disc of radius 1/2 centered on Ω 0 = 0 for ω > 1/2 where the second equation is just the ordinary Taylor expansion centered on Ω 0 = 0 for n ≥ 2. The factor 2ω(ω − 1 ± Ω 0 ) is not yet a polynomial. We utilize the technique of uniformization [83], where we define the complex variable u such that 2ω(ω − 1 ± Ω 0 ) becomes a polynomial in u, This definition turns (D13) into rational functions of u, allowing us to rephrase (C4) as a polynomial eigenvalue problem. The zero eigenvalues now live at the four roots of the equation u 4 = −1, around which we perform a small eigenvalue search using the Krylov subspace methods implemented in Matlab. Once the u eigenvalues and eigenvectors have been computed, we must convert back to check that they correspond to eigenvalues of the original irrational eigenvalue problem. In short, we have reduced the original irrational eigenvalue to a polynomial eigenvalue problem of degree at-least 4, depending on the degree of accuracy one wants to achieve.
Finally, we define the matrix S which encodes the nonderivative term in the Sommerfeld radiation condition, which can be written as a sum of the matricesŜ λ , which are matrices of all zeros except for the lower right entry of the λth diagonal block, which is 1. This entry corresponds to the outer boundary of the grid, with λ ranging from −L to L, and the upper left block corresponding to λ = −L, where L is the chosen order of the Floquet expansion. In this notation, the non-derivative part of the Sommerfeld boundary conditions may be written where M = min(2j max,1 + 3, 4j max,n + 4), where j max,n is the order of the Taylor expansion of the nth momentum eigenvalue.
Collecting terms in (D13), we have the following expressions for the coefficients of theŜ λ matrices, for λ = 1 and λ = n > 1 respectively Finally, we define the matrix of frequencies where the even entries are dropped when V has parity. These versions of Ω and S are not to be confused with those used to solve for the physical quasibreather (B3), as the correct version to use will always be clear from context. With these definitions, the irrational eigenvalue problem (C4) has been reduced to the polynomial eigenvalue problem where the matrices M i are defined To summarize, the precision of this approximation can be increased to the desired level by 1. increasing the resolution of the radial grid by reducing dr, 2. increasing the physical radius of the simulation r out , 3. increasing the number L of Floquet blocks kept in the expansion, 4. increasing the order j max,λ of the momentum expansions, 5. increasing the number of PQB harmonics kept in the background.

Appendix E: Explicit time evolution -numerical methods
Throughout the text, we refer to explicit numerical simulations for validation of our results. Here we outline the numerical setup used to compute the time evolution of the field θ in the equations of motion (2), and the methodology used to measure oscillon frequency ω and radiated power P .

(E1)
We introduce the variable v = rθ, which eliminates the friction term. We now discretize time and space, with time steps dt and radial steps dr, and introduce the notation v(N dt, M dr) = v N (M ) (E2) Finally, we define the ratio ξ ≡ (dt/ dr) 2 . In this notation, the equations of motion lead to the following leading-order finite difference equation Dirichlet boundary conditions are imposed at the origin by fixing v N (0) = 0. The oscillon is assumed to be evolving in empty space, and therefore the box size must be chosen large enough that the radiation from the oscillon reflected off the outer boundary does not propagate backwards and interfere with the oscillon itself. The length scale of the nth harmonic is 2π/ (nω) 2 − m 2 .
During an instance of destructive interference, typically the fifth harmonic dominates, and in rare cases the seventh may contribute significantly. Since 2π/ √ 7 2 − 1 ≈ 0.9, we choose a safe value of dr = 0.1/m, about 10 times smaller than the length scale of radiation at the highest possible frequency. We find that ξ = 1/4 leads to stable evolution for dr of order 0.1/m. To check that this choice of dr is good, we increased the resolution by a factor of 2 and 4 which resulted in marginal discrepancies.
The frequency of the oscillon is then measured by tracking the times at which v N (1) crosses through zero. The outgoing flux is measured outside the oscillon bulk, typically between 20 and 100 in units of the mass. We do not measure the flux too far from the source, since the different frequency modes travel at different velocities, and the PQB formalism does not account for this dispersion.