Symmetric multifield oscillons

Oscillons are long-lived, spatially localized field configurations, which are supported by attractive non-linearities in the scalar potential. We study oscillons comprised of multiple interacting fields, each having an identical potential with quadratic, quartic and sextic terms. We consider quartic interaction terms of either attractive or repulsive nature. In the two-field case, we construct semi-analytical oscillon profiles for different values of the potential parameters and coupling strength using the two-timing small-amplitude formalism. We use analytical and numerical techniques to explore the basin of attraction of stable oscillon solutions and show that, depending on the initial perturbation size, unstable oscillons can either completely disperse or relax to the closest stable configuration. We generalize our analysis to multifield oscillons and show that the governing equations for their shape and stability can be mapped to the ones arising in the two-field case. Finally, we study the emergence of multicomponent oscillons in one and three spatial dimensions, both numerically and through Floquet theory.


I. INTRODUCTION
Oscillons are part of a large family of soliton-like structures. The stability of these solutions can be traced back to non-linear terms in the theory, since a linear theory will in general force rapid dispersion of any localized wave-packet. In general, solitons can be classified according two characteristics. They are either completely static (ignoring propagation) or oscillate during their lifetime. Furthermore, their stability is either due to some conserved charge or to interactions between nonlinearities and dispersive effects.
Among such non-linear dynamics, oscillons [1-10], localized long-lived objects, have attracted significant attention in the early Universe cosmology community. Oscillons can arise in scalar field models where the potential is quadratic near the minimum and becomes shallower than quadratic ("flattens out") at larger field values. Currently, inflationary models that are preferred by observations are models that contain plateau-type potentials [11][12][13][14], which satisfy the necessary condition needed to support oscillons. Furthermore, numerical simulations have shown that the post-inflationary fragmentation of the inflaton in such models leads to copious production of oscillons [15].
Recently oscillons have seen renewed interest. Numerical simulations have revealed that oscillon formation after inflation is accompanied by the generation of significant amounts of gravitational waves [16][17][18][19][20][21][22][23][24]. Since gravitational waves can be one of the very few tools able to probe the end of inflation, the dynamics of oscillon formation, evolution and stability are becoming an essential part of early Universe cosmology.
Our current understanding of fundamental physics suggests that multiple scalar fields are likely to be present at high energies. In spite of this, most of the work on oscillons has focused on singlefield models, ignoring decay channels of oscillons to other fields and -more interesting-the dynamics that can lead to oscillons comprised of multiple fields. Numerical simulations have uncovered twofield oscillons arising after hybrid inflation [25] and in an SU (2) gauged Higgsed model [26], for which a semi-analytical construction of the observed oscillons was given in Ref. [10], along with a detailed study of their parameter dependent stability. The analysis of oscillons in an Abelian Higgs model can be found in Refs. [27,28]. While a type of composite Q-balls have been found in Ref. [29] and further studied in Ref. [30], their comparison to composite oscillons is beyond the scope of the present work. Finally, oscillons were found in an SU (2) × U (1) model [31,32], inspired by the electroweak sector of the Standard Model. In the case of the SU(2) Higgsed oscillon, a particular mass ratio between the Higgs and W fields was required for the oscillon to be stable. In particular, oscillons where the Higgs mass was almost twice the W mass were found to be stable over very long time-scales. This observation (which was examined in detail in Ref. [10]) can be taken to imply that rather special conditions need to exist for multi-field oscillons to arise.
The current work aims to provide a detailed look into the conditions for the existence and stability of multi-component oscillons, albeit in a simplified model. We organize our presentation as follows. In Section II we describe the two-field model that we consider. In Section III we use the small-amplitude two-timing analysis to construct solutions for composite oscillons, comprised of the two fields oscillating in a phase-locked configuration. Section IV contains the stability analysis of oscillons against long-wavelength perturbations, which is an extension of the Vakhitov-Kolokolov criterion in the case of two interacting fields with identical potential structure, going beyond the specific potential that we chose. In section V we find qualitative and quantitative estimates for the lifetime of stable oscillons that were found in previous sections, following the methods developed in Refs. [33][34][35]. In section VI we generalize our analysis to a model containing an arbitrary number of interchangeable interacting fields. Finally, Section VII explores the emergence of multi-field oscillons from a variety of initial conditions. We provide our conclusions and outlook on future work in Section VIII.

II. MODEL
We consider a model of two identical real scalar fields, each with a quadratic-quartic-sextic potential and a two-to-two interaction term, which can be either attractive or repulsive. This form of the potential was first analyzed for single-field oscillons in Ref. [6], where a large sextic term was introduced to stabilize three-dimensional oscillons in a symmetric quadratic-quartic potential.
We ignore the expansion of the Universe 1 and restrict ourselves to an action consisting of the two scalar fields on a Minkowski background with metric signature (−, +, +, +): Before proceeding, it is worth noting both the restrictions that the action of Eq.
where we included terms up to sextic order in the fields and neglected higher order contributions in α in each term. For that we assume that the parameter α, which is inversely related to the field-space curvature, is small. We see that this expansion has the characteristics that we consider in the action of Eq. (1): two scalar fields with identical potential parameters, a negative quartic term and a very large sextic term (for α ≪ 1), as well as a quartic coupling. Ref. [37] reported the formation of oscillons during two-field preheating after inflation on a T-model potential, but did not provide a detailed analysis of their properties. Of course, the true T-model expansion differs from the idealized action of Eq. (1), in that it contains sextic interaction terms and it provides less freedom for choosing the various potential coupling strengths (see Fig. 1). Furthermore, the analysis of the T-model must also take into account the non-canonical kinetic structure, leading to the derivative couplings single-field systems [8]. In the present work, we instead focus on the symmetric sextic potential of Eq. (1) and consider canonical kinetic terms for the two fields. An extension of the current analysis to realistic α-attractor models, taking into account the intricacies of derivative couplings, can reveal interesting phenomenology for this family of well-motivated and observationally relevant models and is left for future work.

