Gravitational wave spectra from oscillon formation after inflation

We systematically investigate the preheating behavior of single field inflation with an oscillon-supporting potential. We compute the properties of the emitted gravitational waves (GWs) and the number density and characteristics of the produced oscillons. By performing numerical simulations for a variety of potential types, we divide the analyzed potentials in two families, each of them containing potentials with varying large- or small-field dependence. We find that the shape and amplitude of the emitted GW spectrum have a universal feature, with the peak around the physical wavenumber $k/a \sim m$ at the inflaton oscillation period, irrespective of the exact potential shape. This can be used as a smoking-gun for deducing the existence of a violent preheating phase and possible oscillon formation after inflation. Despite this apparent universality, we find differences in the shape of the emitted GW spectra between the two potential families, leading to discriminating features between them. In particular, all potentials show the emergence of a two-peak structure in the GW spectrum, arising at the time of oscillon formation. However, potentials exhibiting efficient parametric resonance tend to smear out this structure and by the end of the simulation the GW spectrum exhibits a single broad peak. We further compute the properties of the produced oscillons for each potential, finding differences in the number density and size distribution of stable oscillons and transient overdensities. We perform a linear fluctuation analysis and use Floquet charts to relate the results of our simulations to the structure of parametric resonance. We find that the growth rate of scalar perturbations and the associated oscillon formation time are sensitive to the small-field potential shape while the macroscopic physical properties of oscillons (e.g. total number) depend on the large-field potential shape.


I. INTRODUCTION
Inflation [1][2][3], a period of accelerated expansion in the very early Universe, has been receiving increasingly strong support by several observations. Inflation makes the Universe (almost) spatially flat, and provides a mechanism for generating not only primordial curvature perturbations but also primordial gravitational waves. The amplitude of gravitational waves is model-dependent and in the simplest models it reveals the energy scale of inflation. In fact, the observations of the cosmic microwave background (CMB) anisotropies [4][5][6] detected primordial curvature perturbations, being almost scale invariant and Gaussian, and confirmed the spatial flatness of the Universe, as predicted by inflation. Thus, even though primordial tensor perturbations and small scale curvature perturbations generated during inflation have unfortunately not yet been detected, the slow-roll dynamics during inflation, responsible for the CMB-relevant fluctuations, is well understood and tested.
On the other hand, the transition from the inflationary epoch to the hot big-bang, a radiation dominated epoch, is much less known. Reheating is needed to bring the Universe into a state filled with a thermal plasma, as required by Big Bang Nucleosynthesis. Originally, the reheating process was assumed to be solely controlled by a perturbative decay of the inflaton to radiative degrees of freedom. The importance of non-linear dynamics was later recognized [7][8][9][10] and has since received significant attention, both analytically and numerically (see e.g. Ref. [11] for a review of preheating).
Among such non-linear dynamics, oscillons [12][13][14][15][16][17][18], localized long-lived objects, are now attracting increasing attention, partly because, some inflation models [19][20][21][22] preferred by observations and the axion potentials suggested by the string axiverse [23] easily lead to oscillon formation [24] 1 . Furthermore, the formation of such oscillons can be a powerful source of gravitational waves. In fact, gravitational waves might be the only tool to directly probe the dynamics and the non-linear nature of the reheating epoch.
A lot of studies on gravitational waves emitted during the formation of oscillons can be found in the recent literature [29][30][31][32][33][34][35][36]. In this paper, we try to address the following question: How much do the properties of gravitational waves emitted from the formation of oscillons, such as the shape and the amplitude of the power spectrum, depend on the potential of a source scalar field. If the resulting gravitational wave power spectrum has an almost universal shape, irrespective of the details of the scalar potential, it can be a smoking-gun for gravitational waves associated with efficient preheating and oscillon formation. If, on the other hand, the gravitational waves retain a memory of the inflaton potential, this spectral information will be useful for discriminating between different potential shapes and thus probing the inflaton potential at small field values. A similar idea was recently proposed in Ref. [32], though we systematically examine a larger variety of potential types and classify them based on their small-field and large-field shape.
The organization of this paper is as follows. In the next section, basic equations to describe oscillon formation and to estimate gravitational waves are given. In Section III, the potentials we consider in the paper are listed and classified. In Section IV, our numerical setup is given. In Section V, the results of numerical simulations on oscillon formation and gravitational waves emitted from such processes are given and discussed. In Section VI, we perform a linear analysis of the system, in order to understand the numerical results qualitatively and build physical intuition about the various contributing factors to oscillon formation. We offer our conclusions and prospects for future work in Section VII.

II. BASIC EQUATIONS
We consider a canonical scalar field coupled minimally to gravity. The relevant action is given as with ds 2 = a 2 (−dη 2 + dx 2 ), where η is the conformal time, which is related to cosmic time t as dt = a dη. The scalar field satisfies the Klein-Gordon equation, where H = aH is the reduced Hubble parameter. Throughout this work primes represent derivatives with respect to conformal time η, and = δ ij ∂ i ∂ j is the spatial Laplacian. To remove the first time-derivative, we redefine the field as ψ = aφ, leading to ψ − a a ψ − ψ = −a 3 dV dφ .
The perturbed gravitational field h ij satisfying h ii = 0 = ∂ i h ij obeys the equation, where M −2 pl = 8πG is the reduced Planck mass, Π ij ≡ T ij − g ij p is the anisotropic stress and p denotes the background homogeneous pressure. The superscript TT represents the transverse-traceless part of the anisotropic stress tensor. Defining χ ij = a h ij , we have The energy-momentum tensor of the scalar field is given as The possible components sourcing the gravitational waves are whereas the other terms are dropped when we perform the transverse-traceless projection. The details for the evaluation of the gravitational wave spectrum are given for completeness in Appendix A 2.

