Gravitational wave production from preheating with trilinear interactions

We investigate the production of gravitational waves (GWs) during preheating with monomial/polynomial inflationary potentials, considering a trilinear coupling $\phi\chi^2$ between a daughter field $\chi$ and the inflaton $\phi$. For sufficiently large couplings, the trilinear interaction leads to an exponential production of $\chi$ particles, and as a result, a large stochastic GW background (SGWB) is generated throughout the process. We study the linear and non-linear dynamics of preheating with lattice simulations, following the production of GWs through all relevant stages. We find that large couplings lead to SGWBs with a large amplitude today, of the order of $h^2\Omega_{\rm GW}^{(0)} \simeq 5\cdot10^{-9}$. These backgrounds are however peaked at high frequencies $f_{\rm p} \sim 10^6-10^8$ Hz, which makes them undetectable by current/planned GW observatories. As the amount of GWs produced is in any case remarkable, we discuss the prospects for probing the SGWB indirectly by using constraints on the effective number of relativistic species in the universe $\Delta N_{\rm eff}$.

While many preheating studies have considered four-leg interactions φ 2 χ 2 of the inflaton φ with another scalar field χ, the inflaton can never decay completely in such a case [26,[115][116][117][118]. To allow for a complete inflaton decay, interactions of the form φχ n must be included. Trilinear interactions with n = 2 are simply the most natural choice, as interactions of this type are ubiquitous in particle physics. A paradigmatic example are Yukawa interactions, which represent a coupling between a scalar and two fermions. Three-leg decay via fermions was in fact the first channel of inflaton perturbative decay considered [119,120], though it was realized later on that the inflaton also decays via nonperturbative parametric excitations in this case [121][122][123][124]. Three-leg decay via interactions with bosons leads also to non-perturbative parametric excitations [30], and this channel is actually expected to dominate over Yukawa interactions, simply due to Pauli blocking of the Fermions. Furthermore, theories with spontaneous symmetry breaking or scalar theories charged under gauge symmetries, lead naturally to trilinear vertices between bosonic species, often involving gauge fields. Even if we restrict ourselves to scalar field interactions, trilinear interactions are naturally expected in many contexts. A trilinear interaction φχ 2 between an inflaton φ and a scalar daughter field χ represents therefore a natural coupling between an inflationary sector and a matter sector, and hence we focus on this type of interaction in this paper.
Considering in particular the interaction term 1 2 σφχ 2 , where σ is a coupling of mass dimensions, we can understand the excitation mechanism underlying this interaction as follows. If following the end of inflation the inflaton oscillates coherently around the origin of its potential, the effective oscillatory frequency of a χ-field mode function χ k , is given by ω 2 k (t) = k 2 a 2 + m 2 χ + σ φ (t), with k the modulus of the comoving momentum, and a(t) the scale factor. If the bare mass of χ satisfies m 2 χ σ |φ|, the trilinear coupling will induce a tachyonic mass for χ whenever φ < 0, which happens during half of the time of each inflaton oscillation. Correspondingly, the modes within the infrared (IR) window k 2 < σ |φ| a 2 will be exponentially amplified periodically, every half-period of the inflaton oscillations. This resonant mechanism, dubbed tachyonic resonance in Ref. [30], leads to a very efficient production of χ particles, so preheating concludes within few oscillations of φ.
If the daughter field χ is the Standard Model (SM) Higgs, the question of the stabil-ity of its vacuum becomes of upmost relevance. Even though the presence of a trilinear coupling does not affect the stability of the Higgs vacuum during inflation, it can influence it during preheating [125]. If the trilinear coupling σ is large enough, the exponentially large occupation number of the Higgs created by tachyonic resonance can destabilize the vacuum. Lattice simulations have set upper bounds on the value of the trilinear coupling σ, and on the dimensionless coupling of a scale-free interaction 1 2 λ hφ φ 2 h 2 , based on demanding a stable evolution of the Higgs during preheating. According to Ref. [33], Higgs vacuum stability is assured whenever the couplings respect the bounds 10 −10 < λ hφ < 10 −8 and σ < 2 × 10 8 GeV.
In our present work we consider a daughter field χ with its self-interaction coupling strictly positive 1 . The excitation of the field modes χ k within the aforementioned IR tachyonic window provides a non-trivial time-dependent spatial configuration for χ, which in turn generates gravitational radiation. Here we study for the first time the GW production during preheating with a trilinear interaction 1 2 σφχ 2 . We consider different inflaton potentials V inf (φ) with minimum at the origin, and study the GW production during preheating in each case as a function of σ. As the periodic tachyonic excitation of the IR modes χ k resembles to some extent the excitation of gauge fields during preheating after axion-inflation [53,54], and the latter is known to source GWs very efficiently [55,56,100], we expect as well a very efficient GW production due to the trilinear coupling. We investigate therefore whether we may obtain SGWBs such that their energy density could be constrained by upper bounds on the number of relativistic species in the universe [128][129][130][131][132], hence placing constraints on σ. Furthermore there is also the possibility that due to the nature of the trilinear interaction, the frequencies of the expected SGWB could be shifted towards smaller (more observable) frequencies than the SGWBs produced in standard φ 2 χ 2 preheating [98].
The paper is divided as follows. In Section 2, we analyse the preheating dynamics and GW production for a monomial inflaton potential. In Section 3, we study the preheating dynamics and GW production for a polynomial inflaton potential. In Section 4, we discuss our results and their potential implications. Some technical details are given in the appendices.
Conventions -. We consider a spatially-flat Friedmann-Lemaitre-Robertson-Walker (FLRW) background metric, ds 2 = −dt 2 +a(t) 2 δ ij dx i dx j , where a(t) is the scale factor, and t the cosmic time. We assume summation over repeated indices. On top of the background we consider tensor metric perturbations that satisfy h ii = ∂ i h ij = 0, hence identified with GWs. We denote the reduced Planck mass by m p = 2.435 · 10 18 GeV.

Monomial potential
We start our investigation by considering a quadratic potential for the inflaton φ around the origin. While a purely monomial potential V inf (φ) ∝ φ 2 for inflation is strongly dis-favored 2 by B-mode CMB data [17], we can consider consistent inflaton potentials that flatten out towards large field values by developing a plateau, and then adopt a monomial shape around the origin only after inflation ends. These potentials are inspired by α-attractor models [134] and can be written as V inf (φ) = 1 2 Λ 4 tanh 2 (φ/M ), where M and Λ constants have dimensions of mass. In the limit |φ| M , the potential becomes Compatibility with CMB observations at CMB scales allow us to choose M = 5m p , which then leads to Λ = 0.00564m p , and from there to an inflaton mass m φ 1.6 · 10 13 GeV (see e.g. Section II of Ref. [117] for further details). Finding m φ from fitting CMB anisotropy observables with an inflationary quadratic potential V inf (φ) = 1 2 m 2 φ φ 2 at CMB scales, leads however to a very similar scale for the inflaton mass as in the α-attractor case. So, in practice, we then simply consider a quadratic potential for the inflaton φ, and supplement it with a trilinear interaction φχ 2 with a daughter scalar field χ, for which we allow a self-interaction ∝ χ 4 , where σ is a coupling with dimensions of mass, and λ is the self-interaction coupling of the daughter field. We neglect the term 1 2 m 2 χ χ 2 in Eq. (2.1) as we anticipate that the typical momenta excited k/a σ|φ| will be very large, and hence we simply assume that m χ σ|φ|. It is however necessary to include the χ 4 self-interaction to ensure that the potential is bounded from below and, consequently only stable dynamics are developed. We then identify a critical coupling so that values λ < λ c are not allowed because the potential in Eq. (2.1) would not be bounded from below. The critical value λ c defines in fact a flat direction φ = − σ χ 2 For λ > λ c the flat direction is lifted and the potential has a minimum at φ = χ = 0 where V (φ, χ) = 0. The requirement that the potential has to be bounded from below thus imposes an important relation between the two couplings of the model, σ and λ. Throughout our following analysis, we adopt a value λ = 2λ c so that our numerical dynamics are guaranteed to avoid running-away solutions.

Preheating dynamics. Lattice simulations
The dynamics of preheating depend on the curvature of the potential around its minimum, so our results hold for a class of inflationary models with similar effective inflaton potential around the origin, as long as the leading term is quadratic. As mentioned before, the monomial inflaton potential that we considered in Eq. (2.1), V inf (φ) = 1 2 m 2 φ φ 2 , is simply a very good approximation to the α-attractor potential V inf (φ) = 1 2 Λ 4 tanh 2 (φ/M ) at field amplitudes |φ| M . For simplicity in our analysis, we will fix the initial amplitude Φ i of the inflaton at the beginning of preheating (onset of inflaton oscillations), by determining the end of inflation in a quadratic inflaton potential V inf (φ) ∝ φ 2 . Following the end of inflation, the inflaton (in the form of a homogeneous condensate) starts oscillating around the origin. During the oscillatory stage, the inflaton produces copiously particles of the χ field, transferring exponentially energy into the new quanta created. In our model, preheating occurs via a trilinear interaction, which provides an effective time-dependent mass term for χ with oscillatory sign in time. Whenever φ < 0, which occurs during half the time of each inflaton oscillation, the trilinear coupling provides a negative squared mass for χ. Hence, the corresponding χ k modes will be exponentially enhanced periodically, every time that φ < 0. This mechanism, known as tachyonic resonance, leads to a very efficient production of χ particles after few oscillations of φ [30,33].
In order to understand tachyonic resonance in more detail, we start by analyzing the equations of motion of the fields φ and χ, where a(t) is the scale factor and H ≡ȧ/a the Hubble rate. It is convenient to perform a change of variables and then move into Fourier space, where we can naturally quantize the field χ as usual, writingχ withâ † k ,â k standard creation/annihilation operators that satisfy [â k ,â † k ] = (2π) 3 δ(k − k ). The evolution equation of the mode functions,χ k , can be written as a linear equation using the Hartree approximationχ 3 → 3χ χ 2 in Eq. (2.4), so that where prime denotes derivatives with respect to the new time variable, and the variance χ 2 can be related to the power spectrum of χ as follows We note that in Eq. (2.7) two dimensionless resonance parameters have naturally emerged, with Φ i the initial amplitude of the homogeneous field at the of inflation. These parameters, as we will see later, play a fundamental role in our model, as they control the intensity of the daughter field excitation. Looking at the mode frequency in Eq. (2.7), we can writẽ Figure 1: Left: Time evolution of the inflaton's volume average φ times the resonance parameter q 3 (q 3 eff = q 3 φ ), and the variance of the daughter field χ 2 for q 3 = 7. The horizontal dashed line corresponds to q 3 eff = 1/3, whereas the vertical dashed line represents the time at which the maximum value of the variance is reached, χ 2 max . Right: Time evolution of q 3 φ and the variance of the daughter field χ 2 times 3 q χ for q 3 = 50. We see that once the term 3 q χ χ 2 overtakes q 3 φ , the growth of the variance stops. For both panels, the solid (dashed) blue lines correspond to the positive (negative) part of the inflaton mean value.
As the inflaton field oscillates in a quadratic potential asφ , we see that literally during every half oscillation the inflaton amplitude becomes negative. Hence, from Eq. (2.10), we can see that every time that φ < 0,ω 2 k becomes negative for sufficiently infrared momenta with modulus below the cut-off k 2 < −m 2 χ . This leads to an exponential growth of such infrared modesχ k , which translates at the same time into an exponential growth of the variance χ 2 (or equivalently into an exponential growth of the power spectrum P χ (k) of the IR modes). Such exponential growth must actually happen in a step-like manner, as every time the inflaton amplitude is negative we expect the variance to grow, whereas every time the inflaton is positive we expect the variance to oscillate (simply due to the mode oscillations whenω 2 k > 0 for φ > 0). An overall exponential growth in steps is then expected.
In order to study the non-linear dynamics of the subsequent stages of the field excitation, we make use of lattice simulations. In this work, we use CosmoLattice [136,137]. We have considered the set of equations Eq. (2.3) and Eq. (2.4) and solved numerically a discrete version of them. The initial conditions for the homogeneous part of the inflaton field during inflation is obtained by solving the Friedmann equation together with Eq. (2.3) in the absence of the χ field (as well as in the absence of the ∇ 2 φ term). We obtain the initial amplitude and velocity from the breaking of slow roll condition, (2.11) The initial condition for the χ field are such that the homogeneous amplitude is set to zero, and standard quantum vacuum fluctuations are added on top, following the standard recipe for initial lattice power spectra, see e.g. section 7 of Ref [136] for details.
We have performed a series of simulations with N = 512 lattice sites per dimension and a minimum infrared momentum κ IR ∈ [0.2 : 0.75], where κ IR = 2π L with L the comoving length of the lattice, making sure we capture well all relevant scales of the problem. Simulations were run for different values of q 3 while considering the relation q χ = q 2 3 for q 3 ∈ [1 : 10 4 ]. Such relation between the resonance parameters simply implies that we are using λ = 2λ c , in order to be sufficiently safe from run-away solutions in the numerical dynamics. Simulations were forced to end at time m φ t f = 400, as by this time in all cases (independently of the coupling σ), the system has reached a steady state with no further energy exchange between fields.
We now explain the dynamics we identify in our lattice simulations. An initial stage is characterized by the oscillations of the inflaton, which remains as a homogeneous field for as long as the backreaction of the χ field can be neglected. In this regime, the solution to the equation of motion of the inflaton is φ(t) [26]. The effective frequency of oscillation of the modes χ k , as given in Eq. (2.7), is negative whenever the inflaton is negative, which induces an exponential growth of the modes below the cut-off κ < q 3 |φ| a . (2.12) Looking at the mode frequency squaredω 2 k in Eq. (2.10), we naturally identify an effective tachyonic resonance parameter as q 3,eff = q 3 |φ|, with q 3 the initial resonance parameter defined in Eq. (2.9). Due to the expansion of the universe, the inflaton envelope amplitude Φ(t) decays in time, and hence the effective resonance parameter decreases in time following the same decay as q 3,eff ∝ |Φ| ∝ 1 a 3/2 . Once the value of q eff drops below 1/3, the system switches to a narrow resonance regime [30], where the tachyonic excitation is no longer effective. Hence, depending on the value of the trilinear coupling, two regimes can be identified: i) For 'small' values of q 3 , q 3,eff drops below 1/3 before or soon enough after the first zero-crossing of the inflaton, and hence the growth of the variance of χ shuts down at a maximum value χ 2 max determined by the moment when q 3,eff = 1/3, as shown in the Left panel of Fig. 1. ii) For 'large' values of q 3 , the variance χ 2 grows exponentially in steps (as many as negative semi-cycles of the inflaton) up to a new larger maximum value χ 2 max , which is reached before q 3,eff drops below 1/3. We refer to this new maximum as the critical value of the variance, χ 2 max = χ 2 crit , see Right panel of Fig. 1. In this second case, just about when the variance reaches χ 2 crit , the backreaction of the daughter field into the inflaton cannot be ignored any longer. In particular, the growth of χ 2 towards χ 2 crit displaces the position of the minimum of the inflaton potential, so the inflaton is 'locked' to be around a negative value after the variance of χ reaches χ 2 crit . Numerically, we find the division between the two regimes to be at around q 3 10, i.e. the regime i) of 'small' q 3 values corresponds to q 3 < 10, whereas the regime ii) of 'large' q 3 values corresponds to q 3 > 10. Although the maximum value χ 2 max corresponds in both regimes to the end of the tachyonic resonance, in the regime i) such value is reached when the condition q 3,eff ≤ 1/3 is satisfied, whereas in the regime ii) the maximum value corresponds to the moment when the χ self-interaction overtakes the trilinear interaction for the first time in χ's effective mass mχ = q 3 φ + 3q χ χ 2 . By cancelling out these two Figure 2: Lattice results of χ 2 max for different values of q 3 (black dots). It is clearly observed that for q 3 10 the amplitude of χ 2 max is much smaller than for q 3 10, as the tachyonic growth of χ 2 is simply stopped by the condition q 3,eff = 1/3. We choose χ 2 max as the value of χ 2 averaged over one oscillation when it reaches its maximum value. For comparison, we also show the corresponding theoretical prediction for q 3 10 in Eq. (2.13) (blue dots).
terms with each other (whenφ < 0), we can obtain an analytical estimation for the critical value χ 2 crit in the regime ii) as In Fig. 1, we show the time evolution of the mean value of the inflaton times the resonance parameter, q 3 φ , and the variance of the daughter field χ 2 , for q 3 = 7 (Left panel) and q 3 = 50 (Right panel), and for q χ values that satisfy the relation q χ = q 2 3 . Using the definitions of q 3 and q χ in Eqs. (2.9), we can easily see that the critical relation between σ and λ in Eq. (2.2) translates into q χ = q 2 3 2 . Thus, the relation q χ = q 2 3 used represents simply a conservative choice to avoid a possible runaway unstable solution of the dynamics. As expected, we observe that χ 2 grows every negative semi-oscillation, whereas it oscillates during the positive ones. The variance stops growing once it reaches its maximum value χ 2 max . In Fig. 2, we plot the maximum values χ 2 max extracted from lattice simulations, as a function of the resonance parameter q 3 , and compare them against our theoretical estimation for the critical value χ 2 crit in the regime ii) with q 3 10, c.f. Eq. (2.13). Recalling the relation q χ = q 2 3 we used in our simulations, we can see that the critical variance for q 3 10 decreases inversely proportional to q 3 , as χ 2 φ /(3q 3 ), hence inversely proportional to the coupling constant σ. Although our analytical prediction Eq. (2.13) only captures roughly the real amplitude χ 2 max 3 , it captures exactly the dependence on q 3 .

