Characterizing the post-inflationary reheating history, Part II: Multiple interacting daughter fields

We characterize the post-inflationary dynamics of an inflaton $\phi$ coupled to multiple interacting daughter fields $X_n$ ($n=1,\dots N_d$) through quadratic-quadratic interactions $g_n^ 2\phi^2 X_n^2$. We assume a monomial inflaton potential $V(\phi) \propto |\phi|^p$ ($p \geq 2$) around the minimum. By simulating the system in 2+1-dimensional lattices, we study the post-inflationary evolution of the energy distribution and equation of state, from the end of inflation until a stationary regime is achieved. We show that in this scenario, the energy transferred to the daughter field sector can be larger than $50\%$, surpassing this way the upper bound found previously for single daughter field models. In particular, for $p \geq 4$ the energy at very late times is equally distributed between all fields, and only $100/(N_d + 1) \%$ of the energy remains in the inflaton. We also consider scenarios in which the daughter fields have scale-free interactions $\lambda_{nm} X_n^2 X_m^2$, including the case of quartic daughter field self-interactions (for $n=m$). We show that these interactions trigger a resonance process during the non-linear regime, which in the single daughter field case already allows to deplete more than $50\%$ of the energy from the inflaton for $p\geq 4$.


Introduction
The inflationary paradigm describes an early phase of accelerated expansion of the universe [1][2][3][4]. Inflation provides a solution to the initial condition problems of classical cosmology, and generates a spectrum of quantum fluctuations that seeds the later structure formation. In the simplest model realizations, the accelerated expansion is sourced by the vacuum energy of a scalar field φ (the inflaton) in a slow-roll regime. Cosmic Microwave Background experiments such as Planck [5] or BICEP [6] have been able to significantly constrain the range of observationally-viable scalar potentials.
Inflation must be followed by a reheating stage, during which the universe must transition to a radiation-dominated thermal state before the onset of Big Bang Nucleosynthesis at T ∼ 1 MeV [7][8][9][10]. Relevant questions concerning this phase are how the energy stored in the inflaton is transferred to other light fields and eventually to the Standard Model particles, or the exact post-inflationary evolution of the equation of state. In fact, a complete characterization of the expansion history after inflation is crucial for reducing the theoretical uncertainty in the predictions of CMB observables for specific inflationary models [11][12][13][14][15]. Details of reheating depend strongly on the physics model under consideration, but its early stage is typically characterized by a non-perturbative, out-of-equilibrium excitation of field fluctuations called 'preheating' [16][17][18][19][20][21][22][23][24][25][26][27]. Although these fluctuations can be treated with a linearized approximation immediately after inflation, non-linearities become relevant at later times, so one must resort to real-time lattice simulations in order to fully capture the dynamics of this stage [20,21,26]. For reviews on (p)reheating see e.g. [28][29][30][31][32], while for a review on lattice techniques for the simulation of the early universe field dynamics see [33].
In this work we consider inflaton potentials that can be approximated as a monomial function V (φ) ∝ |φ| p (p ≥ 2) after inflation. In this case, the post-inflationary oscillations of the inflaton homogeneous mode generate an effective equation of statew hom ≡ (p − 2)/(p + 2) during the linear regime [34]. If the inflaton is decoupled from other fields, it fragments through a process of self-resonance for p > 2, which gives rise to a radiation-dominated universe ∼10 e-folds of expansion after the end of inflation [35,36]. If p = 2, the inflaton fragments instead due to gravitational effects at much longer time scales [37].
In any case, a successful reheating process requires an efficient depleting of energy from the inflaton, which can be naturally achieved by coupling it to other fields (which we refer to as 'daughter' fields). One possibility is to couple the inflaton to a daughter scalar field X through a quadratic-quadratic interaction g 2 φ 2 X 2 , with g 2 a dimensionless coupling strength. This is e.g. the leading term in the interaction between charged scalars and gauge fields [38,39], and has the advantage of not introducing new dimensional scales in the theory (which facilitates lattice simulations). In this case, the field fluctuations get excited through a process of broad parametric resonance, and the inflaton fragments instead ∼3−4 e-folds after the end of inflation. The evolution of the energy distribution and equation of state during preheating in these models has been studied in detail with lattice simulations, see [40][41][42][43]. The case of an inflaton coupled to a daughter field through a trilinear interactions has also been studied on the lattice, see [44]. Other works have studied the role of non-minimal kinetic terms during preheating in different scenarios: in the context of DBI inflation in [45], in the context of α-attractor scenarios in [46] (see [47] for a semi-analytical study), and in the context of multifield models with non-minimal couplings to gravity (which appear as non-minimal kinetic terms in the Einstein frame) in [48,49] (see also [50][51][52] for semi-analytical studies). The fields eventually achieve a turbulent regime at late times, see [53,54].
While most lattice studies of (p)reheating have focused on the linear and early non-linear stages, in the Letter [55] we instead characterized the entire evolution of the energy distribution and equation of state after inflation, from the end of inflation until the achievement of a stationary regime. We considered an observationally-viable inflaton potential that behaves as V (|φ|) ∼ |φ| p around the minimum, and coupled the inflaton to one (effectively massless) daughter field through a g 2 φ 2 X 2 term. By simulating the post-inflationary dynamics in 2+1 dimensions, we were able to parametrize how the energy density distributes between its components at late times as a function of p and g 2 , as well as the final values for the equation of state. We found that the fraction of energy transferred to the daughter field is always negligible for p < 4, while it is at most ∼ 50 % for p ≥ 4. We continued this work in Ref. [56] (which we refer to as Part I from now on), in which we expanded the results of our Letter [55], as well as generalised our analysis to a class of inflaton potentials that have a 'displaced' minimum V (|φ|) ∼ |φ − v| p (with v ≥ 0). Moreover, by using our information of the equation of state evolution, we were able to obtain exact predictions for the inflationary observables n s and r in the model under consideration, for those cases in which the universe ends up in a radiation-dominated state.