III. MODELS
A. Models of systematic survey for small-field shape dependence In Ref. [24], the authors studied oscillon formation during preheating in a one-parameter family of models, in which the inflaton potential is This potential class has the necessary feature to allow for the existence of oscillons: it is locally quadratic around the minimum and shallower than quadratic at larger field values. Furthermore, for α A = 1/2, one recovers the well-known axion monodromy potential, which is linear at large field values, V A ∼ m 2 M |φ|. In Ref. [29] this potential was revisited and the authors obtained a gravitational wave spectrum possessing multiple peaks, which are related to the higher harmonics present in the oscillon itself. The higher harmonics of an oscillon are related to the Taylor expansion of the scalar potential (see Ref. [29] and references therein) and thus they are determined from the features of the potential at small-field values. In order for oscillons to form after inflation, parametric resonance must be efficient enough to allow for certain wavenumbers to grow enough to probe the nonlinear structure of the potential. After the oscillons form, the oscillon itself must be supported by non-linear effects. In particular, the effects of dispersion, that would make the oscillon dissipate into radiative modes, are counter-acted by the non-linearity of the potential. In another -albeit equivalent-description, the frequency of the oscillon is smaller than the mass of the free particles in the theory, due to the fact that the potential is flatter than quadratic at larger field values. Hence, the oscillon does not decay into free particles, because they are held together by an interaction energy, which makes the decay kinematically suppressed (see e.g. Ref. [37] for quantum effects on oscillon decay). The size of the resulting oscillons and the peaks of gravitational wave spectrum can be related to both the small-field and the large-field features of the potential.
In our present study, we first investigate how the small-field shape of the scalar potential affects oscillon formation and GW emission. To see this, we introduce a series of model potentials which asymptotically behave similarly to the potential V A (φ) given in Eq. (8), with α A = 1/2, namely, V ∝ φ for φ M and whose structure around the origin is deformed by approximating V A with the Padé approximation starting with the functional form, where x := φ/M and the number in parentheses indicates the order of the approximant. These shapes are shown in Fig. 1. The lowest-order approximation is V . However, as shown in the left panel of Fig. 1, the asymptotic shape B. Models of systematic survey for large field shape dependence Many oscillon-supporting potentials exhibit a shallow growth or a flat "plateau" at large field values and a quadratic minimum, joined together through a transitional regime. These potentials are also observationally favored for inflation, since they lead to small values of the tensor-to-scalar ratio r, as required by the latest CMB measurements. We define three potential types, in order to examine the relation of the exact potential shape to the emitted gravitational wave spectra and the corresponding formation efficiency of large overdensities. We first introduce a generic four-parameter model,   The potential V X (φ) behaves as ∝ |φ| α at φ → 0, and as ∝ |φ| α−γδ at φ → ∞. We fix the behavior around the minimum to that of a massive scalar field and focus our attention on three kinds of one-parameter families in which we restrict the variation of the parameters as (α, β, γ, δ) = (2, 1, α 1 , 1), (2, α 2 , 2, 1), (2, 1, α 3 , 2/α 3 ), namely, The shapes of these potentials are shown in Fig. 2 and are chosen to model basic features of shallow inflationary potentials • V 1 describes potentials with a variety of asymptotic power-law growth at large field values.
• V 2 describes plateau potentials with varying asymptotic amplitude.
• V 3 describes plateau potentials with varying effective width around the minimum at φ = 0, defined as the field amplitude at which the potential approaches the plateau value. In other words, the potentials of the type V 3 correspond to the potentials with the same large-and small-field behavior and differ by the size of the transitional region between the two.
It is important that all potentials are locally quadratic near the origin, hence describe free massive particles at small field values, while they become "flatter" for larger field values, thus in principle supporting the formation of oscillons. Notice that they are related to each other; in fact they coincide for a particular parameter choice V 1 (φ; α 1 = 2) = V 2 (φ; α 2 = 1) = V 3 (φ; α 3 = 2) = 1 2 m 2 M 2 φ 2 /(M 2 + φ 2 ).

IV. NUMERICAL SETUP
We impose periodic boundary conditions on the boundaries of the computational domain, and we choose the initial conditions, 2 where η 0 is the initial conformal time and δφ(x) is a Gaussian random field so that the power spectrum is equivalent to that given in Minkowski spacetime, This is consistent with the Universe at the end of inflation, given that the size of the computational domain is sufficiently less than the horizon scale, so that the relevant quantum fluctuations do not "feel" the space-time curvature, at least initially. The initial scale factor is set to be a(η 0 ) = 1, and its time-evolution is governed by the Friedmann equation, where ρ φ is the averaged energy density of the scalar field φ We impose an initial condition for the amplitude of the scalar field φ 0 such that V i (φ 0 ) = m 2 M 2 /4. For the timederivative, we impose φ (x, η 0 ) = 0 and h ij (x, η 0 ) = 0. Therefore the Hubble parameter at the beginning of simulations We redefine the scalar field φ as φ = φ/M to make it dimensionless, and then the gravitational coupling constant in the right-hand side of Eq. (5) is normalised such that We assume m = 10 −2 M and the coupling parameter G is set to be G = 10 −2 . The field equations are solved with the Leap-frog method in a three-dimensional box with 256 3 grid points whose (comoving) size is chosen as L = 60 m −1 . Note that the box size is less than the initial horizon scale, L/a H −1 in . The time-interval is dη = 0.05 m −1 and we perform simulations until η = 600 m −1 . The spatial derivatives are approximated as the second-order finite differences.
For all cases, we evaluate the gravitational wave spectrum using Eq. (A14) and the power spectrum of the scalar fluctuations using Eq. (A2). Furthermore we compute the time-evolution of number density of oscillons, and their size distribution. To reduce the variance from the initial random field, possibly leading to unphysical artifacts, we perform 10 simulations for each model parameter and average them.