Gravitational wave production
In this section, we study GW production during preheating using the recently released module of CosmoLattice dedicated to the generation of GWs by scalar fields [138]. The dynamics of GWs are governed by the linearized equation of motion of the tensor perturbations h ij ,ḧ The energy density power spectrum of GWs normalized by the critical energy density of the universe is where dΩ k represents a solid angle element in k-space. Lattice simulations follow the procedure described in Ref [139], by solving a discrete version of Eq. (2.14) (see appendix A for further details). In order to compute the GW energy density power spectrum we make use of a discrete version of Eq. (2.16). To obtain the total energy density of GWs, we integrate Eq. (2.16) over momentum k: (2.17) In Fig. 3, we show in the Top panels the matter power spectrum of the daughter field, Pχ, whereas in the Bottom panels we show the GW energy density power spectrum, Ω GW , as a function of κ = k/m φ . Spectra are measured every ∆(m φ t) = 0.2 up to a final time m φ t f = 400, assuring that the system has reached a steady state and the transfer of energy between φ and χ fields is negligible. The gravitational waves are produced in pulses coinciding with the periods of tachyonic instability of χ, with step-like jumps visible in both the matter spectra Pχ(k, t) and in Ω GW (k, t). The production of GWs ceases when the variance of the daughter field stops growing due to the end of tachyonic resonance (from that moment onward the χ k modes are no longer excited, so the spatial gradients ofχ do not evolve any further). In Fig. 4 we show the evolution of the total GW energy density in time.
We observe that it grows exponentially until it saturates to a certain value. The saturation occurs of course when χ 2 stops growing and reaches its maximum value χ 2 max . For instance, we see that GW production for q 3 = 50 stops around m φ t ∼ 15 and comparing this with the Right panel of Fig. 1, we can observe that χ 2 stopped increasing around the same time. Once the GW spectrum saturates, the SGWB redshifts freely. Since the maximum value of χ 2 depends on q 3 , the same dependence is inherited in the maximum value of the GW energy density, see Fig. 4. Another feature worth noticing is that for Figure 3: Top panels: Matter power spectrum of the daughter field χ, P χ (k), as a function of the momentum k/m φ for q 3 = 10 (Left) and for q 3 = 100 (Right). Bottom panels: GW energy density power spectrum Ω GW (k), for q 3 = 10 (Left) and for q 3 = 100 (Right). In all panels, colors run from red (early times) to blue (late times). Each spectra is measured every ∆(m φ t) = 0.2. The SGWB produced during preheating is redshifted freely once its amplitude saturates. Today's amplitude and peak frequency of the SGWB depend on the expansion history of the universe, so in order to estimate them, we need to redshift the background from the final time of our simulations t f , till the present time t 0 . This requires to know the expansion history from t f till the moment t RD when radiation domination (RD) is fully established. If we parametrize such intermediate era with an effective equation of statē w = p/ρ, we can write the redshifted frequency and amplitude today of the SGWB as (see appendix B for further details) where f is defined as with a f ≡ a(t f ) and a RD ≡ a(t RD ). In order to estimate h 2 0 Ω GW and f GW we have characterized the position κ p and the amplitude Ω (f) GW of the peak in the GW spectra at t f . While the position of the peak is characterized by the initial instability tachyonic band given in Eq. (2.12), the final amplitude of the peak Ω (f) GW decreases smoothly with the resonance parameter. We have also measured the values of a f and H f in our simulations. Very importantly, we have observed that for q 3 10 the dominant contribution to the total energy of the system comes from the daughter field, which behaves as radiation. Therefore, the system at t f is already quite close to RD, with an equation of stateω slightly smaller but very close toω 1/3. This implies that we can estimate f to be 1. Putting everything together and taking f = 1 for simplicity, we find the following parametrization for the SGWB peak today (valid only for q 3 10) .