III. TWO TIMING ANALYSIS
In this section we construct approximate profiles of two-field oscillons and compare them to the single-field oscillons found in Ref. [6]. Eq. (1) contains several physical parameters, such as the mass of the two fields m, the coupling Λ and self couplings g and λ. It is convenient to work with dimensionless space-time variables and fields. This is done by the rescalings x µ →x µ = x µ m, For simplicity, in what follows we define {φ 1 ,φ 2 } ≡ {φ,χ} and drop all tildes. The Lagrangian of Eq. (6) leads to a system of two coupled equations of motion for ϕ and χ. In this work we focus our attention on spherically symmetric solutions ϕ(⃗ x, t) → ϕ(r, t), leading to the equations of motion In d spatial dimension the above equations are altered by the substitution 2 The attentive reader will immediately notice that these equations are symmetrical under exchange of the fields, meaning that sending χ ↔ ±ϕ in either of the equations gives back the other. This greatly simplifies the analytical search for oscillons which is more difficult (if at all possible) in general multi-component systems.
It is well established in the literature that oscillons live on long time-and large spatial scales [40][41][42]. This suggests a perturbative approach to extract oscillon solutions from the equations of motion, known as the small-amplitude two-timing analysis. The idea is that oscillons exhibit 2 Additionally, in d dimensions the rescalings of the fields will also be different. behaviour on two time scales, one capturing the natural frequency of the free field theory and one capturing the correction to this frequency characteristic for the non-linear potential. Furthermore, we attempt to capture oscillons which are slowly varying in space, hence are "broad". This behaviour is found by introducing new time and space variables where α is a free parameter and ϵ ≪ 1. While we expect all spatial variation to occur on the scale of ρ, the time variation occurs on two scales: t and τ . The (double) time derivatives in Eqs. (7) must now be interpreted as full time derivatives, leading to the equations Lastly, we assume that the amplitude of the oscillon profile also scales with the expansion parameter ϵ, allowing us to compute the effects of the non-linearities order-by-order. We therefore consider solutions of the form Inserting Eqs. (10) into Eqs. (9), the lowest order equations in ϵ are which are identical to the harmonic oscillator equation and capture the main oscillatory behaviour of the oscillon. The equations in the next non-trivial order in ϵ are Notice that we used the fact that g is large and have therefore written g = 1/ϵ 2 . The solutions to the harmonic oscillator Eqs. (11) are trivial. Therefore we can write ϕ 1 = Re{A(ρ, τ )e −it } and where A(ρ, τ ) and B(ρ, τ ) are complex functions 3 . Inserting this solution 3 Note that this procedure is equivalent to going to the non-relativistic limit, where a real scalar field can be written as the real part of a complex wavefunction Ψ(x, t) exp −imt .
into Eqs. (12) and eliminating secular terms gives us the envelope equations of the system of oscillons These equations are of the Nonlinear Schrödinger type [43]. They control the behaviour of the oscillon on long time-and large spatial scales. To find solutions of these equations that are localized in space we insert the "oscillon ansatz". Since by assumption the oscillon is just some localized structure oscillating in time we should look for solutions of the form, A(ρ, τ ) = a(ρ)e ic 1 τ and B(ρ, τ ) = b(ρ)e ic 2 τ ; where a and b are real functions determining the spatial character of the oscillon. Due to the symmetry of the potential we can also set c 1 = c 2 = c. Inserting this into Eqs.
13, we obtain the profile equations of the oscillons.
where we set c = 1/2. A few remarks about the parameter c are now in order. c is in essence a free parameter, as long as it is not too large. This is because it can always be absorbed into a redefinition of ϵ. We do require it to be positive, since the oscillon should behave as a, b → 0 for ρ → ∞, since they are localized solutions. We can therefore set c = 1/2 without loss of generality.
An interesting property arising due to the symmetry of the system now becomes apparent.
Namely, any localized solution a(ρ) of the equation directly solves the system of Eqs. (14) if we set a = b. We arrive at the same equation if we instead choose a = −b, meaning that the two fields oscillate out of phase. It is worth noting that in the one other case in the literature, where two-field oscillons were constructed using the two-timing analysis, Ref. [10], the masses of the two fields were chosen to have a 2 : 1 ratio, meaning that for half the period of the light field, they were in phase, and for the other half they were out of phase.
The difficulty of finding oscillon solutions in the coupled system of Eqs. (14) is therefore greatly reduced and can be related to the solutions for single-field oscillons in a quartic-sextic potential, which were studied in Ref. [6].
In order to solve Eq. (15), we must first define the boundary conditions. At large distance from the origin, far away from the oscillon core, the non-linear terms become subdominant, hence Eq. (15) can be approximated as ∂ 2 ρ a + (2/ρ)∂ ρ a ≈ αa, leading to a(ρ) ∝ e − √ αρ /ρ. Finally, we require regularity at the origin. Eq. (15) can be re-written as where would be the conserved energy for a one-dimensional system with the same potential. The right hand side of Eq. (16) is non-positive, hence the energy of any solution E ρ will decrease as ρ → ∞.
Since the energy of a localized solution must be zero in the far-distance regime (this can be seen by inspecting Eq. (17)) we conclude that E ρ ≥ 0 at ρ → 0. It is best to think about this in terms of the trajectories of solutions in phase space (a, ∂ ρ a). Since we're interested in functions that are smooth at the origin, all trajectories start on the ∂ ρ a = 0 axis. Any trajectory will continuously intersect curves with E ρ = const, and the constant value of these curves must decrease as ρ → ∞ . The E ρ = 0 curve then defines the boundary between localized and non-localized solutions. Namely, a localized solution must intersect this curve in the origin. The requirement that the oscillon is smooth at the origin, meaning ∂ ρ a| ρ=0 = 0, and the fact that E ρ ≥ 0 there, leads to α < α c ≡ (27/160)(1 + Λ) 2 , generalizing the constraint found in Ref. [6] for single-field oscillons 4 , or equivalently Λ = 0. For α → α c the oscillon becomes infinitely wide with an amplitude a 2 (ρ = 0) → (9/10)(1 + Λ).
It turns out that only a countable set a n (ρ) of localized trajectories can be drawn in phase space [44]. Here n is the amount of nodes of the solution. We conclude that there is exactly one zero-node solution of the profile equation, and thus one zero-mode oscillon for every choice of ϵ and α. We do not pursue solutions with n ≥ 1, since they will have a higher energy than their zero-node counterpart, and thus are expected to be unstable. We used both the shooting and relaxation methods in order to find the corresponding oscillon profiles and check our results. Fig. 2 shows some characteristic profiles for Λ = ±0.5. For the one-dimensional case, Eq. 15 can be solved analytically to give ϕ(r) = 9 10 (1 + Λ) where u 2 = 1 − α/α c . For Λ = 0 we recover the results of Ref. [6]. respectively, meaning that oscillons can acquire larger amplitudes in a system with attractive interactions. Fig. 3 shows the width and height of the oscillons as a function for various interaction strengths between the two fields, attractive or repulsive, rescaled by the small-amplitude parameter ϵ. We see that in the non-interacting case of Λ → 0, we recover the single-field results of Ref. [6]. Furthermore, fields which interact repulsively lead to a smaller range of oscillon amplitudes, ending in no oscillons for Λ → −1, as is also evident by the form of α c ∼ (1 + Λ) 2 . On the contrary, attractive interactions allow for oscillons with a larger amplitude.  narrower oscillons exist. In particular, the minimum width for oscillons with Λ = −0.5 is almost triple the corresponding value for Λ = 0.5.