Part II: Multiple interacting daughter fields
Notably, most (p)reheating studies (including our Letter [55] and Part I [56]) have considered scenarios in which the inflaton is coupled to at most one daughter field. However, the correct physics model at high energies may well contain many scalar fields. For example, one could consider models with multiple inflaton fields. In this work we consider instead a different scenario, in which the stage of inflation is still generated by a single inflaton φ, but the inflaton is coupled to multiple daughter fields during the subsequent phase of reheating. The existence of multiple daughter fields does doubtlessly change the post-inflationary evolution of the energy distribution and equation of state with respect to single daughter field scenarios, as they give rise to additional channels through which the energy of the inflaton can be extracted. Knowing the exact evolution of the equation of state after inflation is essential, for example, to make accurate predictions for the CMB observables n s and r, as shown in Letter [55] and Part I [56] for the single daughter field case. The post-inflationary equation of state is also important for e.g. determining the exact redshift of a GW signal produced during inflation or preheating until today. Moreover, the energy transfer to daughter fields can have a significant influence on the produced baryon asymmetry of the universe from non-thermal leptogenesis, see e.g. [57]. In addition, systems of multiple daughter fields can have further interesting phenomenological consequences. For example, during (p)reheating a daughter field gets excited during the linear regime at a specific momentum scale, set by the strength of its coupling to the inflaton [58]. As shown in [59], in systems with multiple daughter fields these different scales can be imprinted in the produced spectrum of gravitational waves, which features a 'stairway' pattern that potentially allows for particle coupling spectroscopy (see also [60] for a previous study with multiple daughter fields with equal couplings).
The aim of this paper (which we refer to as Part II ) is, therefore, to study the postinflationary dynamics of models with an arbitrary number of daughter fields. This work constitutes a direct continuation of the research carried out in Letter [55] and Part I [56], which considered single-daughter field scenarios. As in these past works, we will characterize the post-inflationary evolution of the energy distribution and equation of state, from the end of inflation until the establishment of a stationary regime. Special emphasis will be put on describing how the energy gets distributed at very late times. Our analysis will be based both on a numerical analysis of the linearized field equations under a Hartree approximation, as well as on 2+1-dimensional lattice simulations of the system carried out with the code CosmoLattice [61]. 1 Regarding our lattice-based analysis, we will start by studying the post-inflationary dynamics of an inflaton coupled to multiple daughter fields X n (n = 1, . . . , N d ) through quadratic-quadratic interactions g 2 n φ 2 X 2 n , where the coupling strength of each daughter field can be different. After that, we study again the case of one daughter field X, but include now a quartic self-interaction λX 4 into our analysis, which was not taken into account in Refs. [55,56]. Finally, we will consider different examples of multidaughter field theories with scale-free interactions of the type λ nm X 2 n X 2 m , which include both quartic self-interactions (when n = m) and quadratic-quadratic interactions between different daughter fields (when n = m).
The structure of this work is as follows. In Sect. 2 we present the details of our model set-up. In Sect. 3 we investigate the early stage of preheating by means of a linearized analysis of the field equations in the Hartree approximation. In Sect. 4 we deploy lattice simulations to study the non-linear regime of the field dynamics in the different scenarios explained above. In Sect. 5 we summarize and discuss our results.

Inflaton potential
Let us consider an inflaton field φ with monomial potential around the minimum, where λ is a dimensionless parameter, µ is an energy scale, and p is an arbitrary coefficient obeying p ≥ 2. Cosmological observations rule out the case of potential (2.1) sustaining inflation at all amplitudes [5]. However, we can still consider observationally-viable potentials that behave like (2.1) around the minimum but flatten at larger amplitudes. One example is the α-attractor T-model [62], where Λ and M have dimensions of energy. The ratio Λ 4 /M p can be fixed so that we recover the monomial potential (2.1) in the limit M → ∞. The observed value of the tensor-to-scalar ratio (r < 0.036 at 95% confidence level [6]) imposes the upper bound M (8.5 − 9.5)m pl when the number of e-folds from the horizon crossing of the pivot scale until the end of inflation is fixed to N k = 60 (the exact bound for M depends on p). Moreover, by fitting the theoretical prediction for the scalar amplitude to the observed value A s = 2.1 × 10 −9 , one can determine a relation between the model parameters as Λ = Λ(p, M, N k ) (the explicit expression is written in Eqs. (A5) and (A6) of Part I [56]).
Slow-roll inflation takes place at large field values, and ends approximately when the condition V ≡ m 2 pl V 2 ,φ /(2V 2 ) = 1 holds, at the amplitude The inflaton then starts oscillating around the minimum of the potential. Note that if M 1.65m pl , we have φ * < φ i for all values of p ≥ 2, where φ i ≡ M arcsinh (p − 1)/2 is the inflection point of the potential. For these values of M , the oscillations always take place in the positive-curvature region of the potential, and can be approximately parametrized during the initial linear regime as [34] where φ ≡ φ(t ) φ * is the field amplitude at a time scale t = t close to the end of inflation, and F(t) is an oscillatory function with oscillation period The period is constant for p = 2, but time-dependent for p = 2. The effective equation of state, i.e. the ratio between the oscillation-averaged pressure and energy densities, is approximatelȳ (2.6)

Daughter field sector
The objective of this work is to study the post-inflationary dynamics of an inflaton φ coupled to multiple daughter scalar fields X n , with n = 1, 2, . . . N d and N d ≥ 1. We will consider different particularizations of the following potential, where V inf (φ) is the potential sustaining inflation (either (2.1) in Sect. 3 or (2.2) in Sect. 4). Each daughter field is coupled to the inflaton and other daughter fields via quadratic-quadratic interactions, with g 2 n and λ nm (= λ mn ) denoting the corresponding dimensionless coupling strengths. The last term also contains quartic self-interactions of the daughter fields when n = m: in this case we will use the short notation λ n ≡ λ nn .
Mimicking the procedure in Part I [56], it is convenient to work in natural variables, defined for field amplitudes and spacetime coordinates as This way, the amplitude and oscillation period of the 'natural' inflaton ϕ during the linear regime are constant in 'natural' time u. The equations of motion can then be written as where we have fixed a(t * ) = 1 at the end of inflation, F (u) ≡ 6(p−4) (p+2) 2 (a /a) 2 − 6 2+p (a /a) ∼ u −2 is a function that becomes subdominant after a few inflaton oscillations, and we have defined the following time-dependent functions, Here, q (n) * is the (initial) resonance parameter of the daughter field X n , andq (n) the corresponding time-dependent (effective) one. We have defined σ (nm) * andσ (nm) for the daughterdaughter interactions in an analogous way. At the end of inflation we haveq (n) = q (n) * and σ (nm) = σ (nm) * , but both functions evolve in different ways depending on p: they decrease for p < 4, grow for p > 4, and remain constant for p = 4.