(2.22)
We observe that the produced SGWB has a rather large amplitude but is peaked at high frequencies. We also note that the exponent of our frequency characterization coincides with that in Eq. (2.12). While a high-frequency detection program at ∼ MHz frequencies and above, has been recently put forward [140], this is yet in its infancy, and the reality of the situation is that neither current nor planned GW detection experiments are so far expected to be able to detect these high frequency SGWBs. Direct detection is therefore challenging, as it is usually the case in φ 2 χ 2 preheating scenarios also with monomial inflaton potentials. As the GW production is in any case very significant in our scenario, we may still look for another manner to probe these backgrounds, namely, using indirect constraining means. In particular, as GWs produced during preheating have shorter wavelengths than the Hubble scale at the time of production, they contribute to the energy density of the universe as radiation [113]. Hence, it might be possible to set bounds on the amount of GW energy density that can be produced, by using cosmological constraints on the total radiation density species allowed in the universe. This is typically parameterized by a deviation of the effective number of relativistic degrees of freedom beyond the Standard Model ∆N eff = N eff − N eff,SM , with N eff,SM = 3.0440 [141]. Assuming that only our SGWB contributes to the extra radiation energy density at the time of CMB decoupling, a bound on the extra radiation degrees of freedom ∆N eff sets a constraint on h 2 0 Ω  where the present energy density of photons is h 2 0 Ω (0) γ = 2.47 × 10 −5 . The Planck observations set a limit |∆N eff | 0.29 at 95% C.L. [128,129], whereas the next-generation CMB-S4 experiments will probe ∆N eff 0.06 at 2σ [130]. Additionally, the next generation of satellite missions like COrE [131] and Euclid [132] will impose bounds at 2σ on ∆N eff 0.013. A hypothetical cosmic variance limited (CVL) CMB polarization experiment could presumably get down to ∆N eff 3 · 10 −6 [142], though this does not seem experimentally realistic. In Fig. 5, we show ∆N eff for different values of q 3 considering two cases: q χ = q 2 3 (blue dots) and q χ = 10 4 (black dots). We found that the coupling that produces the largest amount of GWs is q 3 = 9, yielding ∆N eff ∼ 0.003, which is above the bound expected in a hypothetical CVL experiment, but still below the sensitivity of actually planned realistic experiments. We conclude that neither direct or indirect detection methods can probe the SGWBs predicted in the scenario studied in this section, as either the peak frequency of the GW signals is too high, or the amount of GW energy density produced is simply not enough to be constrained by ∆N eff . In order to improve this situation, two circumstances could be sought for: 1) decreasing the energy scale of inflation, so that we can shift to smaller frequencies the SGWB peak, and/or 2) increasing the amount of GWs produced relative to the field energy budget in the universe. In the next section, we introduce a polynomial potential for the inflaton which, a priori, seems to allow to satisfy both conditions. We therefore explore next to what extent GW production during preheating after polynomial inflation can improve the detectability of the SGWB due to a trilinear interaction.