V. RESULTS
A. Axion monodromy VA and small-field dependence We start by performing the simulation for the axion monodromy potential, V A (φ) with α A = 1/2, in order to make contact with the results of Ref. [29]. In Fig. 3, we show the isosurface of the energy density at η = 400 m −1 with G = 10 −2 , which is the fiducial parameter in Ref. [29]. We find the formation of multiple oscillons in our box, which are stable in the time-scale of the simulation.
The time-evolution of the corresponding GW spectrum is shown in Fig. 4 in which the line color becomes thicker as time goes forward. We confirmed the existence of distinct peaks at the final time as reported in Ref. [29]. The amplitude of the GW spectrum at the end of simulation is given as Ω GW,f ≈ 10 −6 at most. If the early matterdominated phase is terminated at a time corresponding to the end of the simulation, after which the universe quickly transitions to the radiation-dominated epoch, the present value of Ω GW is given by multiplying by Ω r,0 = 0.916 × 10 −4 and (g * ,f /g * ,0 ) 1/3 . Taking these factors into account, we compute the current GW amplitude Ω GW,0 ≈ (10 −10 − 10 −11 ). This amplitude is slightly larger than the value shown in Ref. [29] because m = 10 −2 M is slightly larger than the mass used in Ref. [29] and the initial condition is also slightly different. Given these well understood differences, our results are consistent with those of Ref. [29].  After mostly recovering the results of the axion monodromy potential V A , we focus on models V (n) A deformed from V A . As we mentioned, V (n) A asymptotically behaves as V A , while the small-field shape is deformed from V A by using the Padé approximants. In this section, we investigate the impact of the (small-field) deformation on the oscillon formation and the resulting gravitational wave spectra.
In the left panel of Fig. 5, we show the time-evolution of the number of oscillons. To display them clearly, we omit the error bars resulting from averaging over 10 realisations. The simulation results contain transient objects, local over-densities that do not possess the longevity of oscillons. To remove them from our counting, we used a simple criterion of only considering over-densities whose width w exceeds a cutoff value w > w c . In the left panel of Fig. 5, we set w c = 3.5 m −1 . We discuss the selection criterion as well as the oscillon identification algorithm that we used in Appendix B. We see that for the case V A the number density of oscillons starts growing at η ∼ 110 m −1 monotonically and reach an asymptotic value of n 35 × 10 −5 m −3 . On the other hand, for the approximated cases, the number densities of oscillons grow later as the approximations get farther away from the original monodromy potential. In fact, for the case V (10) A , oscillons starts appearing at η ∼ 110 m −1 , which is almost the same time as for the case V A , and for the case V (6) A oscillons starts appearing slightly later at η ∼ 125 m −1 . In the least approximated case V (4) A ), early oscillon production is severely suppressed until η 250 m −1 , which indicates that the inflaton field does not undergo efficient parametric resonance. The instability bands of this system depend on the details of the functional shape of the given potential, as explained in detail in Section VI A. Hence, in the less approximated cases, fluctuations need more time to grow, enter the non-linear regime and ultimately form oscillons. Then, even in the approximated cases, though the number densities of oscillons slightly oscillate and have spikes, they eventually grow and reach almost the same asymptotic values as that of case V A . It should be noted that the spikes correspond to the amplification of fluctuations, leading to transient inhomogeneities, which are picked up by our oscillon detection algorithm. 3 Thus, even though the formation times were delayed in the approximated cases, the total number of oscillons at late times is largely insensitive to the exact form of the potential that we use. Apparently, at the final time of the simulation, the number of the oscillon in the least approximated case V (4) A is still growing, but we expect that it will saturate around the same number of oscillons with the other cases. This fact implies that the number is insensitive to the small-field shape of potential. In fact, after an oscillon forms, the field value in an oscillon becomes larger and the simple picture of parametric resonance breaks down. The physical properties of an oscillon are controlled by non-linearities, probing the potential beyond its local shape near the minimum.
In the right panel of Fig. 5, we show the gravitational wave spectra for each case evaluated at the end of the simulation. The amplitude and the detailed structure of the gravitational wave spectrum are found to be (slightly) sensitive to the small-field shape of potential though its shape is almost universal. The gravitational waves are produced most efficiently when the oscillons are formed and non-spherical structures appear. After that, the oscillons become spherical, which stops the production of gravitational waves. For example, the slight shift of the peaks and troughs of Ω GW for V (4) A can be attributed to the later fragmentation of the inflaton field, occurring closer to η 250 m −1 , rather than η 100 m −1 , which is the case for the other three potential shapes. Therefore the final amplitude of gravitational waves is sensitive to the time when oscillons form, which is, in turn, sensitive to the small-field shape of the potential.  In summary, in all cases but one oscillons are efficiently formed and the final number density is almost universal though the formation time is sensitive to the shape around the origin and oscillon formation occurs earlier for potentials that contain more terms and more closely approximate the monodromy potential of Eq. (8). This fact indicates that the initial growth of oscillons strongly depends on the instability bands of the given potential which are determined  by the details of the functional form around the origin. These issues will be discussed further in detail in Section VI A. Furthermore, the final GW spectra are similar in all cases (except the slight shift of the peaks and troughs which depends on the time of inflaton fragmentation) and -as expected-approach the form of V A for potentials V (n) A with larger values of n, hence potentials that approximate V A more closely. We leave a more thorough analytical and numerical investigation of the type and longevity of created oscillons in each case for future work.