Linearized analysis of the field equations
The inflaton oscillations may lead to a strong growth of fluctuations for either the inflaton (in a process of self-resonance), the daughter field (in a process of parametric resonance), or both. In order to illustrate this, let us expand the fields as ϕ( y, u) ≡φ(u) + δϕ( y, u) , whereφ is the homogeneous component of the inflaton (note that we haveχ n 0 at the end of inflation). Under the approximation F = 0, the homogeneous part of the inflaton obeys the equationφ + |φ| p−2φ 0, whose solution isφ = cos(u) for p = 2 andφ cos(βφu) with βφ ∼ 1 for p > 2. The fluctuations can be described by their mode equations in Fourier space, is the 'natural' resonance momentum. The effective frequencies ω k,ϕ andω k,χn may vary non-adiabatically when the inflaton homogeneous mode crosses the minimum of the potential, which leads to a strong growth of the field fluctuations. More specifically, the mode equations (3.3) and (3.4) allow for solutions of the form δϕ k ∼ e µ k u and δχ n,k ∼ e ν k u , with µ k ≡ µ k (p) and ν k ≡ ν k (p,q) the corresponding Floquet indices. For certain values of p andq (n) , the real parts of the Floquet indices are positive, thus leading to the following two resonance phenomena: • Inflaton self-resonance: The structure of (narrow) resonance bands is shown as a function of p in the left panel of Fig. 1. For a given choice of p, the dominant band is the one of lowest momenta, with Re[µ k ] 0.036. Note that there is no self-resonance for p = 2.
• Parametric resonance of the daughter field: The structure of resonance bands is shown in the right panel of Fig. 1 (we show the case p = 4, but very similar charts can be depicted for other values of p, see Fig. 7 of Part I [56]). The regime of strongest resonance corresponds toq 1 (broad resonance), while forq 1 the resonance is very weak (narrow resonance). The type of resonance can change as the universe expands: for p < 4 an initially broad resonance becomes narrow at later times, while if p > 4 an initially narrow resonance becomes broad at later times. If p = 4, the type of resonance never changes. The maximum Floquet index for broad parametric resonance is Re[ν k ] 0.26. Thus, parametric resonance is much stronger than inflaton self-resonance.
The linear regime of resonant excitation ends when the energy of the fluctuations becomes comparable with the one of the inflaton homogeneous mode. In Part I [56] we computed an analytical estimate of the backreaction time u br when the equality ρ |φ| ρ δf holds exactly, for both fields f = ϕ, χ.

Termination of resonant growth
Let us perform a Hartree or mean-field approximation [17,21,22,24,63] to the field equations, which allows to partially capture the early backreaction effects of the λ nm X 2 n X 2 m interactions during the linear regime. Following the procedure of Refs. [28,63], we substitute δf 2 → δf 2 and δf 3 → 3 δf 2 δf in the fluctuation equations in position space, where f = ϕ, χ labels the fields and δf 2 ≡ (2π 2 ) −1 dkk 2 |δf k | 2 . This approximation neglects the couplings between modes of different momenta, but captures sufficiently well the field dynamics during the linear regime as we shall see. We obtain the following equations for the field modes which incorporate the corrections generated by the λ nm X 2 n X 2 m interaction terms to Eqs. (3.3) and (3.4). Note that, if we set k = 0 and δϕ 0 ≡φ in Eq. (3.5), we obtain the corrected equation for the inflaton homogeneous mode.
The field variances contribute to the effective masses of both fields, and grow exponentially during the initial inflaton oscillations. If they become large enough, they trigger the decay of the inflaton homogeneous mode due to backreaction effects, which terminates the resonant growth. However, they can also block the growth of the field modes before that happens. As an example, consider the case of one daughter field in broad parametric resonance: the first situation happens when q * σ * (as long asq > 1) via the coupling ∼q δχ 2 in Eq. (3.5), 2 while the second situation happens when q * σ * and the effective mass becomes of the order of the term responsible for the resonance, 3σ δχ 2 ≈qφ 2 . The maximum 2 In the following we will use the shorter notation of χ ≡ χ1, q * ≡ q (1) * and σ * ≡ σ (1) * for single daughter field scenarios. Figure 2. Left: Evolution of the daughter field variance for N d = 1, p = 4 and q * = 30, obtained by 1) solving the linearized equations under the Hartree approximation (continuous lines), and 2) simulating the system in the lattice (dashed lines). We consider three different self-coupling parameters: σ * = 3 · 10 0 , 3 · 10 4 , and 3 · 10 8 . Right panel: Variance of the daughter field attained at late times for different choices of σ * , computed with both the Hartree approximation and with lattice simulations. The black dashed line indicates the estimate δχ 2 variance attained in this case is δχ 2 max ≈ q * /(6σ * ) (this result can be obtained by settinḡ ϕ 2 φ 2 Tφ ≈ 1/2). In these cases, the inflaton homogeneous mode survives and continues to dominate the energy budget (at least as long as the Hartree approximation is valid).
We illustrate this in the left panel of Fig. 2, where we show the evolution of the variance of one daughter field for p = 4, q * = 30, and three different self-coupling parameters: σ * = 3 · 10 0 (< q * ), 3 · 10 4 (> q * ) and 3 · 10 8 (> q * ). The solution is obtained by solving numerically Eq. (3.6) together with (3.5) for the homogeneous mode k = 0 (we have ignored the growth of the inflaton fluctuations δϕ k with k > 0, as inflaton self-resonance is a much weaker effect). We compare each solution with the result from a lattice simulation, which takes all nonlinearities into account (see Sect. 4 for more details). For σ * = 3 · 10 4 and σ * = 3 · 10 8 , the growth of the variance saturates at δχ 2 max = O(0.1)q * /(6σ * ), where we have included a correction factor of order O(0.1) to our naive analytical estimation. For σ * = 3 · 10 0 we have σ * < q * and the results obtained by the Hartree approximation and the lattice differ more strongly (see left panel). In this case, the inflaton homogeneous mode decays due to backreaction effects, which is not fully captured by the Hartree approximation. The dynamics at later times can only be properly studied with lattice simulations. Parametric resonance terminates before the growth of δχ 2 saturates due to its own effective mass, so the estimation The right panel of Fig. 2 shows the saturated variance δχ 2 max for q * = 30 and different choices of σ * , obtained with both a Hartree approximation and lattice simulations. The vertical dash-dotted line indicates the case σ * = q * . For σ * /q * > 1 the saturated variance approximates quite well the estimation δχ 2 = 0.1 · q * /(6σ * ) (in particular for very large ratios). For σ * /q * ≤ 1 the resonance is terminated by backreaction effects from the daughter field modes onto the homogeneous inflaton condensate and the variance saturates roughly at δχ 2 ∼ (p − 1)|φ| p−2 /q * . The results from the lattice and the Hartree approximation differ more strongly.
Analytical estimations for δχ 2 n max (n = 1, . . . N d ) can also be obtained for two or more interacting fields, but expressions become significantly more complicated. For illus- trative purposes, let us consider a two-daughter field scenario. Fig. 3 shows the evolution of the daughter field variances, obtained again by solving the field equations under the Hartree approximation and with lattice simulations. The left panel shows a scenario in which σ In this case, the role of the daughter-daughter interaction is negligible, so the evolution of each individual variance is similar to the single-daughter field case: they saturate at roughly the value δχ 2 (2) * ). In this case, the variance of the first daughter field (the one with largest q * ) saturates roughly at the value , in agreement again with the single-daughter prediction. However, the evolution of the second daughter field gets affected by the fast growth of δχ 2 1 due to the strong coupling ∼ σ (12) * χ 2 1 χ 2 2 . In fact, we find that δχ 2 2 temporarily saturates when the variance of the first daughter field attains the value δχ 2 1 ∼ q (2) * /(2σ (12) * ) ( δχ 2 1 max ). However, we observe that δχ 2 2 starts growing again after a while in both the Hartree approximation and the lattice.