Polynomial potential
In this section, we consider the most general renormalizable single-field model of inflation, where the inflaton potential is a polynomial of degree four [143] where in the second line we have parametrized the mass and trilinear coupling of the inflaton as m φ ≡ 2 λ φ φ 0 and σ φ ≡ − 8 3 λ φ (1 − β)φ 0 . The linear term in φ is absent in the potential of Eq. (3.1) since it can be absorbed by shifting the field, whereas the constant is set to the value of the cosmological constant today, which is negligible when compared to the energy scale of inflation [143]. All three terms in Eq. (3.1) are used to match the observational data, and the existence of an (near) inflection-point at φ 0 guarantees the flatness of the polynomial potential of quartic order. In this setup, λ φ is the inflaton selfcoupling, and β is a parameter that controls the flatness configuration in the vicinity of φ 0 . Notice that, if β < 0, there is a false vacuum for φ > φ 0 , and the inflaton can get stuck in this second minimum [144]. Hence, in this work, we focus on the case where 0 < β 1, and on small inflaton amplitudes, so that the inflaton field value is sub-Planckian at the time when the curvature perturbations observed at the CMB are generated, i.e., we will assume φ 0 ≤ m p [143,144]. Interestingly, the inflection-point, φ 0 , constitutes the only free parameter in this inflationary model, since all the couplings in the potential Eq. (3.1) can be expressed at its expense (for more details of this, see Refs. [143,144]). In particular, by considering the CMB measurements from Planck data [16] on the scalar power spectrum, P ζ = (2.1 ± 0.1) × 10 −9 and the spectral index n s = 0.9649 ± 0.0042, assuming these were generated N e = 65 e-folds before the end of inflation, the quantities λ φ , β and m 2 φ , can be fixed as As a lower bound φ 0 ≥ 3 × 10 −5 m p emerges naturally by demanding the stability of the inflationary scenario against radiative corrections [143,144], the inflection point φ 0 ranges in practice within the interval 3 × 10 −5 m p < φ 0 ≤ m p . The inflaton mass lies therefore in our modeling within the range 10 2 GeV < m φ ≤ 10 11 GeV. The full potential of the system inflaton-daughter field is then again 5) and in order to ensure that this potential is bounded from below, we require λ φ , λ > 0. Identifying a critical value for the self-coupling of χ as (3.6) We observe that for 0 < λ < λ c , the potential (3.5) has two minima at non-zero χ and φ where the potential is negative, V (φ, χ) < 0. For λ ≥ λ c , there is only one minimum at φ = χ = 0 with V (φ, χ) = 0. In the following analysis we use the limiting case λ = λ c , as we will see, it corresponds to the case for which the production of GWs is maximized.