B. large-field dependence: systematic studies on Vn(φ)
We move on to study the impact of the large-field shape of potential on oscillon formation and the resulting gravitational wave spectrum. In order to disentangle the various potential features, we start with the fiducial potential given in Eq. (14) and consider deformed potentials given in Eqs. (14), having for example different asymptotic behavior, while sharing a similar small-field shape.
The isosurface at ρ = 1m 2 M 2 for V 1 (α 1 = 2; φ) is shown in Fig. 6. This corresponds to the fiducial model in the present study, as shown by the orange solid line in Fig. 2. The oscillons do not exhibit a perfectly spherical shape, but we see the appearance of a number of "spikes" on each individual oscillon caused by unstable modes on small scales. The time-evolution of the power spectrum of δφ is shown in the right panel in Fig. 7. There, we initially see a well-defined range of wavenumbers that become unstable, given by k m. However, at late times we see a broad range of wavenumbers growing and the resulting spectrum is featureless. This is reminiscent of preheating in other models, where lattice simulations showed significant re-scattering between the modes, leading to a UV cascade of power (see e.g. Ref. [38] for a recent study, albeit in a different model). As a result, the gravitational wave spectrum shown in the left panel in Fig. 7 is almost flat over an order of magnitude in k space and the features reported in Ref. [29] do not emerge.
For each of the three models (V 1 , V 2 , V 3 ), we repeat the simulations by varying the model parameter α n . We show the gravitational wave spectra Ω GW (k), the time-evolution of the number density of oscillons n, and the size distribution of oscillons at the final time of simulations in the top-left panels in Figs. 8-10. The size distributions are normalised by the total number of the oscillons, N tot = dN dw dw, where w is the physical size of oscillon. As in the case of V A , we only count oscillons whose width exceeds the cutoff value w c = 3.5m −1 .
The gravitational wave spectra are insensitive to the choice of α 1 , α 2 and α 3 , namely, the asymptotic behavior of the inflaton potential does not affect the spectra, though the amplitudes are slightly different for k 0.3 m in the cases of V 2 and V 3 . These differences come mainly from delayed growth of fluctuations and their less redshift in the cases with higher amplitudes. On the other hand, the number of oscillons depends on the exact value of α 2 and α 3 , while there are no differences in the results for different values of α 1 . We find that the oscillons can form more frequently if the potential minimum is shallower (larger α 2 ) and/or wider (smaller α 3 ). This is one of major differences from the cases discussed in the previous subsection where we found that the total number of oscillons is insensitive to the detailed shape of the potential around the origin. Furthermore, we must note that the oscillon formation time is largely unchanged for various choices of α n , while the potential shape around the origin, by using different approximations for the axion monodromy potential V A , strongly affects it. We explain this behavior by using arguments based on linearized analysis of fluctuations in Section VI B.
From these findings, we can conclude that the the potential shape at φ M is highly responsible for the growth of δφ fluctuations and thus determines the oscillon formation time. Hence small changes in the shape can result in delayed or suppressed oscillon formation. Once the oscillons are formed, the amplitude of the scalar field becomes larger, even probing values of φ M , and their macroscopic physical properties such as the total number depend on the large-field shape of the potential. The stability of oscillons would also depend on it (see also Ref. [27]). However, to see the fate of oscillons in our numerical setup, we need extremely long computational time, since oscillons have been seen to survive for thousands of oscillation times. Furthermore, the life-time of oscillons in a realistic set-up would also depend on how the inflaton itself decays into other particles. The existence of multi-components oscillons has been shown in certain cases (see e.g. Refs. [28,39,40]), but cannot be considered to be a generic behavior. We leave a more in-depth study of the individual oscillon properties, including their lifetime, for future work.
The size distribution is similar in all cases, with the oscillon width w ranging between 4 m −1 w 10 m −1 . In the cases of the potentials V 1 and V 2 , the oscillon width w shows little difference between the three values of the parameter α 1 and α 2 used. For the case of the potential V 3 , the case of α 3 = 3 shows a large"spike" at w 4.4 m −1 and a secondary peak at w 6.5 m −1 , compared to the cases α 3 = 1.5 and α 3 = 2, which show a smoother and broader distribution of oscillon widths. This suggests that the width of the transition region between the quadratic minimum and the plateau of the potential strongly affects the oscillon shape.
The scalar field power spectrum P (k) is also similar between the three cases. Slight differences of the amplitudes appear in the tail of the power spectrum k 5 m in the case V 3 , in the middle of the power spectrum k ∼ m in the case V 1 , and in the low-k for V 2 and V 3 . As for the high-k tail, the case of α 3 = 3 shows larger power for k 5 m, which can be responsible for a large"spike" at w 4.4 m −1 and a secondary peak at w 6.5 m −1 in the size distribution. As for the middle-k, small α 1 gives slightly larger power spectrum for 0.5 m k 3 m. This feature can be explained by different growth rates (Floquet exponents), as investigated later in Section VI A. As for the low-k, the case of α 2 = 5/3 shows less power for k 2 m and the cases of α 3 = 2, 3 show less power for k 0.5 m, which comes from the slow growth of the scalar perturbations, again as investigated later in Section VI A. This feature can be also responsible for the slow increase of the number of the oscillons. On the other hand, the existence of the peak of the power spectrum m k 2m appears to be rather robust.