Energy distribution after inflation: lattice simulations
We now investigate the post-inflationary evolution of the energy distribution with lattice simulations, going beyond the linearized analysis of Sect. 3. The simulations have been carried out with the publicly available code CosmoLattice [61]. We have done simulations in 2+1 dimensions, which allow to properly capture the very late-time regime of the field evolution. In Appendix A of Letter [55] we carried out lattice simulations of single field scenarios (described by potential (4.1) below) in both 2+1 and 3+1 dimensions and compared explicitly the output: this way, we showed that (2+1)-D simulations mimic the dynamics in (3+1)-D very well, at the level of both volume-averaged quantities and field spectra 3 . Furthermore, we checked that this holds true for simulations with multiple daughter fields as well. The discrete field equations have mainly been solved with the velocity-verlet integrator of 2nd order of accuracy implemented in the code, although for some model parameters we have required the 4th order one. Depending on the particular scenario, we have used lattices between N 2 = 128 2 and 1024 2 points.
In the following subsections we consider different particularizations of potential (2.7): Sect. 4.1 reviews the scenario of one daughter field coupled to the inflaton through a quadraticquadratic interaction, already studied in detail in Part I [56]. Sect. 4.2 considers the case of multiple daughter fields coupled to the inflaton with different strengths. Sect. 4.3 analyzes the case of one daughter field with a quartic self-interaction ∼λX 4 . Finally, in Sect. 4.4 we consider multi-daughter field scenarios, but including now interactions between them of the type ∼λ mn X 2 n X 2 m (which incorporates quartic self-interactions when n = m). In all cases we take the α-attractor T-model potential (2.2) for the inflaton, and consider different values of p. We have fixed M = 10m pl , so that φ i > φ * and the inflaton oscillates in the positively curved region of the potential (note that this value of M is in slight tension with the upper bound for the tensor-to-scalar ratio of the inflationary perturbations [6], M (8.5 − 9.5)m pl , but the dynamics is very similar for smaller values of M as long as φ i φ * ). The amplitude of the plateau is fixed through the relation Λ = Λ(p, M, N k ) for N k = 60 (see comment after Eq. 2.2). 4 In the following analysis, we will characterize the evolution of the energy distribution in terms of 'energy density ratios' (or simply 'energy ratios') ε i ≡ E i / i E i , defined as the relative contribution of each (volume-averaged) energy density component i to the total (volume-averaged) energy density. Different contributions include the kinetic and gradient energies of each field (defined as ε f k and ε f g for f = ϕ, χ respectively), and the different terms of the potential (defined as ε t p with t labelling the corresponding term). Expressions for the energy ratios ε i and equation of state w are given in Appendix A. Bared quantitiesε i andw denote the corresponding oscillation-averaged expressions. 5 4.1 One daughter field without self-interaction: The case of one daughter field without quartic self-coupling was studied extensively in Part I [56], so here we simply summarize our main results. The values attained by the energy 4 Note that this gives rise to a slight inconsistency, as the exact value of N k can only be determined with knowledge of the full post-inflationary evolution of the equation of state, which we obtain from a lattice simulation in which N k has a priori been fixed. This was solved in Part I [56] with an iterative procedure, in which several lattice simulations allow to determine N k up to a factor O(10 −2 ) (see Section 5.A of that paper for more details). However, we found that the post-inflationary dynamics remains basically unchanged between different simulations, so here we have fixed N k = 60 for simplicity. 5 In the Figures presented in this paper, the oscillation averages of the energy ratios and equation of state have been obtained by means of a mean filter. ratios and effective equation of state at late times are given in Table 1 for different choices of p and q * (by setting N d = 1). Let us briefly consider each case:

Final energy ratios for
The inflaton does not get excited via self-resonance, but the daughter field does via broad parametric resonance as long asq ≡ q * a −3 1. If the stage of broad resonance takes long enough, backreaction effects trigger the decay of the inflaton homogeneous mode, and the equation of state deviates fromw =w hom (= 0) towardsw =w max ( 1/3). In any case, the production of daughter field fluctuations terminates onceq 1. As fluctuations dilute as radiation, the inflaton homogeneous mode (which dilutes as matter in the present case) eventually dominates the energy budget again, and we getw → 0 andε ϕ k ε ϕ 2 p → 1/2 at late times. b) 2 < p < 4: The daughter field gets excited through broad parametric resonance as long asq ≡ q * a 6(p−4) p+2 1. However, unlike the p = 2 case, now the inflaton also develops fluctuations via self-resonance. Although the excitation of the daughter field is initially much stronger than the one of the inflaton (see the typical Floquet indices in Fig. 1), it eventually terminates when the daughter field resonance becomes narrow (q 1), like in the p = 2 case. On the other hand, the strength of the inflaton self-resonance remains constant. Therefore, inflaton fluctuations are continuously produced even during the non-linear regime, and we end up in a universe withε ϕ k ,ε ϕ g → 1/2 andw → 1/3 at late times.
c) p ≥ 4: The resonance parameterq ≡ q * a 6(p−4) p+2 is either constant (for p = 4) or grows with time (for p > 4), so both inflaton self-resonance and broad parametric resonance of the daughter field are present at late times. In this case we end up in an equilibrium regime in which both fields end up completely fragmented and have the same energy: we getε ϕ k ,ε ϕ g , ε χ k ,ε χ g → 1/4 andw → 1/3 at late times.