Preheating dynamics. Lattice simulations
Similarly to the case of the monomial potential, we can perform a change of variables such that we are able to write the equations of motion with dimensionless quantities. We choose conveniently such that the equations of motion for φ and χ becomẽ where denotes derivatives with respect to the new time variable, and we have used again the Hartree approximationχ 3 → 3χ χ 2 . The resonance parameters emerging now in the problem are defined as The relation between, equivalent to λ ≥ λ c so that that there is only one minimum at φ =χ = 0 where V φ ,χ = 0, is given by Introducing the field decomposition from Eq. (2.6) into Eq. (3.9), we find that the mode functions ofχ obeỹ The χ field has therefore an effective time-dependent mass m 2 χ ≈ q 3φ (t) + 3q χ χ 2 , which oscillates between negative and positive values due to the oscillations of the inflaton. Whenever m 2 χ < 0 and κ 2 /a 2 < −m 2 χ , the mode frequencyω 2 k ≡ κ 2 /a 2 + m 2 χ will become negative leading to an exponential growth of the correspondingχ k modes. The power spectrum of the excited modes, Pχ(k, t), and consequently the variance of theχ field, χ 2 , will then grow exponentially in a step-like manner, every time the inflaton amplitude turns negative.
In order to follow in detail the dynamics of tachyonic resonance in this scenario, similarly as in section 2.1, we have considered Eqs. (3.8) and (3.9) and solved numerically a discrete version of these. We obtain the initial amplitude and velocity of the inflaton from the breaking of the slow-roll condition when H = 1, for different values of φ 0 :  We have performed a series of lattice simulations with N = 512 sites per dimension and a minimum infrared momentum κ IR ∈ [0.5 : 6], making sure we capture well all relevant scales. Simulations were run for different values of q 3 considering the relation q χ = q 2 3 /8 for q 3 ∈ [25 : 5 × 10 4 ], finishing at λ φ φ 0 t f = 400, which guarantees that the system has reached a steady state.
The dynamics of preheating in the polynomial scenario can be described as follows. In an initial stage, the inflaton oscillates around its minimum, like in the monomial scenario, except that the nature of the potential of Eq. (3.1) makes the first oscillations to be anharmonic. As the amplitude decays with the expansion of the universe, the oscillations become gradually more and more harmonic, as the potential becomes essentially more and more quadratic around the minimum. As within each oscillation the inflaton becomes negative during half of the oscillation, the effective mass mχ of theχ k modes becomes also negative for κ < σ|φ| a. This induces an exponential growth of such IR modes during the semi oscillations with φ < 0, as can be clearly observed in Fig. 6. As in the previous scenario we identify two regimes depending on the value of q 3 : i) For q 3 100, tachyonic resonance ends due to the expansion of the universe, since the effective resonance parameter q 3,eff = q 3 | φ | decreases with the scale factor. When q 3,eff < 1 is achieved, the system enters into a narrow regime where tachyonic resonance is not efficient anymore. This can be observed clearly in the Left panel of Fig. 6.
ii) For q 3 100, tachyonic resonance rather ends (before it becomes narrow) due to the presence of the self-interaction of the χ field. This happens when the χ self-interaction term in mχ overtakes the trilinear interaction term. In that moment the variance Figure 6: Left: Time evolution of the inflaton's volume average φ times the resonance parameter q 3 (q 3,eff ), and the variance of the daughter field χ 2 for q 3 = 50. The horizontal dashed line represents q 3,eff = 1, whereas the vertical dashed line is the time at which the maximum value of the variance, χ 2 max , is reached. Right: Time evolution of q 3 φ and the variance of the daughter field χ 2 times 3 q χ for q 3 = 500. We see that once the term 3q χ χ 2 overtakes q 3 φ the growth of the variance stops. For both panels, the solid (dashed) blue lines are the positive (negative) part of the inflaton mean value. reaches a maximal critical value, which can be expressed as The end of tachyonic resonance when χ 2 reaches its critical value (3.13) can be observed in the Right panel of Fig. 6. In Fig. 7, we show the theoretical prediction for q 3 100 of the maximum value of χ 2 given by Eq. (3.13) (blue dots) and the actual values obtained in lattice simulations (black dots). We observe that for q 3 100 the variance χ 2 grows with q 3 , while for q 3 100 it decays inversely proportional to q 3 , in correspondence with the theoretical description χ 2 crit q 3 φ 3 qχ = 8 φ 3 q 3 (where in the last identity we have introduced the saturation relation q χ = q 2 3 8 from Eq. (3.11)). While our theoretical estimation in Eq. (3.13) captures only approximately the truly measured maximum variance χ 2 max , it describes nonetheless exactly its dependence on q 3 . After the variance reaches the critical value χ 2 crit , in the following stages the system dynamics becomes fully non-linear due to the backreaction of the daughter field into the inflaton. Similarly as in the scenario with the monomial potential, due to the non-zero value acquired by the variance of χ, the inflaton is forced to oscillate around a new minimum developed in its effective potential, a minimum that corresponds to a negative value of φ.