IV. STABILITY ANALYSIS
It has been shown repeatedly in the literature that no true breather solutions like the oscillon can exist in nonlinear systems (with the notable exception of breathers in the one-dimensional Sine-Gordon model). In general there will always be (classical) outgoing radiation in the tails of the oscillon. Classically, this radiation is exponentially suppressed, which is why oscillons can be extremely long-lived [45]. Quantum mechanical radiation might play a more important role in real physical systems [46], but is beyond the scope of our present work. Radiation will in general perturb the oscillon system and it is therefore necessary to assess the stability of oscillons with respect to small perturbations. In the model we are investigating, let us consider small fluctuations δ(x, t), ∆(x, t) ≪ O(ϵ) added to the oscillon solutions as ϕ(x, t) = ϕ osc (x, t) + δ(x, t) ,  Plugging these into Eqs. (7), the equations of motion for the two fields, and keeping only terms linear in δ and ∆ leads to 5 Since the exact shape of the initial perturbation is in essence unpredictable, a full linear stability analysis would require solving Eq. (20) for arbitrary initial conditions. This requires solving the full Floquet matrix for the coupled nonlinear system. It turns out that, for this system, there is a useful simplification if we only consider perturbations that are about the same size as the oscillon itself. This is in essence an extension of the Vakhitov-Kolokolov criterion that has been used to assess stability of single-component oscillons [6]. 5 The terms that were ignored here (the terms only involving ϕosc and χosc) can by nature of the two-timing analysis not source instabilities. However, in Section V we will show that these terms are a source of outgoing radiation.
In Section III we have already shown that this system supports oscillons of the form where the frequency 6 is ω = 1 − αϵ 2 /2. The symmetry of the system dictates that we have We present the details of the derivation in Appendix A. Fig. 5 shows the regions of stability, based on Eq. (22).

B. Oscillon Dynamics and Decay
Having constructed two-field oscillon solutions and extended the V-K criterion to examine their stability, we move to numerically demonstrate our results and understand the relevant time-scales of interest. As a benchmark case of unstable oscillons, we use profiles corresponding to α = 0.04α c and ϵ = 0.08, for both attractive (Λ = 0.5) and repulsive (Λ = −0.2) interactions. which we perturb using the most unstable mode-functions computed in Section IV. However, while the perturbation shape is given, the perturbation amplitude is a free parameter. We define the perturbation amplitude as the relative size of the oscillon and perturbation mode-function at the origin and consider various values for δ magn. . We investigate both the effect of placing the perturbation on both fields (with equal magnitude), and the effect of placing the perturbation on only one of the two fields and initializing the other field in the "perfect" oscillon configuration. 6 In several oscillon studies, the frequency is taken to be of the form ω = 1 − c · h 2 , where h is the oscillon height and c is some O(1) constant. This matches our expression to lowest order in perturbation theory when h ∝ ϵ ≪ 1, which is true in the context of the small amplitude analysis used to construct oscillon solutions.  wave-packet simply disperses. We must note here, that initializing the two fields out of phase a = −b and using the corresponding perturbation derived using Eqs. (A12) gives indistinguishable curves, and thus we do not plot them for clarity. This however confirms our analysis, that the phase of the two fields in the initial oscillon configuration, 0 or π, does not affect the subsequent Figure 6. Upper Panels: The evolution of the oscillon peak energy density (left) and width (right) as a function of time t(m −1 ) for Λ = 0.5, α = 0.04α c and ϵ = 0.08. We perturb the oscillon by the most unstable mode with amplitude δ magn. = −0.5%, 0.05%, 0.5% (green, magenta and red respectively). The brown curve corresponds to perturbing only one of the two fields with δ magn. = 0.2% and initializing the other field as the "perfect" oscillon profile. The data for both fields is plotted for this curve (dashed and dotted line) but the difference is not perceivable since the fields quickly attract and undergo simultaneous collapse.
Lower Panels: The same simulation for Λ = −0.2 with again α = 0.04α c . Here the curves correspond to δ magn. = −0.1%, 0.01%, 0.1% (green, magenta and red respectively). We again plotted brown curves in a scenario where we only perturbed one of the fields with δ magn. = 0.1%. The fate of the two oscillons differs as one of the field eventually disperses (dotted) while the other settles into a stable single-field oscillon (full).
It is clear however that all of these multi-field oscillons contain unstable modes.
dynamics. We will show numerically in Section VII that the phase relation between the components of the oscillon arises dynamically when they form, leading to phase-locked fields, oscillating either in or out of phase.
Furthermore, one can ask whether having the exact form of the most unstable modes computed in the previous section is a necessary condition for oscillon decay. In single-field oscillons, one can see that any perturbation can be decomposed in a basis of eigen-functions of the operator H 1 H 2 [see Eqs. (A9) and (A10)], meaning a generic perturbation will include a component along the unstable mode. In two-field systems, it is interesting to see how an initial perturbation is transferred from one field to the other. To check the validity of our analysis, we perturbed only the ϕ field, leaving the χ field profile identical to the perfect oscillon solution. Perturbing one of the fields with δ magn.
while leaving the other field intact leads in some cases to a similar evolution as perturbing both fields with δ magn. /2. For the specific cases that we plotted, this occurs for Λ = 0.5 when perturbing one of the fields with δ magn. = 0.2%, leading to almost indistinguishable evolution from perturbing both fields with δ magn. = 0.1%. We do not plot both these simulations in Fig. 6, because the results are visually identical (instead only one of these cases is plotted). This indicates that due to the interaction of the two fields and the symmetry of the system, a small initial perturbation is quickly distributed evenly among the two fields. In other cases we observed the complete dispersion of one field and the relaxation of the other towards a stable single-field oscillon solution (see lower panels of Fig. 6). In particular this occurred for Λ = −0.2, when we perturbed only one of the two fields by δ magn. = 0.1%. From our numerical investigation it remains clear that the symmetric multi-field oscillons failing the VK-criterion contain unstable modes and will in general decay. However, when perturbed outside the ϕ = ±χ Ansatz, the dynamics can either restore or destroy the equality between the two fields. We will see that for Λ > 0 the two-field configuration is a strong attractor of the system. We note that it is not necessary to perturb the oscillon with the exact form of the perturbation given by the stability analysis of this Section, since unstable oscillons will eventually decay if perturbed, albeit on slightly different time-scales and with different final state, depending on the details of the perturbation. Finally, we also perturbed the oscillon with negative δ magn.
(meaning that the unstable mode is "flipped"). We observed that this leads to a collapse instability more often than the positive δ magn. case. This highlights the sensitivity of the fate of the unstable oscillons on the initial conditions. Finally, we perturbed oscillons that are characterized as stable by the two-field version of the V-K criterion in various ways, seeing no decay over the total run of the simulation, lasting several hundred units of "slow time" τ = ϵ 2 t.
It is now important to understand this behaviour. Fig. 7 shows the trajectories of several perturbed oscillons on the height-width curve for oscillons 7 for Λ = 0.5 and Λ = −0.2. We see that unstable oscillons will surely move away from their initial configuration if perturbed. However their ultimate fate is not unique. It is known that localized solutions of the three-dimensional non-linear Schrödinger equation are unstable and undergo a collapse instability, meaning that their width becomes smaller and their amplitude becomes larger. For oscillons that are "mildly perturbed", 7 In this plot, the height and width curve is derived from the energy density, not the profile of the oscillons. It is thus related to the one shown in Fig. 3, but is not identical to it. Eq. (13). Once the oscillon amplitude grows enough, so that the |A| 4 A term becomes important, the oscillon is stabilized, which occurs near the minimum of the height-width curve. Beyond this point, the oscillon solution satisfies the V-K criterion. Because our simulation box does not allow for the energy to escape (it does not possess absorbing boundary conditions), the oscillon will never truly reach the "ideal" shape, and will instead stay in the vicinity of the first stable point on the height-width curve.
Oscillons that are perturbed by a large amount, which puts them too far away from the heightwidth curve, do not undergo this slow collapse instability, which would take them to a stable configuration. Using different perturbation types on the initial oscillon configuration amounts to examining the basis of attraction for two-field oscillons, finding the approximate size of the initial perturbation below which the oscillon is allowed to "rearrange" itself into a stable configuration and above which the oscillon is completely destabilized. An important question remains, as to the oscillon configurations that arise naturally after inflation: it is not known if a universe that is governed by the Lagrangian of Eq. (1) during preheating will be dominated by single or multi-field oscillons. We are currently working on the corresponding three dimensional lattice simulations required to address this point and will present the results in a subsequent publication.

