Tachyonic production of dark relics : classical lattice vs. quantum 2PI in Hartree truncation

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Tachyonic production of dark relics : classical lattice vs. quantum 2PI in Hartree truncation © The Authors. Article funded by SCOAP3 Published version Kainulainen, Kimmo; Nurmi, Sami; Väisänen, Olli Kainulainen, K., Nurmi, S., & Väisänen, O. (2024). Tachyonic production of dark relics : classical lattice vs. quantum 2PI in Hartree truncation. Journal of High Energy Physics, 2024, Article 9. https://doi.org/10.1007/jhep10(2024)009 2024 JHEP10(2024)009 Published for SISSA by Springer Received: July 2, 2024 Accepted: August 29, 2024 Published: October 1, 2024 Tachyonic production of dark relics: classical lattice vs. quantum 2PI in Hartree truncation Kimmo Kainulainen,a,b Sami Nurmia,b and Olli Väisänen a aDepartment of Physics, University of Jyväskylä, PL 35 (YFL), Jyväskylä 40014, Finland bHelsinki Institute of Physics, University of Helsinki, PL 64, Helsinki 00014, Finland E-mail: [email protected],[email protected], [email protected] Abstract: We study the out-of-equilibrium production of non-minimally coupled selfinteracting scalar dark matter during reheating using classical lattice simulations. The outcomes of the classical simulations are in qualitative agreement with the previous results obtained using the quantum 2PI approach in the Hartree truncation. In particular, the novel non-linear resonance found in the 2PI Hartee study is present also in the classical lattice simulations and can dominate the final dark matter yield. For the parameters considered, the difference in final value of the scalar two-point function between the two approaches is a factor of O (1). Keywords: Early Universe Particle Physics, Models for Dark Matter, Non-Equilibrium Field Theory, Nonperturbative Effects ArXiv ePrint: 2406.17468 Open Access,©The Authors. Article funded by SCOAP3.https://doi.org/10.1007/JHEP10(2024)009 JHEP10(2024)009 Contents 1 Introduction 1 2 The setup 2 3 The lattice implementation 3 4 Lattice results and comparison with the 2PI Hatree approach 5 5 Conclusions and discussion 12 1 Introduction Out-of-equilibrium dynamics of quantum matter in the presence of classical scalar fields is involved in many primordial processes. Important examples include resonant phenomena [ 1 – 5 ] and tachyonic instabilities [ 6 – 12 ] in reheating, the electroweak baryogenesis [ 13 – 19 ], the leptogenesis mechanism [ 20 – 29 ], and also many dark matter setups with very weakly coupled fields [ 30 – 37 ]. Resolving the strongly non-linear dynamics encountered in such setups often requires non-perturbative application of field theory methods [ 38 – 43 ]. In this work we study non-linear dynamics in the dark matter scenario proposed in [ 35 , 36 ], where a spectator scalar singlet with a non-minimal coupling ξRχ2 undergoes tachyonic instability during reheating, when the classical Ricci scalar R oscillates from positive to negative values. The singlet has no non-gravitational couplings to the visible sector and the produced excitations constitute a dark matter component [ 35 , 36 ]. In [ 44 ] the quantum dynamics in the setup was investigated using the 2PI approach [ 45 , 46 ] in the Hartree truncation. The results of [ 44 ] indicate that the tachyonic instability is followed by a novel transient resonant stage driven by the two-point function ⟨χ2⟩ in the presence of the self-coupling λχ4 . The resonance can enhance the net particle production by an order of magnitude compared to the semi-analytical estimates of [ 36 ], based on [ 12 ]. The results of [ 44 ] indicate that the transient resonance occurs somewhat after the point when the effective mass contribution ξR falls below λ⟨χ2⟩ , and its strength depends sensitively on the couplings λ and ξ , and on the equation of state of the universe during reheating. In the current work, we run classical lattice simulations in the setup of [ 35 , 36 ] and make a quantitative comparison with the 2PI Hartee results of [ 44 ]. There are two conceptually different effects that can generate differences between the two approaches. First, while [ 44 ] investigates quantum evolution starting from vacuum initial conditions, the classical lattice simulation describes on-shell dynamics starting from an initial field configuration with classical plane waves. Second, the Hartree truncation used in [ 44 ] does not account for momentum exchanging processes and it is a priori unclear to what extent they could affect the resonant growth. All on-shell effects of the mode mixing induced by tree level-lagrangian are included in the classical lattice simulation. A third aspect is that the sources of numerical error in the lattice simulation and in solving the 2PI equations as implemented in [ 44 ] are different and it is not a priori obvious which approach is more efficient in this respect. – 1 – JHEP10(2024)009 We perform the lattice simulations using the CosmoLattice code [ 47 , 48 ]. CosmoLattice and other lattice codes have already been applied to similar non-minimally coupled scalar setups in [ 49 – 52 ], but for a reheating equation of state and spectator couplings different from those used in [ 44 ], and direct comparison of the results is therefore not possible. The results obtained in this work confirm that the transient resonance seen in [ 44 ] is present also in the classical lattice solution and the net particle production is in a broad agreement with [ 44 ]. Details of the resonant stage and the final momentum distribution of the two-point function measured from the classical lattice solutions however differ from the corresponding results in [ 44 ]. This paper is organized as follows. In section 2we briefly specify the setup, in section 3 we discuss relevant aspects of the lattice formulation, and in section 4we present the lattice results and compare them against [ 44 ]. Lastly, section 5presents our conclusions. 2 The setup We study the setup of [ 36 ] with a non-minimally coupled spectator scalar singlet χ , whose action is given by Sχ=Zd4x√−g1 2(∇µχ)(∇µχ)−1 2m2 χχ2+ξ 2Rχ2−λ 4χ4.(2.1) Following [ 44 ], we use the (+ ,−,−,− )sign convention. We study the dynamics during the reheating stage after inflation, assuming the universe is dominated by a homogeneous inflaton field oscillating in a quadratic potential. The classical equations of motion for the setup read ¨χ+ 3H˙χ−1 a2∇2χ+ (ξR +m2 χ)χ+λχ3= 0 ,(2.2) ¨ ϕ+ 3Hϕ +m2 ϕ= 0 ,(2.3) H=1 √6MP˙ ϕ2+m2 ϕϕ21/2,(2.4) and the Ricci scalar is given by R=1 M2 P˙ ϕ2−2m2 ϕϕ2.(2.5) In [ 44 ], the inflaton decay was also accounted for by including an effective decay term Γ ˙ ϕ and a radiation fluid in the analysis. The resonant growth of ⟨χ2⟩ was found to be strongest for Γ=0. Here we will set Γ=0in order to compare the 2PI results of [ 44 ] against classical lattice simulations in the non-trivial limit where the resonance is expected to have a maximal impact on the dynamics. In reality, the inflaton of course eventually needs to decay into radiation and the Γ = 0 case should therefore be physically understood as the limit where the inflaton decay can be neglected in the time-scale of our simulation runs. Throughout this work we focus on the parameter range, where the singlet χ is an energetically subdominant spectator. The energy density of χ with the non-minimal coupling (2.1) is given by ρχ=1 2˙χ2+1 2a2|∇χ|2+1 2m2 χχ2+λ 4χ4+3ξH2χ2+2Hχ ˙χ−2 3a2∇·(χ∇χ),(2.6) – 2 – JHEP10(2024)009 where the brackets denote volume averages in the case of a single classical lattice simulation and ensemble averages in the quantum computation, (for which the non-commuting terms also need to be symmetrised). The last total derivative term is irrelevant in our case, as it vanishes both in a spatially isotropic quantum system and in a classical setup with periodic boundary conditions. The condition for χ being a spectator is that ρχ≪ 3 H2M2 P throughout the computation so that its contribution to the Friedmann equation (2.4) can be neglected. 3 The lattice implementation We perform classical lattice simulation of the spectator field χ during reheating using a modified version of the CosmoLattice code [ 47 , 48 ]. We first solve separately the homogeneous equations (2.3) and (2.4) for the inflaton and the scale factor. The built-in function used by CosmoLattice for a ( t )in the case of a fixed background was then replaced by an interpolant of this solution. The spectator field is evolved in this background using CosmoLattice’s leapfrog evolver with the non-minimal coupling to the Ricci scalar in (2.2) added to the evolution kernel. For a comparison, selected runs were also performed using a fourth-order velocity-verlet evolver, whose output agreed with the leapfrog solutions. The spectra and most averaged quantities of the spectator field solution were extracted using built-in CosmoLattice routines, but as in [ 49 ] the measured field energy was modified to include the corrections from the non-minimal coupling in equation (2.6). The evolution is performed in comoving direct space on a cubic grid of N points per dimension. The lattice parameters were chosen so that the comoving infrared cutoff corresponds to the initial Hubble scale, kIR = Hi (see the text below for the definition of the initial time ti ). The lattice size should be chosen so that the linear ultraviolet cutoff klin UV = NkIR/ 2is well above the sharp drop-offs seen in the spectra presented in [ 44 ]. We find that this is achieved with N = 512, though larger lattices were used to test the solutions for resolution dependence. The maximum momentum the lattice can contain is kUV = √3klin UV . The system is evolved in conformal time with the timestep dτ = 10 −4H−1 i and the calculation is performed in units of the initial Hubble scale H−1 i . Initial conditions. Following [ 44 ], we initialize the homogeneous inflaton sector with slow roll initial conditions at ϕ = 15 MP . Inflation ends at ϵH = −˙ H/H2 = 1, corresponding to a conformal time τ0 . We start the lattice computation at the initial time τi defined slightly before the end of inflation as the moment when Ni = ln ( a0/ai ) = 2 . 485. In [ 44 ] the quantum field χ was initialized in the Bunch Davies vacuum with a vanishing one-point function ⟨χ ( τi, k) ⟩ = 0 and with the spectrum of the two point function, ⟨χ ( τi, k) χ ( τi,k′ ) ⟩ = (2 π ) 3δ (k+ k′ ) Pi ( k ), given by Pi(k) = π 4a2 i e−πIm(ν)(−τi)|H(1) ν(−kτi)|2.(3.1) Here τi = − ( aiHi ) −1 and H(1) ν is the Hankel function of the first kind with the index ν2 = 1 4− 12( ξ−1 6 ). Here we study the limit ξ≫ 1and Hi≫mχ where the effective potential for χ is initially dominated by the non-minimal coupling. Therefore the field is initially effectively massive and the index ν is imaginary. For ξ = 50 used in this work and in [ 44 ], the modulus squared of the Hankel function is approximated to three digit precision by – 3 – JHEP10(2024)009 |H(1) ν ( x ) |2≈eπIm(ν) (2 /π ) /px2+|ν|2 for all argument values. Therefore, the initial spectrum can simply be approximated by the Minkowski result Pi(k) = 1 2Hia2 i 1 q(k/Hi)2+a2 i|ν|2,(3.2) which is supported in CosmoLattice by default. When the corresponding classical system is simulated on a lattice, the initial field configuration χ ( τi, x)is drawn from a distribution which generates the same one- and twopoint functions with ensemble averages replaced by volume averages over the lattice. We neglect all higher order connected correlators at the initial time as the non-minimal coupling ξRχ2 initially dominates the effective potential. The initial field χ ( τi, x)therefore has a Gaussian distribution with zero mean and a spectrum determined by the discretised version of (3.2). We note that CosmoLattice discretises a continuum initial spectrum of the form (3.2) by substituting k2 with the square of kj = njkIR , where nj labels the lattice sites and kIR is the smallest comoving momentum on the lattice [ 47 ]. This does not fully coincide with the actual vacuum spectrum on lattice, (see equation (3.4) below), obtained by using the lattice dispersion relation that follows from discretising the derivative operator. However, we have checked that in our setup the spectrum rapidly evolves to (3.4) well before the onset of tachyonic instability. Subtracting the vacuum. The results for the two-point function in [ 44 ] were given in terms of δ ∆ F = ∆ F− ∆ F0 , where ∆ F is the finite part of the full two-point function and ∆ F0 the finite part of the vacuum two-point function, computed for the effective mass solved from the 2PI gap equation and excluding the tachyonic modes. To compare with these results, we need to perform a similar subtraction of the vacuum contribution from the two-point function measured on the lattice. The lattice vacuum power spectrum can be obtained from equation (3.2) by substituting k2 with k2 lat =4L2 π2X i sin2 πki 2L!,(3.3) which is the Fourier space representation of −∇2 for the symmetric definition of the lattice derivative operator [ 47 ]. Here ki are the Cartesian components of the comoving momentum and L = NkIR/ 2is half the comoving length of the cubic lattice box side. In addition, we must restrict (3.2) to only include modes which fit on the lattice and exclude eventual tachyonic modes as was done in [ 44 ]. The resulting angle-averaged vacuum power spectrum is Pvac(k) = ZdΩ 4π 1 H 1 2a2 Θ(k2 lat +M2 eff) qk2 lat/H2+M2 eff/H2Y i Θ(L−|ki|),(3.4) where Θdenotes a step function. The first step function cuts out tachyonic modes when M2 eff < 0and the latter product of three step functions ensures that only modes that fit on the cubic lattice are included. – 4 – JHEP10(2024)009 We define the effective mass M2 eff as the solution of M2 eff =a2m2 χ−a2ξ−1 6R+ 3λa2⟨χ⟩2+ 3λa2⟨χ2⟩−⟨χ2⟩vac(3.5) +3λ 16π2"M2 eff ln M2 eff a2m2 χ!−M2 eff +a2m2 χ#, where ⟨χ2⟩ is the contact limit of the full two-point function measured from the lattice and ⟨χ2 vac⟩ = R d 3k/ (2 π ) 3Pvac ( k ), with Pvac given by equation (3.4). Equation (3.5), which can be solved iteratively alongside equation (3.4), is just the 2PI gap equation of [ 44 ] written in terms of lattice quantities and tree-level couplings. The gap equation is introduced here only for the purpose of comparing our results with [ 44 ], as to this end we need to define M2 eff and subtract the vacuum part from the lattice two-point function in same manner as was done in [ 44 ]. In particular, the gap equation does not enter the classical lattice equations of motion in any way. However, as we will discuss below, the vacuum spectrum (3.4) with M2 eff solved from (3.5) matches well with the lattice two-point function before the tachyonic or resonant particle production starts. Finally, following [ 44 ] we split M2 eff into different components as M2 R=−a2ξ−1 6R, (3.6) M2 ∆= 3λa2⟨χ2⟩−⟨χ2⟩vac,(3.7) M2 σ=M2 eff −M2 R−M2 ∆.(3.8) We stress that while the effective mass defined by (3.5) allows a sensible definition of the vacuum in the lattice calculation, it differs essentially from its 2PI-counterpart, which, in the 2PI-method, controls the dynamical evolution of the 2-point correlation function and hence of the variance, including the back-reaction from one to another. 4 Lattice results and comparison with the 2PI Hatree approach We ran the lattice simulation for three setups matching those studied in [ 44 ]. The inflaton mass and the initial inflaton field value were chosen to be mϕ = 1 . 5 × 10 −13 GeV and ϕ = 15 MP , and the slow roll initial conditions were used for ˙ ϕ . The lattice simulations were initialized at Ni = 2 . 485 e-folds before the end of inflation (corresponding to the inflaton value ϕ = √2MP ) and we set ai = 1. In the following we denote by a0≈ 12 . 00 the scale factor at the end of inflation defined by ϵH = −˙ H/H2 = 1. We set ξ = 50 in all of our simulations and vary the self-coupling λ . Unless noted otherwise, the runs were performed with N = 512 and a comoving computation box size of H−1 i . Variances, effective masses and energies. The left panel of figure 1shows the lattice results for the variances ⟨χ2⟩−⟨χ2⟩vac , where the subtracted vacuum part is determined by equation (3.4). The right panel shows the corresponding squared effective masses solved using the gap equation (3.5). For comparison, we have also plotted the corresponding 2PI Hartee results of [ 44 ] with transparent lines. The lattice results for the effective mass components, defined in equation (3.6), are shown in figure 2. – 5 – JHEP10(2024)009 100101 a / a 0 100 102 104 106 108 a 2(2 2 vac )/ m 2 = 10 1 = 10 4 = 10 7 100101 a / a 0 102 103 104 105 M 2 eff / m 2 = 10 1 = 10 4 = 10 7 Figure 1. (Left panel.) The measured lattice variances (solid lines) and the corresponding 2PI solutions reproduced from [ 44 ] (transparent lines) for runs with λ = 10 −1 , λ = 10 −4 and λ = 10 −7 . The vacuum has been subtracted off from all lattice variances, and the 2PI-solution for λ = 10 −7 overlaps with the lattice results. (Right panel.) The lattice results for the effective mass calculated using equation (3.5) (solid lines) compared with their 2PI equivalents (transparent lines). The initial growth of the variance is driven by tachyonic periods where the effective mass becomes imaginary. The end of the tachyonic growth depends on the self-coupling: larger values of λ correspond to earlier onset of non-linear dynamics which shuts off the tachyonic growth. For λ = 10 −7 the non-linear region is never reached. Instead, the growth ends as the effective mass decreases and the tachyonic windows become increasingly narrow, as discussed in [ 44 ]. For λ = 10 −4 and λ = 10 −1 , the tachyonic stage ends when the effective mass contribution from the two-point function, M2 ∆ becomes comparable to M2 R . After this the solutions experience a transient resonance during which the lattice results for ⟨χ2⟩−⟨χ2⟩vac grow by factors of 5 and 10, respectively. For λ= 10−1there is a short period of non-linear oscillations between the end of the tachyonic stage and the onset of the resonance. After the end of the resonance, the variances redshift close to a law a−2 . Beyond a/a0≈ 30, our lattice results start to develop numerical resolution related spurious effects which we discuss at the end of this section. The lattice results for the variances and the effective masses closely agree with [ 44 ] in the region where the dynamics is essentially linear: for λ = 10 −7 this is the case throughout the simulation, and for λ = 10 −4 and λ = 10 −1 until the end of the tachyonic growth. The subsequent non-linear evolution for λ = 10 −4 and λ = 10 −1 also share qualitatively similar features between the classical lattice and the quantum 2PI Hartree solutions. In particular, the lattice results confirm the existence of the non-linear transient resonance which was first observed in [ 44 ] and which, within the 2PI formalism, is clearly seen to be driven by the two-point function ⟨χ2⟩ . There are however diffences in the details of the non-linear dynamics. – 6 – JHEP10(2024)009 100101 a / a 0 10 6 10 4 10 2 100 102 104 M 2 i / m 2 = 10 7 M 2 eff M 2 R M 2 M 2 100101 a / a 0 10 3 10 2 10 1 100 101 102 103 104 105 M 2 i / m 2 = 10 4 M 2 eff M 2 R M 2 M 2 100101 a / a 0 101 102 103 104 105 M 2 i / m 2 = 10 1 M 2 eff M 2 R M 2 M 2 Figure 2. The components of the effective mass for runs with λ = 10 −1 , λ = 10 −4 and λ = 10 −7 as defined in equations (3.6)–(3.8). In the lattice results for λ = 10 −4 the resonant-like growth starts almost immediately after the tachyonic stage when the effective mass contribution M2 ∆ quickly overtakes M2 R . This is contrary to [ 44 ], where one sees a plateau of coexistence between M2 R and M2 ∆ at a/a0∼ 5. However, such a plateau is present in the lattice solution for λ = 10 −1 , which is qualitatively very similar to the 2PI runs until the onset of the resonance. In both cases the non-linear growth is more efficient on the lattice and the two-point function saturates to slightly larger final values than in [ 44 ]. It can also be noted that oscillations of the two-point function decay faster on lattice at the end of the resonance than in the solutions of [ 44 ]. We expect the differences arise mostly from the lowest order Hartree truncation used in the 2PI study of [ 44 ]. In particular, the Hartree truncation neglects momentum exchanging scatterings while the lattice simulation contains the full non-linear dynamics at the classical level. Hartree results therefore do not properly capture processes that smear out coherent oscillations through the self coupling. On this basis, one could have expected that the lattice results would show less prominent resonant growth than the Hartree results. However, as discussed above, and seen in figure 1, the situation is actually the opposite at least for the parameter choices studied in this work. Another fundamental difference is that while the classical lattice simulation of course describes only on-shell physics, the 2PI results of [ 44 ] capture, within the Hartree truncation, the full quantum evolution of the system. Also, the dynamical coupling between the effective gap-mass and the 2-point function and hence the variance in the 2PI-approach may play an important role in explaining the difference. 1 We will return to this in more detail in a future work where we plan to extend the 2PI analysis 1 Indeed, we have checked that replacing the effective mass in 2PI-equations with the interpolated effective mass from a corresponding lattice run leads to an inconsistent 2PI-evolution. – 7 – JHEP10(2024)009 [4] J. Braden, L. Kofman and N. Barnaby, Reheating the Universe After Multi-Field Inflation, JCAP 07 (2010) 016 [arXiv:1005.2196] [INSPIRE]. [5] J. Berges and J. Serreau, Parametric resonance in quantum field theory,Phys. Rev. Lett. 91 (2003) 111601 [hep-ph/0208070] [INSPIRE]. [6] E. Calzetta, Spinodal Decomposition in Quantum Field Theory,Annals Phys. 190 (1989) 32 [INSPIRE]. [7] A.H. Guth and S.-Y. Pi, The Quantum Mechanics of the Scalar Field in the New Inflationary Universe,Phys. Rev. D 32 (1985) 1899 [INSPIRE]. [8] E.J. Weinberg and A.-Q. Wu, Understanding Complex Perturbative Effective Potentials,Phys. Rev. D 36 (1987) 2474 [INSPIRE]. [9] B.A. Bassett and S. Liberati, Geometric reheating after inflation,Phys. Rev. D 58 (1998) 021302 [Erratum ibid. 60 (1999) 049902] [hep-ph/9709417] [INSPIRE]. [10] G.N. Felder et al., Dynamics of symmetry breaking and tachyonic preheating,Phys. Rev. Lett. 87 (2001) 011601 [hep-ph/0012142] [INSPIRE]. [11] G.N. Felder, L. Kofman and A.D. Linde, Tachyonic instability and dynamics of spontaneous symmetry breaking,Phys. Rev. D 64 (2001) 123517 [hep-th/0106179] [INSPIRE]. [12] J.F. Dufaux et al., Preheating with trilinear interactions: Tachyonic resonance,JCAP 07 (2006) 006 [hep-ph/0602144] [INSPIRE]. [13] J.M. Cline, M. Joyce and K. Kainulainen, Supersymmetric electroweak baryogenesis,JHEP 07 (2000) 018 [hep-ph/0006119] [INSPIRE]. [14] K. Kainulainen, T. Prokopec, M.G. Schmidt and S. Weinstock, First principle derivation of semiclassical force for electroweak baryogenesis,JHEP 06 (2001) 031 [hep-ph/0105295] [INSPIRE]. [15] K. Kainulainen, T. Prokopec, M.G. Schmidt and S. Weinstock, Semiclassical force for electroweak baryogenesis: Three-dimensional derivation,Phys. Rev. D 66 (2002) 043502 [hep-ph/0202177] [INSPIRE]. [16] J.M. Cline, K. Kainulainen, P. Scott and C. Weniger, Update on scalar singlet dark matter,Phys. Rev. D 88 (2013) 055025 [Erratum ibid. 92 (2015) 039906] [arXiv:1306.4710] [INSPIRE]. [17] J.M. Cline and K. Kainulainen, Electroweak baryogenesis at high bubble wall velocities,Phys. Rev. D 101 (2020) 063525 [arXiv:2001.00568] [INSPIRE]. [18] T. Konstandin, Quantum Transport and Electroweak Baryogenesis,Phys. Usp. 56 (2013) 747 [arXiv:1302.6713] [INSPIRE]. [19] K. Kainulainen, CP-violating transport theory for electroweak baryogenesis with thermal corrections,JCAP 11 (2021) 042 [arXiv:2108.08336] [INSPIRE]. [20] W. Buchmüller and S. Fredenhagen, Quantum mechanics of baryogenesis,Phys. Lett. B 483 (2000) 217 [hep-ph/0004145] [INSPIRE]. [21] M. Beneke et al., Flavoured Leptogenesis in the CTP Formalism,Nucl. Phys. B 843 (2011) 177 [arXiv:1007.4783] [INSPIRE]. [22] A. Anisimov, W. Buchmüller, M. Drewes and S. Mendizabal, Quantum Leptogenesis I,Annals Phys. 326 (2011) 1998 [arXiv:1012.5821] [INSPIRE]. [23] P.S.B. Dev et al., Flavor effects in leptogenesis,Int. J. Mod. Phys. A 33 (2018) 1842001 [arXiv:1711.02861] [INSPIRE]. – 14 – JHEP10(2024)009 [24] A. De Simone and A. Riotto, Quantum Boltzmann Equations and Leptogenesis,JCAP 08 (2007) 002 [hep-ph/0703175] [INSPIRE]. [25] M. Garny, A. Hohenegger, A. Kartavtsev and M. Lindner, Systematic approach to leptogenesis in nonequilibrium QFT: Self-energy contribution to the CP-violating parameter,Phys. Rev. D 81 (2010) 085027 [arXiv:0911.4122] [INSPIRE]. [26] B. Garbrecht and M. Herranen, Effective Theory of Resonant Leptogenesis in the Closed-Time-Path Approach,Nucl. Phys. B 861 (2012) 17 [arXiv:1112.5954] [INSPIRE]. [27] M. Garny, A. Kartavtsev and A. Hohenegger, Leptogenesis from first principles in the resonant regime,Annals Phys. 328 (2013) 26 [arXiv:1112.6428] [INSPIRE]. [28] B. Dev et al., Resonant enhancement in leptogenesis,Int. J. Mod. Phys. A 33 (2018) 1842003 [arXiv:1711.02863] [INSPIRE]. [29] H. Jukkala, K. Kainulainen and P.M. Rahkila, Flavour mixing transport theory and resonant leptogenesis,JHEP 09 (2021) 119 [arXiv:2104.03998] [INSPIRE]. [30] E.W. Kolb, D.J.H. Chung and A. Riotto, WIMPzillas!,AIP Conf. Proc. 484 (1999) 91 [hep-ph/9810361] [INSPIRE]. [31] D.J.H. Chung, E.W. Kolb and A. Riotto, Superheavy dark matter,Phys. Rev. D 59 (1998) 023501 [hep-ph/9802238] [INSPIRE]. [32] M. Garny, M.C. Sandora and M.S. Sloth, Planckian Interacting Massive Particles as Dark Matter,Phys. Rev. Lett. 116 (2016) 101302 [arXiv:1511.03278] [INSPIRE]. [33] M. Garny, A. Palessandro, M.C. Sandora and M.S. Sloth, Theory and Phenomenology of Planckian Interacting Massive Particles as Dark Matter,JCAP 02 (2018) 027 [arXiv:1709.09688] [INSPIRE]. [34] Y. Tang and Y.-L. Wu, Pure Gravitational Dark Matter, Its Mass and Signatures,Phys. Lett. B 758 (2016) 402 [arXiv:1604.04701] [INSPIRE]. [35] T. Markkanen and S. Nurmi, Dark matter from gravitational particle production at reheating, JCAP 02 (2017) 008 [arXiv:1512.07288] [INSPIRE]. [36] M. Fairbairn, K. Kainulainen, T. Markkanen and S. Nurmi, Despicable Dark Relics: generated by gravity with unconstrained masses,JCAP 04 (2019) 005 [arXiv:1808.08236] [INSPIRE]. [37] J.A.R. Cembranos, L.J. Garay and J.M. Sánchez Velázquez, Gravitational production of scalar dark matter,JHEP 06 (2020) 084 [arXiv:1910.13937] [INSPIRE]. [38] D. Boyanovsky and H.J. de Vega, Quantum rolling down out-of-equilibrium,Phys. Rev. D 47 (1993) 2343 [hep-th/9211044] [INSPIRE]. [39] D. Boyanovsky, D.-S. Lee and A. Singh, Phase transitions out-of-equilibrium: Domain formation and growth,Phys. Rev. D 48 (1993) 800 [hep-th/9212083] [INSPIRE]. [40] J. Baacke and S. Michalski, Nonequilibrium evolution in scalar O(n) models with spontaneous symmetry breaking,Phys. Rev. D 65 (2002) 065019 [hep-ph/0109137] [INSPIRE]. [41] A. Arrizabalaga, J. Smit and A. Tranberg, Tachyonic preheating using 2PI-1/N dynamics and the classical approximation,JHEP 10 (2004) 017 [hep-ph/0409177] [INSPIRE]. [42] A. Arrizabalaga, J. Smit and A. Tranberg, Equilibration in φ4theory in 3+1 dimensions,Phys. Rev. D 72 (2005) 025014 [hep-ph/0503287] [INSPIRE]. [43] K. Kainulainen and O. Koskivaara, Non-equilibrium dynamics of a scalar field with quantum backreaction,JHEP 12 (2021) 190 [arXiv:2105.09598] [INSPIRE]. – 15 – JHEP10(2024)009 [44] K. Kainulainen, O. Koskivaara and S. Nurmi, Tachyonic production of dark relics: a non-perturbative quantum study,JHEP 04 (2023) 043 [arXiv:2209.10945] [INSPIRE]. [45] J.M. Cornwall, R. Jackiw and E. Tomboulis, Effective Action for Composite Operators,Phys. Rev. D 10 (1974) 2428 [INSPIRE]. [46] J. Berges, Introduction to nonequilibrium quantum field theory,AIP Conf. Proc. 739 (2004) 3 [hep-ph/0409233] [INSPIRE]. [47] D.G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, The art of simulating the early Universe — Part I,JCAP 04 (2021) 035 [arXiv:2006.15122] [INSPIRE]. [48] D.G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, CosmoLattice: A modern code for lattice simulations of scalar and gauge field dynamics in an expanding universe,Comput. Phys. Commun. 283 (2023) 108586 [arXiv:2102.01031] [INSPIRE]. [49] D.G. Figueroa, A. Florio, T. Opferkuch and B.A. Stefanek, Lattice simulations of non-minimally coupled scalar fields in the Jordan frame,SciPost Phys. 15 (2023) 077 [arXiv:2112.08388] [INSPIRE]. [50] G. Laverda and J. Rubio, Ricci reheating reloaded,JCAP 03 (2024) 033 [Erratum ibid. 06 (2024) E01] [arXiv:2307.03774] [INSPIRE]. [51] D. Bettoni, A. Lopez-Eiguren and J. Rubio, Hubble-induced phase transitions on the lattice with applications to Ricci reheating,JCAP 01 (2022) 002 [arXiv:2107.09671] [INSPIRE]. [52] G. Laverda and J. Rubio, The rise and fall of the Standard-Model Higgs: electroweak vacuum stability during kination,JHEP 05 (2024) 339 [arXiv:2402.06000] [INSPIRE]. – 16 –