Gravitational wave production
In this section, we study the production of GWs during preheating after polynomial inflation. We follow the same procedure explained in Sec. 2.2, tracking the real time dynamics of the GWs and measuring their energy density power spectrum Ω GW (k, t), c.f. Eq (2.16). Figure 7: Lattice results of χ 2 max for different values of q 3 (black dots). It is clearly observed that for q 3 100 the amplitude of χ 2 max is much smaller than for q 3 100, as the tachyonic growth of χ 2 is simply stopped by the condition q 3,eff = 1. We choose χ 2 max as the value of χ 2 averaged over one oscillation when it reaches its maximum value. For comparison we also show the corresponding theoretical prediction for q 3 100 in Eq. (3.13) (blue dots).

���
In Fig. 8, we show Ω GW (k, t) and the matter power spectrum of the daughter field Pχ(k). Spectra are measured every ∆( λ φ φ 0 t) = 0.2 intervals, up to a final time λ φ φ 0 t f = 400. We see that GWs are produced in pulses coinciding with the growth of the variance χ 2 of the daughter field. The production of GW stops around the time when χ 2 reaches its maximum critical value χ 2 crit . At this moment, the growth of Ω GW (k, t) ceases and the non-linearities developed in the dynamics of the inflaton-daughter system leads to a transfer of power into more ultraviolet (UV) modes in the GW spectrum.
Fixing q χ = q 2 3 /8 according to the saturation relation, we have characterized the amplitude and the peak position of the GW background as a function of φ 0 and q 3 . In the Left panel of Fig. 9, we plot the time evolution of the GW total energy density for different values of φ 0 , and observe that the total GW energy density decreases with φ 0 . On the other hand, decreasing φ 0 implies decreasing the energy scale of inflation and hence also the peak frequency of the signal today f (p,0) GW . In the Right panel of Fig. 9, we plot the time evolution of the GW total energy density for different values of q 3 , and observe that the final amplitude decreases mildly when increasing this parameter. At the same time, the dependence of f (p,0) GW on the resonance parameter q 3 is determined by the initial instability band, κ < q 3 |φ| a. Regarding the expansion history of the universe, which is relevant for the determination the f parameter to redshift correctly the signal, our lattice simulations indicate that for sufficiently large couplings q 3 10 3 the energy budget of the system is dominated by kinetic and gradient energy of the daughter field, which behaves as radiation. The universe is therefore in RD when we end our simulations, and therefore we take f = 1. Using Eqs. (2.18) and (2.19), we obtain a parametrization of the SGWB peak position and From the parametrization we observe that the maximum amplitude for the SGWB is obtained for φ 0 m p , but in that case the background is peaked at very high frequencies f GW (φ 0 = 10 −5 m p ) 10 −15 . Furthermore, the dependence on the resonance parameter q 3 just follows a similar behavior as other in parametric excitation mechanism [98], so that increasing q 3 parametrically shifts the signal to higher frequencies while decreasing its amplitude.
In conclusion, despite the fact that in the polynomial modeling we achieve partially our intention to shift towards smaller frequencies the SGWB, this can only be done at the expense of suppressing the amplitude of the signal. In order to have a large amplitude background we are forced to have it peaked, as usual, at very high frequencies, where direct detection experiments are yet far from the desired capabilities to detect a background like this one. Given the inability to detect the SGWB directly also in this model, we can still see whether constraints on a maximum trilinear interaction coupling could be placed using current and planned constraints on ∆N eff . Recalling Eq. (2.23), we translate the energy density of the GW background into ∆N eff and compare with the present/projected bounds on ∆N eff , as we can see in Fig. 10. We find that for q 3 = 100 and φ 0 = m p , we obtain a maximum value as ∆N eff 0.001, which is well above a futuristic CVL experiment, but still below the more realistic expected bounds in near future experiments. We conclude therefore that it is not possible to probe directly or indirectly the predicted signal with upcoming experiments.