V. OSCILLON LIFETIME
It is a well known fact that oscillons are not exactly stable, but eventually lose their energy through an exponentially suppressed radiative tail [47,48]. To make an assessment of the (cosmological) importance of the stable oscillons that were found in the previous sections we wish to find an estimate of their lifetimes. To do this we follow the analyses performed in Refs. [33][34][35].
The core reason why oscillon solutions are in general not stable (classically), even though they are long-lived, is that plugging in an oscillating solution ∼ Φ 0 cos(ωt) into an equation of motion with nonlinear terms will in general lead to terms that are proportional to higher harmonics ∝ cos(jωt) with j > 1. In general, these terms can be made to cancel out using some perturbative procedure (e.g. the two-timing analysis); but since in practice these perturbative techniques must end at some finite order and do not converge, there will always be residual terms of higher harmonics in the equations of motion. These terms will act as a source for radiative modes of the oscillon solution. 8 In order to estimate the decay rate of the oscillons constructed in Section III we follow the methods of Refs. [33][34][35]. To study the relevant radiative modes we perturb our oscillon solutions ϕ = ϕ osc + δ, χ = χ osc + ξ, plug these into the equations of motion and linearize to obtain In the above equations we use the notation V ,ϕ (ϕ osc , χ osc ) ≡ ∂ ϕ V (ϕ, χ)| ϕ→ϕosc,χ→χosc for derivatives of the potential with respect to the fields. We now plug in our oscillon solutions ϕ osc = Φ 0 cos(ωt) and χ osc ∼ X 0 cos(ωt), where ω = 1 − αϵ 2 /2. By virtue of the two-timing analysis, all terms proportional to cos(ωt) in the first square bracket of Eqs. (24) cancel and we obtain Focusing on the solutions where the two fields within the oscillon oscillate in phase Φ 0 = X 0 (the case where Φ 0 = −X 0 is analogous); we add Eqs. (25) and introduce the new radiative mode Ψ = δ + ξ (in the case where the oscillons oscillate out of phase we would introduce Ψ = δ − ξ) to obtain the single equation We are thus interested in wave-like solutions Ψ rad (r, t) of Eq. (26). Using the symmetry of our system we can then deduce that half the amplitude of this wave comes from the oscillon in the ϕ-field and the other half from the oscillon in the χ-field. Notice however that this calculation is equivalent to calculating the radiative mode of an oscillon in a single field with a slightly altered quartic term in the potential As suggested in Ref. [33], when calculating the outgoing radiation of an oscillon in this type of potential, ignoring the spacetime-dependent effective mass term in the equation of motion for the perturbative mode will not alter the final result significantly. Thus we ignore the time-dependent coupling of Ψ on the left-hand-side of Eq. (26) and obtain the simpler equation As was shown in Ref. [33], this equation can be solved by expanding Ψ(r, t) = ∞ j=1 Ψ j (r) cos(jωt).
Obviously only Ψ 3 (r) and Ψ 5 (r) will not vanish in this specific case. The rate of energy loss of the oscillon is dE/dt = 4πr 2 T 0r , where T 0r is the Poynting vector in our spherically symmetric set-up.
By averaging over time in the far distance regime radii where k j = (jω) 2 − 1 and S 3 (k) and S 5 (k) are the Fourier transforms of the terms on the right hand side of Eq. (28) proportional to the third and fifth harmonic respectively. Thus and We find, similarly to Ref. [9], that S 3 ≫ S 5 . The analytic expression for the radiative modes of the oscillon given by Eq. (30) leads to the following expression for the decay rate of the oscillon By using the oscillon amplitude Φ 0 ≡ ϵa(ϵr), which corresponds to the solution of Eq. (15), we computed the value of Γ osc using the method described above. We confirmed that the multi-peak structure of Γ found in Ref. [9] persists for different values of Λ, either positive or negative. In general, the oscillon lifetime increases as one decreases the small parameter ϵ. In our construction, this corresponds to increasing the strength of the sextic term. Fig. 8 shows the results derived by numerically integrating Eq. (33) for three-dimensional oscillons, using the value ϵ = 0.1. For Λ = 0 we recover the results of Ref. [5]. We see that the decay rate increases for Λ > 0 compared to the single field case (Λ = 0) and correspondingly it decreases for Λ < 0. The maximum value of Γ increases from Λ = 0 to Λ = 0.2 and Λ = 0.2 to Λ = 0.5 by around one order of magnitude.