Multiple daughter fields without (self-)interactions:
Let us now consider the case of an inflaton coupled to multiple daughter fields via quadraticquadratic interactions, described by potential (4.2). We consider the power-law coefficients: a) p ≥ 4, b) p = 2, and c) 2 < p < 4.
a) p ≥ 4: In Fig. 4 we show the post-inflationary evolution of the energy ratios for p = 4 and different number of daughter fields (N d = 2, 5, 10 and 20). Each daughter field is coupled Figure 4. [Potential (4.2), p = 4] Evolution of the energy ratios for different number of daughter fields: N d = 2, 5, 10 and 20. All of them are coupled to the inflaton with the same resonance parameter q * = 10 4 . We have plotted the sum of the kinetic and gradient contributions of the inflaton (blue) and the different daughter fields (red lines). The potential contributions, which become negligible at late times, are depicted in green and orange. The dashed line indicates the predicted value (4.5) for the energy ratios at late times. Figure 5.    to the inflaton with the same resonance parameter, q (n) * = 10 4 . For each field we depict the sum of its kinetic and gradient energy ratios: at late times this represents the total energy fraction stored in the field, as the potential energy contributions become negligible.
During the linear regime (for times u ∼ 0−100), the excitation strength for each daughter field is characterized by its resonance parameter q (n) * . However, although all daughter fields in Fig. 4 have the same q (n) * , we observe that the amount of energy transferred during this regime is larger for some daughter fields than for others. This effect can be attributed to the randomness of the initial fluctuations: varying the seed of the random generator changes which particular fields receive more energy.
However, in the deep non-linear stage we always get an equilibrium regime between the inflaton and the daughter fields, analogous to the one in the single daughter field case discussed in Sect. 4.1. More specifically, the energy is equally distributed at very late times between the inflaton and all daughter fields (as long as they are coupled in broad resonance), with the energy ratios satisfyinḡ For example, for N d = 2, 5, 10 and 20, each of the N d + 1 fields of the system get 33%, 17%, 9.1% and 4.8% of the total energy respectively. Remarkably, this equilibrium regime between the different daughter fields is achieved without including explicit interactions between them, such as X 2 i X 2 j . A similar equilibration regime is achieved even if the daughter fields are coupled to the inflaton with different strengths. This can be seen in the left panel of Fig. 5, where we depict the case N d = 2, q (1) * = 2000 and q (2) * = 1000. Both daughter fields end up with the same energy despite having different resonance parameters. However, the smaller the ratio q (2) * /q (1) * is, the later the equilibration between all fields is achieved. We illustrate this in the right panel of Fig. 5, where we depict the fraction of energy transferred to χ 2 as a function of time, for different ratios q (2) * /q (1) * . For small ratios q (2) * /q (1) * 1/10, we observe that the fraction of energy stored in χ 2 decreases for some time during the early non-linear regime, develops a local minimum and starts growing again. Eventually it attains the same energy as χ 1 at late times. In reality, for these ratios we are unable to observe the complete achievement of equilibration on the lattice, but the long-term trend can be extrapolated in all cases.
These examples illustrate that a significant depletion of the inflaton energy can be achieved in multi-field scenarios like the ones considered here, without relying on perturbative decay channels like in combined preheating scenarios [64][65][66][67].
It is also interesting to show how the equilibration regime between the different daughter fields is achieved in momentum space. For this purpose we inspect the power spectrum of the fields in natural variables, defined as f 2 ≡ d log κ P f (κ) for f = {ϕ, χ i }. In the left panel of Fig. 6 we show the time-evolution of the daughter field spectra for N d = 2 and p = 4 (i.e. the same case as in the left panel of Fig. 5). We observe that different ranges of momenta are populated during the linear regime for both fields, which are characterized by q (1) * and q (2) * for χ 1 and χ 2 respectively. However, the spectra of both fields converge at later times.
Finally, let us mention that we have simulated multi daughter field systems for values p > 4 (such as p = 4.5 or 5), and observed that a similar equilibration regime between the field emerges, i.e. the energy ratios at late times are also given by Eq. (4.5). Similarly to the single daughter field case, during the initial stage of broad parametric resonance, the energy ratios of all daughter fields increase exponentially, until backreaction effects trigger the decay of the inflaton homogeneous mode and the nonlinear regime starts. Remarkably, the amount of energy transferred to each daughter field during this stage is different even if they have the same q (n) * . As in the p = 4 case, this effect is induced by the initial random fluctuations, as different realizations change which particular fields get more energy. In any case, the broad resonance terminates onceq 1, so the inflaton homogeneous mode eventually recovers 100% of the total energy, and the energy transferred to the daughter field sector becomes subdominant. Similarly, the fluctuations of the inflaton dilute faster than its homogeneous mode, so the equation of state goes again tow →w hom = 0 at late times. Remarkably, we do not observe an equipartition regime between the different daughter fields during the non-linear regime: the field that has (randomly) received more energy during the linear stage will dominate the energy budget of the daughter field sector later on.
Although adding more daughter fields does not change the final energy distribution and equation of state for p = 2, it can have relevant effects at intermediate times. In order to illustrate this, we consider systems where the N d daughter fields have all the same q (n) * . In the left panel of Fig. 7 we show the maximum fraction of energy attained by the daughter field sector during the simulation as a function of N d and three different choices of q * (more specifically, we depict the sum n (ε χn k +ε χn g )| max ). For fixed q * , the larger the number of daughter fields the greater the energy transfer is (though the average energy transferred to each individual daughter field decreases for larger N d ).
As the energy transferred to field gradients increases with N d , so does the transitory behaviour of the equation of state. In the right panel of Fig. 8 we compare the evolution of the (effective) equation of state for N d = 1 and 10, for the same three choices of q * . As expected, the deviation fromw =w hom (= 0) towardsw =w max ( 1/3) is stronger in the N d = 10 case than in the N d = 1 one, although at late times we always recover the matter-dominated statew → 0. Figure 7.  Finally, we show the evolution of the daughter field spectra in the right panel of Fig. 6, for N d = 2, q (1) * = 10 5 and q (2) * = 10 4 . Initially, a narrow infrared band of modes is strongly amplified for each daughter field, characterized by q (1) * and q (2) * respectively. Once backreaction effects become relevant, the daughter fields populate a wider range of momenta. However, onceq (1) ,q (2) 1, the exchange of energy ceases and the spectra freeze. Therefore, the daughter field spectra do not converge and end up with different shapes, unlike in the p = 4 case. c) 2 < p < 4: The dynamics of the daughter field sector are, for these power-law coefficients, similar to the p = 2 case discussed above. In particular, as their resonance also becomes narrow at late times, the fraction of energy stored in the daughter field sector eventually becomes negligible,ε χ → 0. However, the inflaton now fragments due to its self-resonance, so we haveε ϕ k ,ε ϕ g → 1/2 andw → 1/3 at late times.
If σ * /q * = λ/g 2 1, the self-interaction can significantly affect the field dynamics during both the linear and non-linear regimes (this was already noted and investigated for the early preheating phase in [21]). During the linear regime, the self-interaction gives an effective mass to the daughter field and suppresses its resonant growth, as discussed in Sect. 3.1. However, during the non-linear regime the self-interaction also triggers a self-resonant excitation, as we shall see. We will consider the power-law coefficients: a) p = 2, a) 2 < p < 4, and c) p ≥ 4. a) p = 2: In Fig. 9 we show the evolution of the energy distribution for q * = 10 4 and two choices of σ * : σ * /q * = 0.1(<1) (top-left) and σ * /q * = 2 · 10 2 (>1) (top-right). In the case σ * /q * = 0.1, the effect of the self-interaction in the post-inflationary dynamics is negligible, so the energy distribution evolves in a very similar way as the σ * = 0 case discussed in Sect. 4.1. More specifically, the daughter field gets excited through broad parametric resonance, which triggers the decay of the inflaton homogeneous mode through backreaction effects at the time u ∼ 10 2 . On the other hand, for σ * /q * = 2 · 10 2 we observe that, although the energy of the daughter field grows initially through broad resonance as well, this growth is slowed down due to the effective mass (m eff χ ) 2 ∼σ χ 2 and never reaches a relevant magnitude before the resonance parameter falls belowq = 1 (indicated by the vertical dash-dotted line).
The energy distribution also evolves in different ways during the non-linear regime. In particular, in the case σ * /q * = 2 · 10 2 , the daughter field develops a relevant homogeneous mode during the early stage of broad resonance due to its quartic potential. This triggers a late growth of the daughter field fluctuations through a process of self-resonance, analogous to the one experienced by the inflaton for p > 2. This leads to an exponential growth ofε χ g at times u ∼ 4 · 10 2 − 10 3 , which can be observed in the top-right panel of Fig. 9. This effect can also be clearly observed in the spectral evolution of the fields, depicted in Fig. 10. At these times, the spectrum of the daughter field shows a distinct structure of narrow peaks, reminiscent of the narrow bands that appeared in the Floquet diagram of inflaton self-resonance (see left panel of Fig. 1). These peaks also get imprinted on the inflaton due to their interaction to the daughter field. Remarkably, backreaction effects lead to a wash out of these peaks in the daughter field spectrum, while they stay imprinted in the inflaton spectrum due to the lack of inflaton self-interactions.
Although the inflaton homogeneous mode recovers 100% of the energy at late times independently of the strength of the self-interaction, it can strongly affect the field dynamics at intermediate times. In order to illustrate this, in the bottom-left panel of Fig. 9 we have depicted the maximum value attained by the energy ratios during the simulation, for different ratios σ * /q * . The (temporary) energy transfer to the daughter field gets maximized for intermediate ratios σ * /q * ∼ 1 − 50, while it gets strongly suppressed for σ * /q * 1 due to the suppression of the resonance effects by the effective mass. Similarly, the self-interaction affects the evolution of the equation of state. Its qualitative evolution is similar to the σ * = 0 case discussed in Sect. 4.1: the production of field fluctuations triggers a transitory deviation fromw = 0 tow =w max < 1/3, which then relaxes back tow = 0 at late times. We illustrate this in the bottom-right panel of Fig. 9, which showsw max for different choices of σ * /q * . The maximum deviation towards radiation-domination takes place again for intermediate values Figure 9. [Potential (4.3), p = 2] Top panels: Evolution of energy ratios for q * = 10 4 and two different self-interaction strengths: σ * /q * = 0.1 (top-left) and σ * /q * = 2 · 10 2 (top-right). The vertical lines indicate whenq = 1 (dash-dotted) andσ = 1 (dashed). Bottom-left panel: Maximum value attained by the energy ratiosε ϕ g ,ε χ g ,ε χ k ,ε χ 4 p andε ϕ 2 χ 2 p , extracted from simulations with q * = 10 4 and different values of σ * . Bottom-right: Maximal value attained by the effective equation of state as a function of the ratio σ * /q * , for p = 2 and three different choices of q * . σ * /q * ∼ 1 − 50, while the deviation is minimal for very large ratios σ * /q * 1.
b) 2 < p < 4: The evolution of the daughter field energies is similar to the p = 2 case. However, the inflaton now fragments due to self-resonance, so we haveε ϕ k ε ϕ g → 0.5 at late times (and consequentlyε χ k ε χ g → 0). The system will eventually arrive at a radiation dominated state as well. However, for large enough σ * the transition phase can span over several e-folds, while it happens rather fast in the case of negligible self-interaction. c) p ≥ 4: In the left panel of Fig. 11 we show the fraction of energy stored in the daughter field, for p = 4, q * = 20, and different choices of σ * . For σ * /q * 1, the effect of the daughter field's quartic self-interaction in the post-inflationary dynamics is negligible, so the energy distribution evolves in a similar way as in the σ * = 0 case discussed in Sect. 4.1.
In particular, at late times the energy is equally distributed between the inflaton and the daughter field, according to Eq. (4.5).
On the other hand, for σ * /q * 1 the daughter field fluctuations experience a selfresonance process during the late-time regime, as described above for the p = 2 case. Due to this, the energy from the inflaton is transferred much more efficiently and the fraction of energy stored in the daughter field at late times can be larger than 50%. In fact, as seen in the right panel of Fig. 11, the larger the ratio σ * /q * , the larger the amount of transferred energy. We have observed the same behaviour for lattice simulations of the p = 4.5 case, and Figure 10.
we expect it to happen for all values of p > 4. These results show that even in the case of one daughter field, a significant amount of energy density can be extracted from the inflaton for power-law coefficients p ≥ 4, as long as a quartic self-interaction with λ g 2 is present in the theory.
For p > 4 the late excitation of fluctuations for large ratios σ * /q * , similar to the case seen in Fig. 11, has a relevant influence on the equation of state. While for σ * 0 the transition from the initial homogeneous to radiation dominated averaged equation of state happens rather fast, it appears for σ * q * as a smooth transition spanning over several e-folds.