Discussion
In this work, we have studied the GW production during preheating due to a trilinear interaction φχ 2 between the inflaton φ and a daughter field χ. Considering first a preheating scenario with the inflaton having a quadratic monomial potential 1 2 m 2 φ φ 2 , we added a selfinteraction term for the daughter field 1 4 λχ 4 , to ensure the total potential to be bounded from below. This imposes a relation between the resonance parameters q 3 and q χ , which control the intensity of the daughter field excitation. Only values for which the relation q χ > q 2 3 2 is satisfied avoid runaway solutions of χ. Using CosmoLattice to perform lattice simulations of the system, we find that for q 3 < 10 the growth of the daughter field variance χ 2 stops at an early moment due to the expansion of the universe, determined by when the system enters into narrow resonance. For q 3 > 10 instead, the variance grows up to a critical value χ 2 crit , determined by when the χ self-interaction overtakes the trilinear interaction. GW production follows the step-like growth of the daughter field power spectrum, and ceases once the variance of χ reaches its critical value experiments to probe this background. A large amount of GWs is generated nevertheless in this scenario, which can yield a value in terms of an effective number of relativistic species as large as ∆N eff ∼ 0.003 (corresponding to q 3 = 9). Unfortunately, this is still below the projected sensitivity on ∆N eff of next-generation CMB experiments.
We have also studied the case where the inflaton has a polynomial potential (both during preheating and inflation), whose shape is controlled by a single parameter, the value of a near inflection point φ 0 in the potential, which can range from 3 × 10 −5 m p < φ 0 < m p . This scenario might improve the situation in two directions: a) decreasing the energy scale of inflation, hence shifting to smaller frequencies the SGWB peak; b) increasing the amount of GW energy density produced relative to the field energy budget of the system. In this scenario, runaway solutions of χ are avoided if the resonance parameters satisfy the relation q χ > q 2 3 8 . Using CosmoLattice to perform lattice simulations of the system, we find that for q 3 100, tachyonic resonance ends soon after the onset of inflaton oscillations, when the system enters into narrow resonance due to the expansion of the universe. For q 3 100 instead, tachyonic resonance stops when the variance of the daughter field reaches a critical value χ 2 crit , which occurs when the χ self-interaction term overtakes the trilinear interaction term. Regarding GW production, we find that choosing the highest possible value of the inflection point, φ 0 = m p , leads to the largest GW production, with a SGWB amplitude as large as h 2 0 Ω (0) GW 5 · 10 −9 . This signal is however still peaked at large f GW 5 · 10 6 Hz, which cannot be probed by planned GW experiments. While we can shift the signal down to ∼ 10 5 Hz frequencies by decreasing φ 0 down to its minimum value compatible with CMB constraints, we can only do it at the expense of suppressing the signal amplitude down to h 2 Ω (p,0) GW 10 −15 . In terms of radiation degrees of freedom represented by the SGWB in this model, we find a maximum value as ∆N eff 0.001, which is still out of the reach of near future next-generation CMB experiments.
In summary, preheating scenarios with an inflaton monomial or polynomial potential and a trilinear coupling between the inflaton a daughter field, can lead to very efficient GW production, but the resulting GW signal is still out of the reach of current and planned direct and indirect probes of SGWBs. One direction we may explore in the future is to consider the same type of coupling in lower scale inflationary scenarios, for which the SGWB frequency peak should be shifted to smaller values, while there is no particular reason to expect that the efficiency of GW production should be decreased.