VI. GENERALIZATION TO MANY FIELDS
In this Section we generalize our model to a system of N interacting fields. Following the two-field analysis, we choose the following action where the fields exhibit an "exchange" symmetry: We briefly repeat the derivations presented above for two-field oscillons, showing that that the profile equations are solved by zero-modes of Eq. (15) with an altered cubic term. Furthermore, we show how to extend the V-K criterion to apply to the N -field system.

A. Two-timing analysis
After rescaling the fields and spatial-temporal dimensions in the same way as was done in section II, the action of Eq. (34) results in N equations of motion for the N fields. The equation for the i-th field (assuming spherical symmetry) is Introducing the by now familiar two-timing analysis variables ρ = ϵr and τ = αϵ 2 t, the equations become Finally looking for solutions ϕ i = ϵϕ 1,i + ϵ 2 ϕ 2,i + ... we can write down the O(ϵ) and O(ϵ 3 ) equations where we've again used g = 1/ϵ 2 . The rest of the derivation is entirely analogous to the one in Section III so we won't repeat it here. The procedure results in N profile equations We thus need to find a set of N functions a i (ρ) that solve this system of equations. As before, the symmetry in the model greatly simplifies this question. We need to find localized solutions of the equation Then, setting any n < N number of a i (ρ) = a(ρ) and the remaining a j (ρ) = −a(ρ) solves the system in Eq. (39). Notice that Eq. (40) has the same form as Eq. (15) with an altered cubic term.
It will therefore have localized solutions, which are related to the ones for the two-field model by substituting the coupling strength Λ in the two-field case with the "effective" coupling strength Λ(N − 1) in the multi-field case. This is a proof of existence for such a solution, not uniqueness.

B. Stability analysis
For brevity, we move the full derivation of the stability of these N -field oscillons to long wavelength perturbations to Appendix A. Note however that performing the two-timing analysis and using the ansatz ϕ 0,i = Re A i (x, t)e −it is equivalent to working in the nonrelativistic limit. In this limit, the Lagrangian has a conserved charge proportional to which reduces to the following after inserting the symmetric oscillon solution Since the oscillon has a conserved charge, we can apply the VK-criterion directly. To summarize, the oscillon is stable when ∂Q ∂αϵ 2 > 0. This is just an alternate and equivalent form of the VKcriterion that was given earlier in Eq. (22). The criterion for stability therefore still holds after a generalization of our model to N fields.
To conclude, we constructed oscillons in a model with N fields. Since we work in the nonrelativistic limit it is simple to show that the VK criterion still applies. Although we have not performed a detailed analysis of the lifetime of these oscillons, the result is expected to be analogous to the two-field case. Namely, the computation of the radiation tail of these oscillons can be mapped to an equivalent calculation for single-field oscillons with an altered cubic term in the potential, as done in Section V for the two-field case.

VII. OSCILLON EMERGENCE AND DYNAMICS
Despite the mathematical construction and stability analysis of two-field oscillons, the question of their emergence in realistic scenarios and hence their cosmological consequences remains. In this section we provide a combination of simulations in one and three spatial dimensions, in both an expanding and a static background. The full three-dimensional simulations of oscillon formation after inflation and their GW signatures will be presented in a subsequent publication.

A. Dynamical emergence in 1+1D
To investigate the formation of oscillons in this model we performed a series of one-dimensional simulations of preheating-like scenarios on an expanding background. One of the fields, ϕ, was initialized as a homogeneous condensate, akin to the inflaton at the end of inflation. The other fields of the model χ, θ, ... were then initiated in their vacuum, akin to spectator fields at the end of inflation. The simulations were strongly inspired by those performed in Ref. [49]. We see a clear sequence of events, leading from an oscillating ϕ condensate to a collection of composite oscillons.
• Parametric self-resonance leads to a fragmentation of the inflaton condensate and the formation of single-field oscillons.
• Parametric resonance of the spectator field(s) by the single-field oscillons leads to exponential growth of the spectator, akin to preheating and Floquet theory.
• Once the spectator field acquires a large enough value, non-linear interactions between the two (or more) fields lead to phase-locked long-lived configurations: composite oscillons. positioned at the same locations as ϕ field oscillons. Interestingly, not all ϕ oscillons lead to an efficient amplification of the χ field. Furthermore, we see that the energy density spikes of the two fields are not equal in size. In fact, two-field composite oscillons in one spatial dimension exhibit a persistent energy exchange between the two fields, which we will describe more later.  χ field to be excited in regions of large ϕ field, leading to a more efficient production of large field two-component oscillons. We will return to this point in Sections VII B, VII C and VII D. Finally, we checked the relative phase of the two fields in the two-field oscillons and found no preference for either in-phase or out-of-phase configurations. is nearly the same and the total number of multi-component oscillons is nearly equal to the number of two-component oscillons that emerged in the two-field simulation.