C. From monodromy to plateau potentials
Until now we have examined a variety of potential shapes and deformation, mostly focusing on two similarly disjoint families: the axion monodromy potential V A and its small-field deformations V (n) A and plateau potentials V n with varying large-field characteristics. However, these two families can be related to each other. Fig. 2 shows that the potential type V 1 does not asymptote at a finite value for φ → ∞, but rather grows as V 1 ∼ |φ| 2−α1 . For values of α 1 around 2 we see similar behavior to the "true" plateau potentials V 2 and V 3 . However, for α 1 → 1 this model resembles the axion monodromy potential of Eq. (8), in the case of α A = 1/2, which is where we focused our attention on. Fig. 11 shows the time evolution of the gravitational wave spectra Ω GW and the scalar power spectrum of the inflaton field φ. We consider several values of α 1 ranging from α 1 = 2, which is the fiducial model, to α 1 = 1, which exhibits linear growth at large field values, as in V A . We must note here that the linear growth of V 1 is not the same as that of V A , namely V 1 ∼ 1 2 m 2 M |φ| as opposed to V A ∼ m 2 M |φ|. This is necessary in order for the two models to   have the same mass at small field values, m 2 = ∂ 2 V | φ=0 . The first observation is in all cases of V 1 , both the scalar and GW spectra start growing significantly earlier than in the case of V A . By η = 100 m −1 the scalar power spectra for V 1 have largely equilibrated, regardless of the value of α 1 . The same is seen for V A at η = 150 m −1 . At the end of the simulation, η = 600 m −1 the scalar power spectra for V 1 are very similar to each other. At low k, they are identical. However a sifference is visible at k ∼ 10 m, where we see more power for larger values of α 1 . On the contrary the scalar power spectrum for V A exhibits more power at low k compared to V 1 . At large k, the scalar power spectrum for V 1 with α 1 = 1 is closer to that of V A than to that of V 1 with α 1 = 2. A more interesting observation arises for the produced GW spectrum. We see that all GW spectra for V 1 exhibit a two-peak structure, similar to that of V A , albeit at earlier times. Specifically, the two peaks are clearly visible when Ω GW peaks around 10 −6 Ω GW | max 10 −5 , while they are smooth for smaller values. The small values of Ω GW , for example in the case of V 1 with low α 1 at η = 50 m −1 correspond to GW emission by inflaton fluctuations still being in the linear regime. When non-linearities become important and oscillons are formed, the two-peak structure appears. The relation of the wavenumbers of the peaks and dips of the GW spectrum to the frequency content of oscillons was explained in Ref. [29]. However, when non-linearities occur early in our simulation, the GW spectrum evolves further towards a featureless "equilibrated" distribution, reaching a maximum value of Ω GW | max 10 −4 . This occurs later for V 1 with smaller values of α 1 . For α 1 = 1 the two peaks are clearly visible at η = 150 m −1 but have largely disappeared for η = 600 m −1 . On the contrary, the scalar power spectrum for V A enters the non-linear regime later and the features in the GW spectum appear later and remain there until the end of the simulation. The two-peak structure of the GW spectrum for V A has largely equilibrated by η 400 m −1 , so we do not expect that running the simulation for much longer will result in a complete "smearing" of the two peaks. Figure 12 shows the time evolution of the gravitational wave spectra for the three plateau potentials V n , each simulated for three values of the corresponding parameter α n . We distinctly see three largely identical regions in all cases • During the early period η 35 m −1 the gravitational wave spectrum closely follows the linear scalar power spectrum of P(k), which only shows excitation for modes with comoving wavenumbers k m.
• After the field excitations reach the point, where non-linear effects become important, we see a cascade of power in the scalar spectrum towards the UV, at η 35 m −1 . The timing is different for each pair {V n , a n }, but the overall behavior remains. The timing difference will be explained in Section VI B with the use of Floquet theory and the analysis of linearized perturbations.
• When oscillons are formed, their internal frequency content leads to a power deficiency (a "dip") in the GW spectrum, leading to a two-peak shape, as discussed in Ref. [29]. We see a two-peak structure appearing for all case of V n for 45 m −1 η 55 m −1 .
• After this rather brief period of time, the GW spectra lose memory of the two-peak shape and exhibit a single broad peak for k > 10 m −1 . This shape is largely universal, with some minor differences between some cases, discussed in the previous section.

VI. LINEAR ANALYSIS
The time-evolution of the power spectrum of the scalar field φ points towards the existence of instability bands, causing certain wavenumbers to undergo parametric resonance and exponential enhancement. In order to better understand the numerical results of Section V, we perform a Floquet analysis, by neglecting the expansion of the universe and approximating the motion of the background inflaton field as being purely periodic, without any redshifting due to Hubble drag. Since we are interested in sub-horizon scales, the static universe approximation will capture the essential dynamics.
In the static universe approximation, the equation of motion for the scalar perturbations δφ in Fourier space is given asδ where the dots represent derivatives with respect to cosmic time t and the background fieldφ satisfies Eq. (21) can be written as a matrix first-order equation This equation is of the formẋ where P(t) is a periodic matrix, whose period is controlled by the background field motionφ(t). According to Floquet's theorem, the solution of the above equation is of the form x(t) = e µt Q(t), where Q(t) is also periodic with period T. The quantity µ is the Floquet exponent. When it has a positive real part, it causes exponential enhancement of the relevant mode. In what follows, we will only focus on the real part of µ and call this the Floquet exponent. We compute the instability chart by using the algorithm presented in Ref. [11]. We solve Eq. Finally, we can extract the largest Floquet exponent, which signals the existence of instability bands when it has a positive real part. The instability is controlled by two parameters, the amplitude of the background fieldφ(t) and the physical wavenumber k.

A. Small-field dependence
We start by exploring the parametric resonance behavior of the axion monodromy potential V A with α = 1/2, along with the approximations V A . In all three cases we choose M = 10 −2 M pl ( G = 10 −2 ) and inflation ends at φ 0.4M pl . Fig. 13 shows the Floquet instability charts for each potential, as a function of wavenumber k and background field amplitude φ. As expected, the instability bands look identical for the exact and approximate potential for large field values, since all three potentials asymptote to V ∼ m 2 M |φ| for φ M . However, for φ M , the main instability bands depends on the exact potential shape. For the background field amplitude taken from the relation V = m 2 M 2 /4, the main instability band is larger for the exact potential than for the approximate ones. This leads to a stronger instability for V A as opposed to V A . This explains our finding that oscillon creation occurs earlier for potentials that are closer in shape to V A . However, the oscillons -when formed-probe the potential beyond the minimum, hence all three potentials provide identical number density of oscillons, within the accuracy limits of our simulations.
Going one step further, we solve the linear fluctuation equation for δφ on a self-consistently expanding background {φ, H}, by neglecting non-linearities and back-reaction effects. Other than that, we are using the initial conditions and parameter values that were used in the full lattice simulation. Comparing the linear fluctuation spectra to the ones obtained from lattice simulations, we can see the effects of back-reaction and oscillon formation through a deviation of the full numerical spectra from the linearized ones. Fig. 14 shows the two sets of power spectra for the three cases V A , V (6) A and V (4) A . Before going into details, we can immediately see excellent agreement between the initial and final power spectra. For the full monodromy potential V A , we see that the full numerical spectrum deviates from the linear approximation for η 110 m −1 . Referring back to Fig. 5, we see that this is approximately the time at which oscillons emerge. Soon after that point, the scalar power spectrum loses all similarity with the linearized approximation, which shows that it is governed by non-linear structure formation and field self-interactions and not by Floquet theory. Furthermore, while the spectrum seems to evolve when plotted in the axes of Fig. 14, it does not when plotted as a function of the physical wavenumber k/a. This indicates that the peak of the distribution is governed by the typical scale of inflaton fragmentation, which is related to the typical oscillon size. Similarly the peak of the distribution does not change at late times, when rescaled by the scale-factor cubed. We must note, that the range of time-slices plotted in Fig. 13 is smaller than the one plotted in Fig. 4, because here we are interested in the beginning of oscillon formation, not their late-time behavior. The case of V shows significantly different dynamics, albeit reaching a similar final state. The initial parametric resonance is weaker and the range of excited wavenumbers is smaller. However, we see a deviation of the scalar power spectrum from the linear result for η 200 m −1 , when the peak power spectrum of φ fluctuations reaches values of P(k) = O(0.01), at which time the field has a large enough amplitude to start probing non-linearities 4 . While oscillon formation in this case is not as robust as in the cases of V A and V (6) A , the inflaton field fragments in a similar way, albeit at a later time. This leads to the final scalar and GW spectra being very similar, even though the initial Floquet charts are different in three cases, exhibiting a maximum Floquet exponent that is more than 50% larger for V A than for V (4) A .
In order to make sure that the difference in Floquet charts is not merely an artifact due to the different initial value of the field, we computed the maximal Floquet exponent for 0.4 ≤ φ/M ≤ 0.75 and found that the maximum value of µ for V A is consistently 50% or more larger than the maximum value of µ for V (4) A .