Acknowledgments
We thank Djuna Croon for interesting discussions and collaboration at early stages of the project. Our work is supported by the Generalitat Valenciana grant PROMETEO/2021/083, and by the Spanish Ministerio de Ciencia e Innovacion grant PID2020-113644GB-I00. C.C. is also supported by the grant PROMETEO-2019-083, the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 860881-HIDDeN, and partially by the FCT project Grant No. CERN/FIS-PAR/0027/2021.

A Gravitational wave evolution in the lattice
In order to evolve the h ij fields in the lattice, we follow the method introduced in [139]. We evolve a set of fields u ij with the equation of motion u ij + 3Hu ij − ∇ 2 a 2 u ij = 2 m 2 p a 2 {∂ i φ∂ j φ + ∂ i χ∂ j χ} . (A.1) The physical h ij tensor perturbations are obtained through h ij (k, t) = Λ ijkl (k) u kl (k, t) , where Λ ijkl is defined by The method relies on the fact that at any time we need to compute the GW energy density power spectrum, we fourier transform u ij fields and then project them using Eq. (A.2). This is possible because both operations commute, therefore it is not necessary to have the source projected to TT space at every time step of evolution, which would take a very long time for the simulations as the TT projection is done in Fourier space and one would need to go back and forth between real and Fourier space.
In the lattice, the TT projection is done with an equivalent lattice projector Λ (L) ij,lm written in terms of a lattice momentum k (L) i which depends on the choice of spatial-derivative, see [145] for a discussion. We use the nearest-neighbor derivative of equation (71) in [136], for which the lattice momenta is given by The SGWB spectrum actually peaks at some sub-horizon scale k p = β p a f H f , with a f = a(t f ) and H f = H(t f ). We expect the SGWB spectrum to peak therefore at a frequency f GW today with g s,t and g t the entropic and energy density relativistic degrees of freedom. We can characterize the factor G 1/4 RD by taking into account that the SM degrees of freedom above the electroweak scale amount to g s,RD = g RD = 106.75, so we write G Hz . (B.5) The redshifted GW spectrum amplitude is rad Ω GW (k) , (B.6) and using the value of h 2 0 Ω (0) rad ≈ 4 × 10 −5 , the peak amplitude of the spectrum today is