B. Single to multi-field oscillons in 1+1D
The results of Section VII A showcased the two-stage formation of composite oscillons on an expanding background for the one-dimensional case. In order to better understand the transition from single to multi-field oscillons, we focus on simulating single-field oscillons on a static background, which allows us to fully control the initial conditions and compare the numerical results to theoretical predictions. The decay into a multi-field configuration can be understood in the context of Floquet theory. Although obtaining analytical results can be challenging in a non-homogeneous background such as an oscillon, we show in Section VII D that significant intuition is gained by working in the homogeneous limit.
We initialize the ϕ field as an oscillon that is stable in the single field theory (with Λ = 0) and draw the fluctuations of the χ field from a Gaussian distribution with ⟨χ 2 k ⟩ = 1/2k, as in our lattice simulations. Alternatively, initializing the secondary field as a small perturbation of the same shape  as the oscillon does not change the outcome of the simulation (although it somewhat accelerates the decay into a multi-field configuration). We will demonstrate this last point later on in this section. Fig. 12 shows the evolution of the system for Λ = 0.5 and g = 1/ϵ 2 = 5: from a single field oscillon with α = 0.99α c to a phase-locked two-field configuration. The upper panels show the evolution of the width and height (in energy density space) of the two fields. We clearly see the growth of a localized χ configuration, which backreacts on the ϕ-oscillon, leading to a composite oscillon. We see that after t ≃ 1500(m −1 ) the amplitude and width of the two fields oscillate around an average value, which is common for both. Inspecting the data we can postulate the existence of "breather" modes: stable perturbations atop the perfect oscillon solution that oscillate between the two fields (similarly to how we construct unstable modes in App. A). It seems that fewer of such breather modes exist in higher dimensions (see Section VII C).
The oscillation frequency behaves similarly. Initially, since the χ field is comprised of random fluctuations, one cannot speak of a single "coherent" frequency. As instabilities build up over time however, the two fields "lock" up into a single shared frequency as can be seen in the bottom panels. The localized energy in each field also exhibits strong oscillations of up to 40% around the average value. While the two components of the oscillon exchange energy with each other, the total configuration is still stable. This is strongly reminiscent of the behavior exhibited by oscillons in the SU (2) gauged Higgs model [10], where stable two-field oscillons discovered in the reduced spherical Ansatz exhibit similar exchange of energy with a period much larger than the natural period of oscillation of the fields themselves.
Before concluding this section, we would like to explore the dependence of the evolution on the initial conditions for χ, by altering them in two ways: by drawing a different sample from the same aforementioned Gaussian distribution on the one hand and by initializing χ as a small "bump" with similar width to the ϕ oscillon, but very small amplitude.
We see in Fig. 13 that drawing a different sample from the same Gaussian distribution leads to a similar evolution with a phase-locked final state, albeit with the fields oscillating out of phase (in contrast to the situation of Fig. 12). This verifies our assumption in the analytical construction of the oscillon, where we distinguished between in and out of phase oscillons, while at the same time showing that these two phase configurations are the attractors of the system, depending only on the initial conditions.
On the other hand, by closely looking at the initial (linear) evolution of the χ field we see two stages: a small bump arising from the noise and the bump growing exponentially. This can be seen clearly in Fig. 14; where we note similar behaviour for the 1D and 3D systems. By initializing the χ field as a bump, we circumvent the first stage, without losing any information about the dynamical evolution of the two-field system. In what follows, this is exactly what we did in order to increase numerical efficiency.
Finally, we performed a numerical simulation for three fields. In Fig. 15 we show the evolution of a system initialized as a single field oscillon in the primary field ϕ, interacting with two fields initialized as small perturbations. We again see that at late times the three fields oscillate around Figure 13. Time slices of the fields ϕ (full) and χ (dotted) at different times t = 1900, 6000, 8000(m −1 ) (blue, red and green respectively). The parameters are the same as those in Fig. 12, but with a different Gaussian initialization for the field χ. Contrary to Fig. 12 the fields condense into an out-of-phase configuration. Figure 14. The evolution of the energy density (time moves from bottom to top) of the secondary field starting from Gaussian initial conditions. It is clear that, although the initial configuration is very noisy, the field aligns into a smooth configuration as the relevant modes get amplified due to Floquet instabilites.
Left: 1-dimension with α init = 0.99α c , ϵ = 0.45. Right: 3-dimensions with α init = 0.5α c and ϵ = 0.2. In both cases Λ = 0.5 a common value for the width and height, forming a phase-locked three-field oscillon with internal breathing modes between the fields. In the particular example plotted we observed a phase-locked state after t ∼ 5000(m −1 ) where initially one of the fields is oscillating out of phase with the other two fields. However we should report that due to the presence of the breathing modes a spontaneous change in behaviour is also possible, where a different type of phase-locking is suddenly reached. Figure 15. The decay of a single field oscillon to a three field composite oscillon. Although the presence of breathing modes doesn't allow the system to settle into a perfect symmetric state, all fields seem to oscillate around a common configuration. Top: The evolution of the height and width of the primary (blue) ϕ oscillon and the two secondary oscillons (green and red) that are formed after being initialized as a small perturbation. Bottom: the field configurations (both in field space and energy density) at t = 10000(m −1 ) of the primary field (blue) and the two secondary fields (green and red). Here we used α = 0.99α c as initialization of the primary field and ϵ = 0.2.

C. Single to multi-field oscillons in 3+1D
We now move to the the analogous three-dimensional case, anticipating full lattice simulation of preheating in three dimensions. Deferring the full lattice simulation for a future publication, we instead start with a single field spherically symmetric oscillon and study the transition to the multi-field configuration, similarly to Section VII B.
Specifically we initialize the ϕ field in the single oscillon configuration with ϵ = 0.1, 0.2 and several values of α and the χ fields as a small localized overdensity (bump), as described in Section VII B. In Fig. 16 we show the evolution of the width and height of the single field oscillon towards the two-field solution at the end of the simulation. Clearly, similarly to the 1-dimensional case, the oscillon moves towards a composite structure. As can be seen from the error bars, we see that the resulting configuration is much closer to the analytically constructed form of a twofield oscillon than its one-dimensional counterpart. In other words, the internal breathing mode exists, but is suppressed. We believe that the larger size of the internal breathing mode for one spatial dimension is connected to the increased stability for one-dimensional oscillons, compared to their three-dimensional counterparts. Otherwise, the results are qualitatively equivalent to the 1-dimensional case as represented in Fig. 12. It is noteworthy that we need to initialize the oscillon with a value of α much closer to α c in order to see decay to a two-field solution in 1D: α i ∼ 0.95α c in 1D and α i ∼ 0.45α c in 3D; where the subscript i indicates the initial value of α of the background oscillon. This fact can be understood heuristically from the point of view of Floquet theory.