Multiple daughter fields with (self-)interactions:
Finally, let us discuss the case of multiple daughter fields with interactions of the type λ nm X 2 n X 2 m , represented by potential (4.4). Due to the large number of free parameters, a detailed parametric analysis of the post-inflationary dynamics is not possible in this case. However, we will consider three particularizations of potential (4.4), that allow us to learn about the generic features of the post-inflationary energy distribution in this model. We fix the number of daughter fields to two (N d = 2) for simplicity, but our results can be easily generalized to N d > 2.
a) λ 11 > 0, λ 22 > 0, λ 12 = 0: We first consider a direct combination of the scenarios discussed in Sects. 4.2 and 4.3: an inflaton coupled to two daughter fields with quartic selfinteractions, 1, the role of both self-interactions is negligible and we recover the results of Sect. 4.2: the energy is equally distributed between the three fields at very late times if p ≥ 4, while the inflaton eventually recovers 100% of the energy for 2 ≤ p < 4.
On the other hand, if σ (i) * /q (i) * 1 for any of the daughter fields i = 1, 2, the energy transferred to that field during the non-linear regime gets enhanced due to the self-resonance induced by the quartic self-interaction (the same effect discussed in Sect. 4.3). If 2 ≤ p < 4, the inflaton will still recover 100% of the energy at very late times, as the effective parameter decreases with time. More interesting is the p ≥ 4 case: in this case the three fields no longer equilibrate with the same energy, and in particular the daughter fields with larger ratio σ (i) * /q (i) * will get a larger percentage of the energy at late times. The amount of energy that remains in the inflaton depends strongly on the particular ratios, but it is always less than 33% because more energy gets extracted from it.
We present two illustrative examples for p = 4 and N d = 2 in Fig. 12. The left panel shows a case in which σ (1) * /q (1) * = 0.01 (< 1) and σ (2) * /q (2) * = 2 (> 1), so only X 2 experiences significant self-resonance during the non-linear regime. We observe that the field X 2 gets more than 33% of the energy at very late times (approximately 38%), while both φ and X 1 equilibrate with less than that (each one gets 31% of the energy). On the other hand, the right panel shows a case in which σ (1) * /q (1) * = 2 (> 1) and σ (2) * /q (2) * = 3 (> 1), so both X 1 and X 2 experience self-resonance: in this case both daughter fields get more energy than φ at late times, and in fact X 2 gets more energy than X 1 due to its larger coupling ratio.
b) λ 11 = λ 22 = 0, λ 12 > 0: Let us now consider a scenario in which both daughter fields are coupled to each other through a quadratic-quadratic interaction, but neither of them Figure 13. [Potential (4.7), p = 4] Left: Evolution of the fraction of energy stored in the daughter fields (i.e. the sum ε χ1 k+g + ε χ2 k+g = ε χ1 k + ε χ1 g + ε χ2 k + ε χ2 g ), for q * = 20 and different ratios σ /q * at late times, when ε have quartic self-interactions. The potential reads as During the linear regime, the strength of the resonance and the momenta excited are characterized, for each of the two fields i = 1, 2, by the corresponding resonance parameter q (i) * ≡ g 2 i φ 2 * /ω 2 * . However, during the later non-linear regime, the quadratic-quadratic interaction λ 12 X 2 1 X 2 2 triggers a resonant process for both fields X 1 and X 2 , in a similar way as the quartic term λ i X 4 i did for the field i in the previous example. This allows for an efficient distribution of energy in the daughter field sector.
In order to illustrate this, let us consider first the case p = 4. Remarkably, in the simulations we have seen that both daughter fields equilibrate very quickly during the nonlinear regime (much faster than without such an interaction), and get the same fraction of energy at late times, i.e.ε χ 1 k +ε χ 1 g ≈ε χ 2 k +ε χ 2 g . This happens also when q (1) * = q (2) * , as well as when q (1) * , q (2) * σ (12) * . Due to this, in the left panel of Fig. 13 we have depicted the (at late times) total fraction of energy attained by the daughter field sector as a function of time (i.e. the sumε χ 1 k +ε χ 1 g +ε χ 2 k +ε χ 2 g ), for q . We observe that for σ (12) * /q * 1, the effect of the λ 12 X 2 1 X 2 2 interaction is negligible, so the daughter field sector gets 66% of the energy at late times, in agreement with Eq. (4.5). However, as σ (12) * /q * increases, the transfer of energy to the daughter field sector gets larger. For very large ratios it is difficult to simulate the system long enough to observe the achievement of the stationary regime. Thus, we have depicted in the right panel of Fig. 13 the summed energies ε ϕ k+g ≡ ε ϕ k + ε ϕ g and ε χ 1 A similar behavior has been observed in lattice simulations of scenarios with p > 4.
For p < 4, the inflaton recovers all the energy at late times as expected, while the energy stored in the daughter field sector becomes subdominant (for p = 2 we getε ϕ k ,ε ϕ 2 p → 1/2, and for 2 < p < 4 we getε ϕ k ,ε ϕ g → 1/2). This happens because the effective self-coupling parameter decreases with time asσ (12) p+2 , so the self-resonance triggered by it becomes eventually too weak. However, an interesting effect takes place in the daughter field sector for these power-law coefficients if q * σ (12) * : the daughter fields tend to equilibrate at late times, in contrast to the case of very small σ (12) * . In Fig. 14 we show the evolution of the spectra of both daughter fields for the cases: q (1) * = 7 · 10 3 , q (2) * = 3.5 · 10 5 and σ (1) * = 7 · 10 3 , q (2) * = 3.5 · 10 4 and σ (12) * = 7 · 10 5 (right panel). We can see that in the first case (for which q (2) * > σ (12) ) the two spectra do not equilibrate at late times (similar to the case discussed in Fig. 10), while in the second case the two daughter fields equilibrate, i.e. both spectra match at late times. c) g 2 1 > 0, g 2 2 = 0, λ 12 > 0, λ 11 = λ 22 = 0: Finally, let us consider a 'chain' scenario in which only one daughter field (X 1 ) is coupled directly to the inflaton. The second daughter field (X 2 ) is coupled to X 1 through a quadratic-quadratic interaction, which allows a transfer of energy from the inflaton to X 2 in a two-step process. The potential reads as There are two unspecified parameters to fix: q (1) * and σ (1) * = 10 (right panel). As expected, X 1 is excited during the initial linear stage through a process of broad parametric resonance in both cases. Later, we get a transitory equilibration phase during the early non-linear regime, in which the inflaton and X 1 each hold approximately 50% of the total energy (this happens for time scales u ∼ 10 2 − 10 3 in both cases). At these times, the energy of X 2 remains subdominant because it is not coupled directly to the inflaton, and hence it is not excited via parametric resonance. However, in a later stage X 1 starts to transfer energy into X 2 through the λ 12 X 2 1 X 2 2 interaction. This is reminiscent of the case studied in Ref. [68], where parametric resonance was induced by an inhomogeneous field. Eventually, we end in a situation in which X 1 and X 2 have equilibrated with the same energy, despite X 2 not being coupled directly to the inflaton. This happens for both ratios of σ (1) * we find that each of the three fields get 33% of the energy at late times, but if σ (12) * > q (1) * much more energy is eventually transferred to the daughter field sector. In the depicted case of σ (12) * /q (1) * = 10 only ∼ 18% of the energy remains in the inflaton. We observe the same qualitative behavior for values of p > 4. On the other hand, for p < 4 the energy ratios of both daughter fields attain a maximum whenq (1) ∼ 1 and then start decreasing again. In the later case, the second daughter field usually stays subdominant (except for very large values of q (1) * and σ (12) * ), as there is not enough time to enhance it significantly beforeσ 12 falls below unity.