B. Large-field dependence
We now examine parametric resonance in the three model potentials V n (φ) described in Section III B, in order to disentangle the contribution of the different potential features, such as the height of the asymptotic plateau. The field amplitude at the end of inflation φ end , defined as the time when = 1, can be analytically computed using the slow-roll equations of motion given in Appendix C, where φ end is also computed numerically and shown in Fig. 20. However, our lattice simulations are initialized at a later time, when the potential equals V = m 2 M 2 /4. For V 1 this corresponds to φ = M regardless of the value of the parameter α 1 . For V 2 the background field value at the start of our simulations is φ = M/ √ 2 − α 2 and ranges from 0.9M to 1.7M for the values of α 2 shown in Fig. 2. For the potential V 3 the corresponding field value is φ = M/(2 α3/2 − 1) 1/α3 and ranges from 0.8M to 1.2M for the parameter values shown in Fig. 2. By using the algorithm described below Eq. (23), we compute the Floquet charts for the three potentials and the three parameter values used for each potential. We must note again that there is a "crossover" point, where the three potentials V n have the same form V 1 (φ; . This can be used as the "prototype" potential, against which to compare any modifications, according to the parameters α 1 , α 2 , α 3 . The density plot of the instability bands for this value is shown in the far left panel of Fig. 15. The qualitative form of the 2-D Floquet charts for the other cases is similar. Instead we show the Floquet exponent as a function of wavenumber for the starting value of the background field φ at the start of our simulations. For most cases, this is close to φ = M , except in the case of V 2 with α 2 = 5/3, where the starting value is φ = √ 3M 1.7M . Fig. 15 shows that the Floquet exponent for almost all cases that we simulated is similar, leading to similar initial enhancement of the fluctuations δφ and a similar time of emergence for the produced oscillons. Furthermore, the A . Note that here we do not plot times much after oscillon formation, which are shown in Fig. 4 for VA. k is the comoving wavenumber with a(η0) = 1, as in the figures showing the results of the full lattice simulation.
Floquet exponents are larger than those in the case of the axion monodromy potential V A , leading to an earlier emergence of non-linear effects, inflaton fragmentation and oscillon formation. The case V 2 (α 2 = 5/3) is different, because the φ value is significantly larger. As the universe expands and φ red-shifts, the parametric resonance structure for this case becomes similar to the others 5 . Fig. 16 shows the evolution of inflaton fluctuations using the linear approximation for the case where all three potentials V n overlap. We see that initially the two calculations agree very well. This starts to change at η 25 m −1 , where an increase in the power spectrum at k m becomes visible. This signals the onset of oscillon formation, which has a characteristic scale of O(m). This "bump" of P(k) grows with time and eventually the true spectrum, calculated using lattice results, exhibits a broad peak centered at k m. An important factor for oscillon formation is the existence of large initial inhomogeneities, allowing the field to locally probe the non-linearities in the potential. As far as the emission of GW's is concerned, the formation of true oscillon is not important (see Ref. [36]). The important factor for GW emission is the emergence of large inhomogeneities during the preheating process. Fig. 17 shows the corresponding linear and lattice spectra for the three potentials V n and the different values of α 1 , α 2 , α 3 that we used. For V 1 we see a similar growth of fluctuations among the different choices of α 1 . The slight difference in the growth rate is exactly in line with the slight difference in µ (max) k /m which varies between 0.22 and 0.24, as shown in Fig. 15. Thus we expect the field to enter the non-linear regime in all three cases at a similar time. This is however not the case for the potentials V 2 and V 3 , where we see different linear power spectra for some parameter choices.
For V 2 we see that for α 2 = 5/3 the emergence of non-linearity occurs at η 45 m −1 . This can be attributed to the overall smaller initial instability bands. However, this is not a result of the potential itself, but rather of the initial (larger) field amplitude, which means that initially the inflaton fluctuations probe the narrower part of the Floquet chart. As the universe expands and the field amplitude red-shifts, the instability bands grow and thus the system evolves similarly for α 2 = 5/3, like it does for α 2 = 1 and α 2 = 5/7. Again we see that the onset of back-reaction and non-linearities appears when the peak amplitude of the scalar power spectrum approximately equals unity.
A similar behavior is seen for V 3 , where the case of α 3 = 1.5 shows that the true power spectrum deviates from the linear approximation later than in the case of α 3 = 3. This can again be attributed to the smaller initial Floquet bands for α 3 = 1.5, shown in Fig. 13. The instability chart for each potential, the three values of αn used for our simulations, M = 10 −2 M pl and the field amplitude φ given by V = m 2 M 2 /4. Color-coding goes as follows: α1 = 2, 1.6, 1.8 blue-dotted, red-dashed and green-solid respectively; α2 = 1, 5/7, 5/3 blue-dotted, red-dashed and green-solid respectively; α3 = 2, 3, 1.5 blue-dotted, red-dashed and green-solid respectively. The blue-dotted curves in all panels correspond to the case where all there potentials Vn overlap. We see that the Floquet exponents for Vn are significantly higher than those of VA, shown in Fig. 13, while they are similar between the various cases with the exception of V2(α2 = 5/3) and V3(α3 = 1.5). This is further discussed in the main text.

C. Monodromy and plateau potentials
Before we conclude the analysis of linearized fluctuations, let us discuss the parametric resonance behavior of the potentials V 1 (α 1 = 1) and V A , since both potentials show linear growth at large distances V ∝ |φ|. for α 1 = 2, while the overall shape remains the same. For comparison, the maximum instability exponent for V A is µ max 0.066 m, three times smaller than that for V 1 .
Comparing these results to the evolution of the GW spectrum for each case, show in Fig. 11, we immediately make the following correlation. In systems with large Floquet exponents, the GW spectrum evolves even after oscillon formation. This erases any memory of the oscillon structure that was present in the GW spectrum. On the contrary, systems where the Floquet exponents are lower, but still important enough to lead to inflaton fragmentation, result in a GW spectrum, which does not evolve significantly past the time of oscillon formation. This can be understood as follows, In the case of V n , we see that a significant fraction of the detected overdensities change their characteristics as time progresses. There are two possible explanations for this. Either the overdensities are not true oscillons, and thus they decay or fragment, or they are oscillons which move by interacting with neighboring oscillons or overdensities. In both cases, this leads to GW production long after the oscillon formation time. Furthermore, the emitted GW's are not related to the structure of the oscillons, but to random motion and thus they have a characteristic scale of k/a = O(m) and no other features. On the contrary, for V A , we find overdensities that are compatible with an oscillon height-width distribution (see Appendix B) and no other significant overdensities. Given this result, it is expected that the bulk of the energy density in the system that can source GW's is "locked" in stationary oscillons and thus GW emission will cease after oscillon formation.

VII. CONCLUSION AND PROSPECTS
Oscillon production is ubiquitous after inflation in models with a plateau-or monodromy-type potential, as is the gravitational wave production that is associated with the fragmentation of the inflaton condensate after inflation, due to efficient self-resonance. We performed a systematic study of oscillon potentials, computing the corresponding emergence of oscillons and the shape and amplitude of the emitted GW spectrum. More concretely, we explored the dependence of a potential on the properties of the gravitational waves during the oscillon formation processes such as their amplitudes and shapes, the power spectra of scalar perturbations, and the properties of the resultant produced oscillons such as the number densities and the size distributions. For this purpose, we have performed not only numerical simulations in the expanding universe but also a linear fluctuation analysis (Floquet analysis) in the static universe approximation to interpret the results of our simulations in terms of the corresponding Floquet charts. We focused on two main families of potentials (axion-monodromy and plateau potentials) and arranged a variety of potential types, where the small field dependence is changed with the same (large field) asymptotic behavior and the large field dependence is changed with the same small field behavior.
We also confirmed that the growth rate of the scalar perturbations and the associated oscillon formation time are sensitive to the small-field shape of a potential. In fact, as the potential gets far away from the axion-monodromy potential in the small field region (with keeping the same asymptotic behavior), the formation times were delayed more and more while the total number of oscillons at late times is largely insensitive to the small field behavior of potentials. On the other hand, the macroscopic physical properties of oscillons such as the total number depend on the large-field shape of a potential.
Though we have introduced a simple criterion to discriminate true oscillons from transient objects, it discrimination is quite subtle. In order to introduce more clear criterion, the lifetime of an oscillon might be the key feature. We leave a more thorough analytical and numerical investigation on longevity of created oscillons as well as on other features of oscillons through a new criterion for future work.
Finally, we have found that the shape of the spectrum and the amplitude of emitted gravitational waves are almost universal, irrespectively of the detail of potential shape. This can be used as a smoking-gun for deducing the existence of a violent preheating phase and possible oscillon formation after inflation. However, there are significant subtleties related to this issue. In both potential families, plateau potentials with varying large-field dependence and axion monodromy potentials with different small-field shape, the GW spectrum exhibits two peaks around the time of oscillon formation. In the case of monodromy potentials, the two-peak shape persists until the end of our simulations. In the case of plateau potentials, the two-peak structure is "smeared" at late times and replaced by a broad peak in momentum-space. From that we conclude that potentials exhibiting efficient self-resonance will tend to give a featureless GW spectrum, while potentials that exhibit a weaker self-resonance, and thus delayed oscillon production, will tend to give a GW spectrum that encodes the characteristics of the produced oscillons, like the internal frequencies, leading to peaks and dips in power at specific wavenumbers. If this behavior is verified for more potentials, it can act as a smoking-gun not only for deducing the existence of a violent preheating phase after inflation, but for inferring the strength of inflaton self-resonance and the oscillon frequency content.
Unfortunately, the typical frequencies of the emitted gravitational waves are around GHz and hence cannot be detected by the planned experiments like LISA. However, this frequency range has received increasing interest recently and hence new gravitational wave detectors for such frequency range have been proposed and developed [41][42][43][44][45].
That being said, oscillon formation can be the by-product not only of preheating, but of any fragmentation process of an oscillating massive scalar field with a shallow potential. This can be a modulus field or an axion in the later universe. Such a process would shift the frequency of the GW signature, possibly bringing it into the interferometer range (see e.g. Ref. [34]).
in part by JSPS Grant-in-Aid for Scientific Research Number 18K18764 and JSPS Bilateral Open Partnership Joint Research Projects. This work was supported by Mitsubishi Foundation. from which we can read off the energy spectrum of the produced gravitational waves, The dimensionless energy spectrum during the simulation is given by where ρ c = 3H 2 M 2 pl /a 2 is the critical density. The last integral is replaced by the angular average of discrete data of γ ij , where N k is the number of elements satisfying k = |k|. Eq. (A11) becomes where we used Eq. (A7). We evaluate this equation at a discrete wavenumber, k p = 2πp/L, so it can be more reduced to Assuming that the generation of gravitational waves is terminated at the end of simulation, the energy spectrum at the present time is obtained by multiplying the computed GW spectrum with a damping factor that encodes the subsequent cosmic expansion, Ω GW,0 = Ω r,0 g * ,0 g * ,f where Ω r,0 = 0.916 × 10 −4 , the subscript 'f' represents the quantity evaluated at the end of simulation, η = η f , and we used In the present time, the frequency corresponding to the box size is given by where k phys,f is the physical wavenumber at the end of simulations, a i , a f , a 0 are the scale factor at the initial time, the end of simulations and the present time, respectively. In addition, we assume that the inflaton decays into radiation at a = a f , leading to a reheating temperature T R .
Reinstating the constants and c, the GW frequency is written as which falls outside of the observable frequency range of ground-based interferometers for m ∼ 10 −5 M Pl . Future attempts of high frequency GW detection are becoming increasingly interesting, as they would open a unique observational window into preheating. In fact, there are proposed experiments for the detection of GWs with such high frequencies [41][42][43][44][45], even though their sensitivity needs to be improved significantly in order to provide feasible detection opportunities.

Appendix B: Oscillon identification
We identify oscillons in the computational domain using the following algorithm: 1. Once we obtain the energy density, ρ(x i , η), from the simulations, the average energy density in the computational box is computed as where N is the number of grid.
3. For each region, we find out the point, x c , where the energy density is maximum, and define ρ core = ρ(x c , η).

4.
For each region, we again find out the region centered around the peak location, x c , where ρ(x, η) ≥ ρ core /e, and then we regard the region as an oscillon with the comoving volume V osc .
5. Finally, we define the effective physical size of the oscillon as w = a(6V osc /π) 1/3 , which is the diameter of a sphere whose physical volume is a 3 V osc .
Using this algorithm and given enough data, we can construct a frequency distribution (histogram) for the number of oscillons and for their energy. The distribution functions of the oscillons' size and energy can be estimated from the computed frequency distributions. Given a distribution function f (x), the number f i in a finite bin x i ≤ x < x i+1 , is calculated as If the size of bin is sufficiently small, we estimate The above method will identify any sufficiently large local overdensity. These might or might not be oscillons: localized long-lived scalar field structures that oscillate in time. Fig. 18 shows the overdensities identified by the procedure described above for the case of the monodromy potential. We see that the vast majority of overdensities follow a clear relation between the oscillon width and height, where wider oscillons tend to have a smaller central density. This is reminiscent of models, where the height-width relation was derived semi-analytically. We do however see a significant amount of overdensities that have a small width and a small height, falling clearly below the "main sequence" of oscillons. These are transient overdensities and should not be counted as oscillons. By eliminating these overdensities, whose height-width characteristics put them below the main sequence of oscillons, we significantly reduce the numerical noise in our calculation of oscillon number density. Fig. 5 shows the resulting oscillon number, after such small overdensities have been eliminated from our count.
For our calculation the exact stability properties of the emerging oscillons are not important, since we only require the oscillons to form and live long enough for GW's to be emitted [36]. In some cases of potentials V n , the distribution of oscillons that we discovered showed two somewhat disjointed regions, making the identification of real oscillons and transient overdensities even more difficult. Possible features, such as fragmentation of unstable oscillons into smaller, more stable ones, can lead to a late-time increase in oscillon number density. Furthermore, overdensities starting close to the theoretical oscillon height width curve can relax to a stable oscillon configuration over time, as was shown e.g. in Ref. [18].
Given the rich dynamics of non-linear field theories, We leave a detailed investigation of the properties of the produced oscillons and other localized over-densities in the various potentials that we examined for future work (for recent work on the lifetime of oscillons see Refs. [27,46]).
A . While the vast majority of oscillons falls around a well defined height-width relation, we see the existence of several overdensities with small width and height, especially for V (4) A . These are numerical artifacts and have been removed before we compute the number density of oscillons. The evolution of the number density of oscillons for the case of the potential V2 with α2 = 5/3. The solid (dashed) curves correspond to counting all over-densities as oscillons (applying the selection criterion w > wc = 3.5 m −1 ). The three different colors correspond to different initial field amplitude φ0. Right: The normalized distribution of oscillon widths for the three cases of φ0.
Before we conclude, we revisit Figures 9 and 10, which show the time evolution of the comoving oscillon number density n for V 2 and V 3 respectively. In most cases the oscillon number density asymptotes to a constant value by the end of the simulation at η = 600 m −1 . However, two cases show a late-time growth of the oscillon number density: the case of V 2 with α 2 = 5/3 and V 3 with α 3 = 1.5. The common feature of these two cases is that the field value at the start of the simulation, defined through V = m 2 M 2 /4, is higher than the corresponding initial field value for all other cases (φ 0 ∼ M ). This means that the system must redshift more before entering the main instability band. Two interesting observations arise. One is that when the simulation starts at a lower field value φ = M or φ = 0.8M the time evolution of the oscillon number density exhibits a smaller initial growth of n (which is due to less transient overdensities). The second concerns the final state of the oscillons; even though the total number density is different, the distribution of oscillon widths is nearly identical for all cases. Furthermore, the solid and dashed curves, corresponding to counting all over-densities or using a width cutoff criterion, converge at late times. Overall, while the time evolution of oscillons depends both on the cutoff criterion and on the initial field value, their distribution is identical, since it depends solely on the structure of the scalar potential. We can now summarize the three potentials V 1 , V 2 , V 3 . Fig. 20 shows the field and the potential at the end of inflation, as computed by numerically solving the background equation of motion.