D. Analysis of results in terms of Floquet theory
In this subsection we explore the emergence of the multi-field oscillons as observed in the previous sections using semi-analytic techniques. We do this explicitly for two fields, but the intuition naturally extends to more fields. To first order in perturbations, the equation of motion of the secondary field χ isχ Since χ is a small perturbation at early times, this equation governs its dynamics. Also, to first order in ϵ we write ϕ = ϵΦ 0 (x) cos(ωt). Switching to Fourier space and applying the convolution theorem we obtainχ where Ψ(k) is the Fourier transform of [Φ 0 (x)] 2 . This equation can in principle be solved using Floquet theory, where the solutions are of the form χ k (t) = e µ k t P (t), where P (t) is a periodic function. However the inhomogeneity of the background oscillon couples all different k-modes through the convolution term which makes the challenge of finding the Floquet exponents impossible. To understand the basic dynamics, we investigate the constant amplitude limit Φ 0 (x) = Φ 0 (0), and incorporate the effect of the spatial extent of the oscillon heuristically. In this limit the wavenumbers decouple and it is a simple exercise to find the Floquet exponents µ k . In appendix B we show how to find an approximate solution for the functions χ k (t). They are of the form χ k (t) ∝ Re{cos(ωt)e iµ k τ } where τ = ϵ 2 t is slow time and where Φ 0 is the rescaled height of the background oscillon at the origin and ω = 1 − αϵ 2 /2 its frequency. If Eq. (45) has imaginary solutions the corresponding mode will grow exponentially.
Eq. (45) makes it explicit that a two-field configuration will not emerge for repulsive interactions (negative Λ), so in what follows we focus on attractive interactions. Note that if we find an unstable mode for some ϵ 1 (g = 1/ϵ 2 is the parameter multiplying the sextic term in our Lagrangian) at a wavenumber k 1 , a different choice of ϵ 2 will have an unstable mode with the same Floquet exponent at k 2 = ϵ 2 ϵ 1 k 1 . Since the oscillon width scales like R osc ∝ 1/ϵ, the range of validity of the homogeneous approximation also changes with ϵ. This explains why the decay of a single field oscillon into a two-field configuration does not seem to depend on the exact value of the small parameter 9 ϵ (see e.g. Fig. 16); though it might depend α and Φ 0 (these quantities are linked via the oscillon profile equation; Eq. (14)). This line of argument is also valid in the full inhomogeneous Floquet analysis, although it is impossible to find exact analytical expressions for the Floquet exponents in this case. The reasoning is confirmed in Fig. 17 where we plot the amplification of the energy in the χ field at early times for some of the simulations we performed in three dimensions. There is a dependence on α which is allowed by Eq. (B7). Here we used Λ = 0.5 As expected, the amplification of the field with respect to slow time τ is independent of ϵ. We can use Eq. (B7) to find a crude estimate of the amplification of the localized energy which goes like happens for α ≥ 0.45α c . This is merely a reflection of the invalidity of the approximation we have taken here. Since the oscillon has a width of order R osc ∼ 1/ϵ, this approximation breaks down for some k ≪ ϵ. Correspondingly, when for a given α the oscillon solution is very wide, the approximation is valid for a larger range of k−modes. This helps us understand the difference between the 1D and 3D simulations. In Fig. 19 we plot the dependence of the width and height of the oscillon solution on α. Clearly, the three-dimensional solution has a larger width than the one dimensional solution for α ∼ 0.45α c (also see Fig. 3). In Fig. 18 this information is conveyed by the black line going through the instability band: the approximation is only valid for wavenumbers that are to the right of this line. Note that the line drawn here is merely schematic.

VIII. SUMMARY AND DISCUSSION
Despite the oscillons' ubiquity in non-linear scalar field theories and their assumed presence in the early Universe, very few multi-field oscillons have been found and studied. In the present work we explored a symmetric system comprised of two interacting scalar fields and showed how genuine two-field oscillon solutions can be constructed.
We found two-field oscillon solutions in potentials with either an attractive or a repulsive interaction term. They are qualitatively similar to "flat-top" oscillons found in Ref. [6], with quantitative differences that depend on the sign and strength of the interaction term. The oscillons that emerge in the presence of an attractive non-linearity can be both taller (reaching higher central field values) and also narrower (having smaller width) than their single-field counterparts. In the case of repul- sively interacting fields, the range of possible oscillon amplitudes shrinks, until the repulsive term is strong enough to completely forbid the existence of oscillons, at least within the small amplitude two-timing framework.
The stability of non-linearity-supported localized structures can be assessed through the Vakhitov-Kolokolov (V-K) stability criterion. We formally extended the VK criterion, in order to be used in two-field systems. We checked the VK criterion against numerical simulations (assuming spherical symmetry throughout the evolution), finding excellent agreement between semi-analytical and fully numerical results. Our current proof holds for symmetric potentials and the generalization to arbitrary multi-field oscillons is left for future work. Furthermore, we explored the basin of attraction of stable oscillons, by perturbing unstable initial configurations. We found that, depending on the size of the initial perturbation, the unstable oscillons can either completely disperse, or relax to a stable oscillon configuration with a smaller width and larger height.
Since oscillons are long-lived configurations albeit with a finite lifetime (even classically), we adapted the methods used in the literature to compute the emitted scalar radiation and thus estimate the oscillon lifetime. The results we found are qualitatively and quantitatively similar to the single-field case, showing that the longevity of multi-field oscillons can be -at least-comparable to their single-field counterparts. We showed that the results for the two-field system can be naturally extended to models with an arbitrary number of fields, where each field has the same quadratic -quartic -sextic potential and all fields are pairwise coupled. We were able to construct multi-field oscillons, under certain conditions for the couplings, and prove their stability, by drawing formal analogies to the single-field case. This opens the way for arbitrarily dense oscillons, comprised of multiple fields.
Finally, we studied the emergence of these multi-field oscillons. We performed simulations on an expanding lattice in one spacial dimensions, in a preheating-like scenario. We showed that the fragmentation of the inflaton proceeds in three distinct stages: the inflaton fragments into single field oscillons; the single-field oscillons amplify fluctuations in the spectator field(s); composite oscillons emerged. Prompted by these observations we analysed the decay of single field oscillons into multi-field oscillons, both in one and three spatial dimensions, in Minkowski space. The decay can be understood semi-quantitatively using floquet theory. Since the emergence of multi-field oscillons in three spatial dimensions follows the one-dimensional behaviour in Minkowski space, we expect that multi-field oscillons will also form in full three-dimensional simulations on an expanding lattice, which is left for a future publication.
Overall, our current work provides analytical tools for studying multi-component oscillons and opens up several avenues for future work, such as relaxing the symmetry structure of the potential and the assumption of spherical symmetry. We present the derivation of the V-K stability criterion, applied first to the case of two-field oscillons and afterwards generalized to N fields.