Summary and discussion
The aim of this work has been to study the post-inflationary dynamics of the universe when the inflaton is coupled to multiple interacting daughter fields. As a set-up we have considered an inflaton that oscillates after inflation in a potential that is monomial around the minimum V (φ) ∝ |φ| p , and is coupled to daughter fields through scale-free, quadratic-quadratic interactions g 2 n φ 2 X 2 n . If q (n) * ≡ g 2 n φ 2 * /ω 2 * 1, the daughter fields get excited during the linear regime through a process of broad parametric resonance. We have also included additional scale-free interactions of the type λ nm X 2 n X 2 m , including also quartic self-interactions when n = m. We have first studied the excitation of the field fluctuations during the linear regime by means of a Hartree approximation, see Sect. 3. This analysis has allowed us to obtain a prediction for the maximum variance attained by the daughter field as a function of the particle couplings, see e.g. Fig. 2.
We have then simulated the post-inflationary dynamics of the system with (2+1)dimensional lattice simulations, see Sect. 4, which have allowed to capture the full dynamics from the end of inflation until the achievement of a stationary regime for a large set of model parameters. For clarity purposes, we have divided our results in four subsections 4.1 -4.4, which cover different particularizations of the most generic potential, see Eqs. (4.1)-(4.4). As a reference, we have used the results for one daughter field without quartic self-interaction.
This case was considered extensively in our Letter [55] and in Part I [56], and the results have been summarized in Sect. 4.1. There we observed that the amount of energy transferred to the daughter field after equilibration depends on p in the following way: i) for p = 2, 100% of the energy ends up in the homogeneous mode of the inflaton, ii) for 2 < p < 4, 100% of the energy ends up in a fragmented inflaton, and iii) for p ≥ 4, the inflaton and the daughter field end up fragmented with each field getting 50% of the energy. In this paper, we have investigated how these results change when two or more daughter fields are included in the theory, as well as scale-free λ nm X 2 n X 2 m and λ n X 4 n interactions: • In Sec. 4.2 we have considered the case of multiple daughter fields without X 2 m X 2 n interactions. For p ≤ 4, the results for the energy distribution at very late times remain unchanged with respect to the single daughter field case; however, a larger transfer of energy to the daughter field sector takes place at intermediate times. This has noticeable effects on the equation of state for the p = 2 case, see Fig. 8, and therefore on the predictions of the CMB observables n s and r. On the other hand, for p ≥ 4 the energy ends up equally distributed between all fields at very late times, see Eq. (4.5) and Fig. 4. This is true even if the daughter fields are coupled to the inflaton with different strengths, as long as they are excited in broad resonance, see Fig. 5. Remarkably, this means that in systems with many daughter fields, a significant amount of energy can be depleted from the inflaton, without relying on perturbative decay channels or other elements to the theory.
• In Sec. 4.3 we have explored how the results for one daughter field change when a selfinteraction λX 4 is included in the theory. Results are basically unchanged if λ/g 2 1. However, if λ/g 2 1, the self-interaction does significantly affect the evolution of the energy distribution after inflation. First, during the linear regime it appears as an effective mass to the daughter field fluctuations, which suppresses the excitation process. However, during the non-linear regime it triggers a self-resonance process, which enhances the amount of energy transferred to the daughter field. Due to this, for p ≥ 4 the daughter field gets more than 50% of the energy at late times, see Fig. 11. In fact, the larger the ratio λ/g 2 , the more energy is transferred. Correspondingly, the inflaton retains less than 50% of the total energy density. For p < 4, the inflaton recovers 100% of the energy at very late times as in the λ = 0 case, but a larger amount of energy gets temporally transferred at intermediate times during the early non-linear stage, see Fig. 9. For all power-law coefficients (except p = 4), a large enough ratio λ/g 2 significantly affects the evolution of the equation of state, see e.g. the bottom right panel of Fig. 9 for p = 2. For 2 < p < 4 and p > 4, the large ratio also gives rise to a prolonged transition to radiation domination. These changes also affect the predictions for the CMB observables n s and r.
• Finally, in Sec. 4.4 we have considered different scenarios involving multiple daughter fields and scale-free interactions λ nm X 2 n X 2 m and λ n X 4 n . We have observed that a quadratic-quadratic interaction between two daughter fields, such as e.g. λ 12 X 2 1 X 2 2 , allows for an efficient exchange of energy between the X 1 and X 2 fields, and if λ 12 g 2 1 , g 2 2 , they can indeed achieve an equilibration regime with the same energy at late times, see Fig. 14. Moreover, this term can also induce a resonance during the non-linear regime to both fields X 1 and X 2 , similar to the one induced by the λ 1 X 4 1 interaction to X 1 . For p ≥ 4, this interaction can increase the amount of energy transferred from the inflaton to the daughter fields at late times, see Fig. 13. In fact, even when a given daughter field (say X 2 ) does not have a direct interaction to the inflaton, the resonance induced by the λ 12 X 2 1 X 2 2 interaction does indeed allow for an efficient transfer of energy to X 2 by means of a two-step process: first there is an energy transfer from φ to X 1 through the φ 2 X 2 1 interaction, and later on from X 1 to X 2 through the X 2 1 X 2 2 interaction. An example of this is shown in Fig. 15.
Furthermore we like to note that in our set-up, a significant amount of energy remains on the inflaton at late times (one exception is the hypothetical scenario p ≥ 4 and N d → ∞, where all energy is depleted from the inflaton). The energy must be transferred to light fields via some other mechanisms, such as perturbative decay channels. A similar mechanism is necessary to achieve a radiation-dominated stage in the p = 2 case. Small extensions of our scenario are required for this purpose, such as additional couplings that only become relevant at a later stage. However, let us emphasize that for p > 2, the transition towards a radiation-dominated state is completed after few e-folds (which we can capture with lattice simulations), and the equation of state will not change anymore until BBN.
In this work we have assumed that during the post-inflationary stage of inflaton oscillations, the inflaton potential can be approximated by a simple monomial function V (φ) = |φ| p with p ≥ 2 around the minimum, and that only quadratic-quadratic interactions exist between the different fields. However, it would be very interesting to explore multi-daughter field theories beyond these two assumptions. For example, one could consider lower-scale inflaton potentials in which the minimum of the potential can be expanded as V (φ) = |φ| p around the minimum, but the inflaton does indeed oscillate over flatter-than-quadratic regions (e.g. the case of potential (2.2) with M m pl ). Furthermore, studying the effect of trilinear interactions between the different fields in the post-inflationary energy distribution and equation of state would also be very interesting. One could incorporate terms of the type φX 2 n like in Ref. [44], which induce a stage of tachyonic resonance during the linear regime. One could also incorporate trilinear interactions between the different daughter fields such as X 2 m X n or X m X n X p . However, unlike quadratic-quadratic interactions, these terms generate additional mass scales into the theory which could be difficult to properly capture on the lattice.
Another interesting extension of our work could be to study the evolution of the metric perturbations in set-ups with multiple daughter fields. It has been shown that the postinflationary oscillations of the homogeneous inflaton can trigger a resonant growth of metric perturbations at sub-Hubble scales, a process known as 'metric preheating' [69][70][71][72][73][74]. This may lead to interesting phenomenology such as production of black holes or gravitational waves. Metric preheating has been recently studied in the presence of interactions of the inflaton to cosmological fluids in the form of perturbative decays [75] (see also [76] for an study in multi-field inflation). In the scenario considered in our work, the inflaton is coupled to (one or multiple) daughter fields through quadratic-quadratic interactions, which cause its fragmentation few e-folds after the end of inflation (for p > 2, inflaton fragmentation happens via self-resonance even in the absence of such interactions). Therefore, in order to fully understand of the fate of metric preheating in our set-up, one would need to study the dynamics of the system with numerical relativity simulations, which allow for metric perturbations.
Lattice studies have shown that this kind of field systems virialize very quickly [35,36,41,77], with the fields obeying the following 'equipartition' identities when averaged over both volume and oscillations, These identities can be evaluated for the multi field model under consideration, described by potential (2.7). If we take the monomial approximation (2.1) for the inflaton potential, these can be expressed in terms of energy contributions as follows, Typically, the sum E tot ≡ i E i does not change significantly over one oscillation, thus one can also express these identities in terms of the (oscillation-averaged) energy density ratios by replacing E i Tφ →ε i .