Two fields
For clarity, we will separately analyze the two cases of oscillons where the two fields oscillate either in phase (a = b) or out of phase (a = −b).
a. Two-field system oscillating in phase For oscillons comprised of the two fields oscillating in phase, a = b, the system of Eqs. (20) becomes While Eqs. (A1) is rather complicated, the symmetries of the initial Lagrangian allow us to simplify the stability analysis, by introducing the new variable ξ(x, t) = δ(x, t) + ∆(x, t), which reduces Eqs. (A1) to This equation only depends on ξ(x, t), greatly simplifying our calculations. Since we are interested in perturbations that are about the same size as the oscillon itself we perform the same change of variable as before r → ρ = ϵr. Furthermore, we expect perturbations to oscillate near the oscillon frequency. To capture the growth or decay of the perturbation we should also introduce a new "slow" time variable τ = ϵ 2 t. We also introduce a time-scale related to the corrected oscillation frequency of the oscillon T = ωt. As before we expect behaviour on both time scales, so ξ(x, t) → ξ(ρ, T, τ ), and partial time-derivatives become full derivatives. Looking for solutions that can be written as a power series in ϵ ≪ 1 allows for a perturbative analysis of Eq. A2. The zeroth and first order equations are respectively and The general solution of Eq. (A4) has the form where the functions u(ρ, τ ) and v(ρ, τ ) capture the potential growth of the perturbation on timescales of order τ . The V-K criterion provides a simple way to distinguish cases where perturbation grow exponentially based on the form of the oscillon itself. Namely, inserting Eq. (A6) into Eq. (A5) and eliminating secular terms on the right-hand side leads to equations for u and v where H 1 and H 2 are Hermitian, linear operators defined as Separating variables as u(ρ, τ ) → u(ρ)e Ωτ and v(ρ, τ ) → v(ρ)e Ωτ the problem reduces to the linear equation The question of whether the oscillon is stable to general long-wavelength perturbations is thus reduced to an eigenvalue problem. If the operator −H 1 H 2 has at least one positive eigenvalue Ω 2 > 0, perturbations can grow and the oscillon will in general be unstable. If not, perturbations simply oscillate, leading to an oscillon that is stable, within the limits of the perturbative expansion used to construct it. The problem can be solved using a similar procedure as Vakhitov and Kolokolov [50]. We do not present the entire proof here, as its intricacies will add little to the main goal of analyzing two-field oscillons. The criterion for stability states that max(Ω 2 ) < 0 if and only if dN/dα > 0, where N is defined in Eq. (22).
Before proceeding, we must make a final remark about the validity of this derivation. The criterion we presented above only states whether perturbation of the form ξ(x, t) = δ(x, t) + ∆(x, t) will grow. In principle, the growth of the perturbation might only be present in δ or ∆. However, since the system is completely symmetric under exchange of the fields, we can only conclude that, if ξ(x, t) grows, both δ(x, t) and ∆(x, t) will grow exponentially. The criterion should therefore be valid for the full two-field system and can be used as an indicator for the stability of two-field oscillons, at least ones that are comprised of interchangeable fields.
We see that the equation of motion for ψ in the case of a = −b is identical to the one for ξ for a = b.
The rest of the derivation follows exactly the same steps as before so we won't repeat them here.
Hence the stability of the oscillon will be independent of the phase (0 or π) between the two fields, resulting in exactly the same V-K criterion as before. Namely, there are unstable perturbations of the oscillon if and only if dN/dα > 0, where N is defined in Eq. (22).
Which is just the limit of Eq. (44) in which the oscillon is approximated as a homogeneous background. We can find perturbative solutions by introducing two time scales: the period T = ωt = 1 − αϵ 2 /2t and the slow time τ = ϵ 2 t. We also assume that the solution has characteristic behaviour on both timescales, similar to the two-timing analysis, meaning χ k (t) = χ k (T, τ ).
Finally, we are interested in modes for which k = κϵ where κ = O(1). Eq. (B1) becomes where a dot denotes derivation with respect to T and a prime denotes derivation with respect to τ .
We look for a perturbative solution χ k = χ By plugging in Eq. (B3) into the right hand side of Eq. (B4) and requiring that the 0th-order perturbation doesn't induce resonances in the 2nd-order perturbation we obtain the following constraints on A(τ ) and B(τ ) This system of equations can be trivially solved by applying ∂ τ on one of the two equations and using the other to elliminate the first derivative term, leading to with B obeying the same equation. The solutions are A(τ ), B(τ ) ∼ cos(µ k τ ), sin(µ k τ ). When µ k is imaginary, the trigonometric functions become hyperbolic trigonometric functions, which for late times grow as e Im[µ k ]τ . We thus arrive at the result expected by Floquet theory, where the Floquet exponents are given by This is the expression given in the main body of text. In Fig. 21 we compare the instability bands computed in this approximation with the full numerical Floquet analysis of Eq. (B1) for the case of ϵ = 0.1 in three dimensions. As one can see, the perturbative method highlighted in this section agrees extremely well with the full computation.