Perturbative reheating and thermalization of pure Yang-Mills plasma

We investigate the thermalization of high-energy particles injected from the perturbative decay of inflaton during the pre-thermal phase of reheating in detail. In general, thermalization takes a relatively long time in a low-temperature plasma; therefore, the instantaneous thermalization approximation is not justified, even for the reheating of the Standard Model (SM) sector. We consider a pure Yang-Mills (YM) theory as an approximation of the SM sector or a possible dark sector, considering the Landau-Pomeranchuk-Migdal effect, a quantum interference effect in a finite temperature plasma. We perform the first numerical calculation to solve the time evolution of the system, including the redshift due to the expansion of the Universe, and show the details of the temperature evolution near the maximum and the behavior of the quasi-attractors at later times. The maximal temperature $T_\text{max}$ and time scale $t_\text{max}$ are determined quantitatively, such as $T_\text{max} \simeq 0.05 \times (\Gamma_I M_\text{Pl}^2/m_I^3)^{2/5} m_I$ and $t_\text{max} \simeq 2 \times 10^3 \times (\Gamma_I M_\text{Pl}^2/m_I^3)^{-3/5} m_I^{-1}$ in the SM-like system, where $m_I$ and $\Gamma_I$ are the mass and decay rate of inflaton. We also provide a similar formula for pure $\operatorname*{SU}(N)$ and $\operatorname*{SO}(N)$ YM theories for general values of $N$ and coupling constant $\alpha$, including $T_\text{max} \propto \alpha^{4/5}$ and $t_\text{max} \propto N^{-2} \alpha^{-16/5}$ behaviors and their numerical coefficients. The thermalization occurs in a finite time scale, resulting in a lower maximal temperature of the Universe after inflation than that under the instantaneous thermalization approximation.


Introduction
The Big Bang cosmology has been confirmed by the success of the Big Bang nucleosynthesis (BBN) theory [1,2] and the observation of cosmic microwave background (CMB) [3,4].After the Big Bang, the temperature of the background plasma decreases via the redshift due to the expansion of the Universe.The Universe might experience several phase transitions during cosmological thermal history, including the electroweak and QCD phase transitions.Several models of physics beyond the Standard Model (SM) also predict other cosmological phase transitions.In general, phase transitions lead to rich phenomenology in cosmology, such as the emission of gravitational waves (GWs) [5][6][7][8] and formation of topological defects [9].The dynamics of topological defects also lead to the emission of GWs [9][10][11][12][13][14][15][16][17][18][19][20][21].These GW signals are outstanding signatures of physics beyond the SM that is not accessible by collider experiments, and therefore, understanding the detailed thermal history is important to reveal the physics beyond the SM.For this purpose, the maximal temperature of the Universe is an important quantity.If the maximal temperature is higher than the critical temperature of a phase transition, one can expect that the Universe experiences the phase transition as it cools down due to the redshift.
The several problems of Big Bang cosmology, related to its initial condition, can be addressed by inflation [22][23][24].After inflation, the Universe is reheated by the decay of inflaton into radiation.The temperature of radiation reaches its maximum at a specific time and then cools down due to the redshift.One may assume instantaneous reheating, in which the radiation-dominated era is followed by inflation without an inflaton-dominated era.However, this scenario is considerably simplified in several cases.If the inflaton couples to radiation without any suppression, its decay is considerably fast and might cause a non-perturbative process called preheating [25,26].On the contrary, if the inflaton couples to radiation feebly, then the inflaton tends to dominate the Universe for a while and then slowly decays into radiation perturbatively.In this paper, we focus on the latter scenario.
The slow decay of the inflaton into radiation does not imply that the radiation is thermalized instantaneously after production.Namely, instantaneous thermalization is not generally justified and the thermal history during the pre-thermal phase is not trivial.Let us specifically consider a suppressed decay of inflation into gluons, e.g., via a (Planck-suppressed) higher-dimensional operator.For a typical inflation model, the inflaton mass is very large.♮1 The decay of such a heavy particle produces very high energy particles, such as gluons, in background plasma.The thermalization of high-energy gluons is quite non-trivial because they should lose significant amount of energy to be thermalized and absorbed into thermal plasma [30][31][32][33].In fact, our previous studies have shown that the thermalization-time scale is generally significantly longer than the Hubble-time scale in the early stage of reheating, and the maximal temperature of the Universe may be significantly suppressed [34][35][36][37][38][39].Detailed thermalization also provides a novel dark matter (DM) production mechanism in the pre-thermal phase [39][40][41][42].Several previous studies have focused on the DM production, in which case the particle distribution is reduced to a stationary solution, and the analysis can be significantly simplified.♮1 Although the decay rate is sufficiently small to justify the perturbative analysis, the temperature of plasma can be greater than the inflaton mass.In this case, the decay rate gets modified by the thermal effects [27][28][29].
The bottleneck process for thermalization is the splitting of a high-energy gluon into low-energy daughter particles.The splitting process experiences an interference between wavepackets, known as the Landau-Pomeranchuk-Migdal (LPM) effect [43][44][45][46][47][48], and therefore has a strongly suppressed rate of Γ split =  2 √︁  3 /, where  is the energy of the parent particle, and  is the (effective) temperature of the background plasma.♮2 Comparing this with the Hubble expansion rate, the maximal temperature of the Universe can be estimated as  max ∼  4/5 (Γ   2  Pl / 3  ) 2/5   , where  Pl is the reduced Planck mass, Γ  is the inflaton decay rate, and   is the inflaton mass [34,35].This can be several orders of magnitude smaller than the naive result under the instantaneous thermalization approximation.This qualitative discussion demonstrates the importance of a detailed investigation of thermalization for inflaton-decay products to pin down the maximal temperature of the Universe.
In this paper, we numerically solve the detailed Boltzmann equation to calculate the time dependence of temperature during the pre-thermal phase, considering the LPM effect.To minimize the numerical cost, a pure Yang-Mills (YM) theory is considered.This is a good approximation for the SM sector and also motivated by a thermalization of the dark sector that may explain DM.We provide a quantitative formula for the maximal temperature and thermalization time scale.This is the first work on the quantitative time evolution of thermal plasma during the pre-thermal phase.
The organization of this paper is as follows.In Sec. 2, we briefly review the thermalization process of a high-energy particle under the cosmological expansion.Moreover, a qualitative estimation of thermal history is explained.In Sec. 3, we explain our numerical method and show our results.The results are consistent with the qualitative discussion, but provide more information, including numerical prefactors.We consider a model that mimics the SM sector and a model of the pure YM dark sector.In the latter case, the gauge group  is assumed to be SU() and SO().A formula for maximal temperature is also shown in a large  and small gauge coupling limit.Section 4 is devoted to discussion and conclusions.

Kinetic equations 2.1 Warmup
Let us start our discussion by neglecting the thermalization of the produced particles as a warmup.In this case, the relevant Boltzmann equation is given as follows: where  is the physical momentum, and  is the Hubble parameter.The distribution function of gluon per one degree of freedom is denoted by  , that is, the total number density of gluon is obtained by multiplying the degrees of freedom.The source term is (approximately) given by a delta function at  =  0 .Throughout this paper, we consider the case where the primary particles originate from the two-body decay of a heavy particle, e.g., inflaton, with number density   (), mass   , and decay rate Γ  .
The source term is expressed as follows: where   is the degrees of freedom of the gluon.If the heavy particle behaves as pressureless matter, we have where () is the scale factor, and  0 is the reference time.
Eq. (2.1) can be solved, such as ) where  = / and ()/ 0 = (/ 0 )  .In this paper, we focus on the thermalization in the matter-dominated epoch, in which  = 2/3.Hereafter, we neglect the exponential factor in (2.5) because we are primarily interested in the thermalization that occurs in the regime of Γ  <  ().The momentum ( inf )  0 /() in the second Heaviside theta function in (2.5) represents the redshifted momentum of gluons generated at the end of inflation  inf .We take ( inf ) → 0 for simplicity throughout this paper because its precise value does not affect our result qualitatively.

Qualitative discussion
We move on to the discussion on the thermalization of pure YM plasma produced by the decay of a heavy particle, e.g., inflaton.Before presenting the concrete kinetic equations, we briefly summarize the basic assumptions and qualitative discussion of reheating and thermalization.See also Ref. [35].As emphasized in the introduction, we focus on the case in which the decay rate of inflaton, Γ  , is sufficiently small such that the reheating can be well approximated by the perturbative inflaton decay.It is convenient to use the following combination for the decay rate Γ  : which is assumed to be smaller than unity.This is typically expected at least when the inflaton decay rate is smaller than that induced by the dimension-five Planck-suppressed operators, such as  5      / Pl and  5    G / Pl with  5 ≪ 1, where  represents the inflaton.In this case, the reheating temperature is significantly smaller than the inflaton mass.Hence, we generically expect that the temperature of the "dilute" plasma, which already exists before the completion of reheating, is smaller than the typical energy of gluons just after their production, at least in the later stage of pre-thermal phase.
The bottleneck process of thermalization in such cases is given by the splittings of high-energy gluons into the background plasma.In general, a splitter of momentum  can emit a splittee at a scale of  <  only after  ⩾ / 2 ⊥ with  ⊥ being a momentum of  transverse to the direction of  because otherwise, quantum mechanical interference between the splitter and splittee prevents splittee formation.Suppose there exist soft-thermalized populations of gluons with a temperature  s , whose condition is specified later.The interactions with the soft thermal plasma induce random diffusions of the transverse momentum, which is estimated as  2 ⊥ ∼ q ∼    2  3 s .♮3 Here, we have included the factor  s that depends on the degrees of freedom responsible for the transverse diffusion.One may estimate its dependence as where the summation is taken over species  contributing to the transverse diffusion, the degrees of freedom for  are denoted by   , the dimension and normalization of representation for  are   and   , respectively, and the quadratic Casimir for the adjoint representation is  A .For pure YM plasma,  s ∼  2 A .Combining these two estimations, we find the formation momentum below which the splittees can be emitted for a given (2.9) Equivalently, the formation time before which the splittees of momentum  cannot be emitted is given as follows: Once this condition is met, the splittee of momentum  can be formed with a probability of , which leads to the LPM-suppressed splitting rate, Γ LPM () ∼  A / form (). On the other hand, for  >  form (), the interactions with the medium are not sufficient to build up the transverse momentum, whereas a large-angle emission of  2 vac ⩾ 1/() is always allowed quantum-mechanically with a probability of  (up to a logarithmic factor), that is, Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) vacuum shower.Hence, the splitting rate, where a splitter of momentum  splits into daughters of  and (  − ), can be obtained as (see also Ref. [49]) where logarithmic factors are dropped for simplicity.
The emitted splittees of momentum  immediately cascades and participate in the soft thermal plasma within the Hubble time if Γ LPM () >  ∼ 1/, which reads (2.12) The energy conservation implies  s () ∼  split () Γ     () with  s being the energy density of the soft sector.Assuming the thermalization of the soft sector, i.e.,  s ∼  ★  4 s with  ★ being the effective relativistic ♮3 Here, we neglect the running of  for simplicity.See the later discussion on how the running modifies the estimation.
degrees of freedom in the soft thermal plasma, e.g.,  ★ =   = 2 A for pure YM plasma, we obtain Now, the condition for the soft-sector thermalization can be derived.To maintain the thermal distribution, the large-angle scatterings among the soft populations should be much faster than the cosmic expansion Subsequently, the gluons generated by the decay of inflaton immediately break up into the soft thermal plasma.The temperature of the soft sector is then given by  s ∼   Γ    , which yields for  max <  < Γ −1  .The reheating is completed at  ∼ Γ −1  , and the temperature at that time is often referred to as the reheating temperature,  R ∼  −1/4 ★ √  Pl Γ  .To sum up, in addition to the hard distribution given in Eq. (2.5), the splittings of high-energy gluons yield the cascades towards the soft thermal plasma.Its distribution can be estimated as ) Note that the hard population given in Eq. (2.5) vanishes for  >  max .The temperature of the soft sector is given by ♮4 Here we write the dependence of degrees of freedom for the large-angle scattering rate as   for notational simplicity.
Although this is correct at least for pure YM plasma, it is different from   for a general case.
which is maximized at  ∼  max as The rest of this paper is devoted to refining this qualitative understanding through numerical simulations.

In-medium splitting function
The kinetic equation that describes the reheating and thermalization via the production of high-energy gluons after inflation is summarized as follows [39,50]: where the splittings of high-energy gluons are described by C 1↔2 , and the elastic scatterings are represented by C 2↔2 , which lead to the thermalization of the soft sector.Throughout this paper, we focus on the time scale much longer than the elastic scatterings, i.e.,  ≫  soft .Hence, we simply assume that the splittees of momentum  immediately get thermalized as long as  IR >  and impose the thermalization by hand without explicitly including C 2↔2 in the numerical simulations.The results of our numerical simulations should be insensitive to the choice of  IR as long as  s () ≲  IR ≪  split ().
Throughout this paper, we consider the case where the cascades of high-energy gluons are dominated by the splittings into low-energy gluons.This restriction is trivially fulfilled for weakly coupled pure YM theories and holds for the high-energy gluons in the SM.The splitting of gluons is encoded in C 1↔2 as (2.22) The splitting function of gluons in the presence of thermal plasma is ) ↔ () ≡ ) where the fine structure constant  is evaluated at a scale , the well-known DGLAP splitting function for gluons is denoted as  (vac)  ↔ (), the degree of freedom is   = 2 A , the dimension of the adjoint representation is  A , and its quadratic Casimir is  A ; for SU() and SO() gauge theories, ( A ,  A ) = ( 2 − 1, ) and ( ( − 1)/2, 2 − 4), respectively.The transverse momentum squared at the formation time is denoted by  2 ⊥ , which is developed by interactions with the soft thermal plasma of gluons.It can be computed by solving a self-consistent equation, as done in Refs.[51], whose result at the leading log is given by with   being the one-loop -function coefficient, e.g.,   = −11 A /3 in the pure YM theory.The factor N is proportional to the number density responsible for the transverse momentum diffusions and is given as follows: where the summation is taken over species  contributing to the transverse diffusion, the distribution function of the soft population is denoted as  s , and the Riemann zeta function at 3 is  (3) ≃ 1.20206.Note that the transverse momentum diffusion is dominated by the soft thermal plasma for  ≫  soft .In the second line, we used the fact that the normalization of the adjoint representation fulfills  A =  A .We also define the following quantities: ) In the case of the SM plasma, all the SM particle contributions can be included in the soft sector to the high-energy gluon cascades such that   = −7, N = 15 (3) 3 / 2 , and  2  = 8 2 .Finally, we briefly discuss how the equation given in this section is related to the qualitative discussion in the previous section.For  ≪ , Eq. (2.24) can be approximated with  2 ⊥ ( ; ,  − ) ∼ √︁  q.It is consistent with the diffused transverse momentum at the formation time given in Eq. (2.10), i.e.,  2 ⊥ |  form ∼ q form ∼ √︁  q.The splitting function is expressed as  ↔ ( ; ,  − )/  ∼  A  √︁  q for  ≪ .The corresponding splitting rate is then given as Γ LPM () ∼  ↔ ( ; ,  − )/(  ) ∼  A  √︁ q/, which is consistent with Eq. (2.11).

Numerical simulations 3.1 Numerical method
We numerically solve Eq. (2.21) without the C 2↔2 term, under the assumption that the splittees of momentum  with  <  IR are thermalized immediately via 2 → 2 elastic scattering process.The comoving temperature of the plasma increases as the energy is injected into the modes  <  IR .The physical temperature is calculated as follows: where  ★ is the relativistic degrees of freedom in the thermal plasma.We use  ★ = 106.75 for the SM sector and  ★ =   for the pure YM theory.The total energy injected into the radiation,  r, tot (), is given by We define the energy density for non-thermal contributions as We take  IR = 3 s () in our numerical simulations.The results do not change qualitatively when  IR is varied by a factor of O(1).We rescale the dimensionful parameters and use  s ( 0 ) ≡  0 = 1 in numerical simulations without any loss of generality.
The source term S in Eq. (2.21) can be represented by the initial distribution (2.5).Because the numerical time is limited, our numerical calculation starts from a finite time scale  0 .This implies that the distribution at a small momentum scale is deformed as the second line of Eq. (2.17) (∝  −7/2 ), ♮5 whereas the hard mode at a high momentum remains in the form of Eq. (2.5) (∝  −3/2 ).To include this fact, we adopt the following algorithm in our numerical simulations.We first take an ansatz  h ( ,  0 ) =  h (  0 )( / 0 ) −3/2 for the initial distribution and evolve it by a first time step,  =  0 +.Then we determine  th,0 as a minimal value of momentum satisfying ln [  h ( ,  0 + )/  h ( ,  0 )] < 0.1.Then we perform the main numerical simulations by replacing the initial distribution  h ( ) such that for  ≥  th,0 .

(3.5)
This initial condition allows the distribution function not to change significantly for the first time step, which is important to ensure the stability of numerical simulations.
The splitting term C 1↔2 is given by (2.22).We take an IR cutoff for the splitting term as the initial temperature of the system  0 (= 1), assuming that temperature does not change by many orders of magnitude during the simulation time scale.It is convenient to use a comoving momentum p ≡ () / 0 and temperature Ts () ≡ () s ()/ 0 in numerical simulations.We denote pmax () ≡ () ( 0 )  0 , (3.6) ♮5 We neglect the deformation of the distribution of the soft particles for  >  form (∝  −3 ) for simplicity.

Results of the SM
To describe the thermalization in the SM sector, we consider the thermalization of gluon (i.e., SU(3) gauge field) into the SM thermal plasma.Namely, we take  ★ = 106.75,G = SU(3),   = 16,   = 8,   = 3, , and (  ) = 0.118 with   being the -boson mass.We take  0 = 10 3 GeV = 1 in the numerical simulations, where the dependence on this choice is only because of the renormalization group (RG) running of the gauge coupling constant.
In this case, there remains three parameters that should be specified to perform the numerical simulations: (3.9) Here,  h (  0 ) should be consistently determined by the first line of (2.19) with  =  0 .This implies that We are interested in the regime around  ∼  max ; thus, we take  0 to be of order but smaller than  max .Moreover, we need to ensure  soft ≪  0 .Specifically, we take ,  h (  0 ) 2 0  3 0 ) = (5 × 10 3 , 2, 10 5 ), (10 4 , 1, 10 5 ) respectively.The brown solid curve corresponds to (5 × 10 3 , 5, 10 5 ).from (2.12), where we use (2.19) to eliminate Γ or  h (  0 ).The lower bound on   originates from the condition of  soft ≪  0 .In the numerical simulations, we primarily take and  grid = 2 × 10 4 .We also perform numerical calculations by changing  0 and  0 by a factor of a few to determine whether the result depends only on physical parameters.
Figure 1 shows the temperature as a function of time, where the variables are rescaled by These dependencies on the initial parameters are expected from the analytic estimations of Eqs.(2.20), (2.15), and Eq.(2.7).The blue solid curve represents the case with the initial condition mentioned above.The brown solid curve represents the case with  0 = 5 × 10 3 ,  0 = 5 1/2 0 ,  h (  0 ) = 10 5  −2 0  −3 0 , and  grid = 2 × 10 4 .The dashed curve represents the case with  0 = 10 4 ,  0 =  1/2 0 ,  h (  0 ) = 10 5  −2 0  −3 0 , and  grid = 4 × 10 4 .All results are in good agreement with each other, except for the regime around  =  0 .The discrepancy results from the fact that the initial distribution (and  h (  0 )) is chosen by hand for soft modes.However, because the energy injection into thermal plasma is dominated by the contribution from hard modes, the late-time distribution is not affected, i.e., implying the basin of attraction.The results for other initial conditions are provided in the Appendix.The red and green dotted lines represent the fitting functions for  (): All results are in good agreement with the analytic estimations.The asymptotic behavior for  ≫  max is further confirmed by simulations with different initial conditions as we show in Fig. 4 in Appendix.From Fig. 1, the numerical coefficients of the maximal temperature and corresponding time can be determined as follows: These results confirm the analytic estimations (2.20) and (2.15), and the numerical prefactors are determined from our numerical simulations.Note that  and  ★ are absorbed into the numerical factors because we substitute the SM values for them.The blue curves in Fig. 2 show the spectra at the end of numerical simulations.In the left panel, the initial spectrum (3.5) is represented by the blue dashed curve.The solid (dashed) brown curve represents the thermal spectrum at the end (beginning) of numerical simulations.See Eqs.(2.17) and (2.18) for analytic estimations.In the right panel, the distributions are rescaled by Γ( / 0 ) −7/2 , where Γ ≡ 4 2 Γ    () is defined in our previous paper [39].The three curves correspond to the ones used in Fig. 1 and almost overlap with each other.At a late time, Γ split ≫ , in which case the spectrum can be represented by a stationary solution of for  ≪  0 .The red dotted curve, which is also overlapped with our numerical results, is the result for the stationary solution of the Boltzmann equation in the limit of Γ split ≫ .All results are consistent with the stationary solution, which justifies  ≫  max at the end of numerical simulations.The deviation from the dependence of ∝  −7/2 at a small / 0 may be due to the IR cutoff used in our numerical simulations.
For the stationary solution, the IR cutoff need not to be introduced, as explained in Ref [39,40]; therefore, the stationary solution is almost exactly ∝  −7/2 for a small / 0 .
In the numerical simulations, we take ) and  grid = 2 × 10 4 , for (, ) = (3, 0.1), where  SM = 16 and  * ,SM = 106.75.Here, we define For the case with  ∼ 0.1, including the SM QCD, the difference of the gauge coupling constants at  0 and at  0 due to the RG running is not negligible.For example, for (, ) = (3, 10 −1 ) with G = SU().In this case,  in  soft is relatively larger than that used in  max .This results in a relatively weaker condition on the initial time, (3.11), by a factor of a few.In contrast, for a significantly smaller gauge coupling constants, the RG running is negligible.In this case, the lower bound on the initial time, given by the first inequality of (3.11), is more severe than the case of the SM.Therefore, we instead take  0 = 5 ×  (bm) 0 for (, ) = (5, 10 −2 ), (5, 10 −3 ), (5, 10 −4 ), and (10, 10 −3 ). Figure 3 shows the temperature as a function of time in pure YM theories.To show our results, we rescale parameters such as (3.28) We plot the results of SU() (left panel) and SO() (right panel) gauge theory with (, ) = (3, 10 −1 ) as blue solid curve and those with (5, 10 −2 ), (5, 10 −3 ), (5, 10 −4 ), and (10, 10 −3 ) as brown solid curves.The brown solid curves overlap and cannot be distinguished from each other.The dotted lines are the fitting functions of (3.17).The figure shows that all results are consistent with (3.17) after the rescalings.Moreover, the results of (3.18) and (3.19) can be applied to pure YM theories after the correction from the difference of the splitting rate, Γ split (  0 ), is included: where we substitute the SM parameters, such as Γ split (  0 ) ≃ 5.5 × 10 −4 ×  3/2 s  −1/2 0 , and  0 =   /2.This is a general result that should be applicable to any pure YM(-like) theories.
We also provide useful formulas for  max and  max in large  and small  limits.If  ≫ 1 and  ≪ 0.1, we can approximate where the prefactors depend on  0 logarithmically.Substituting these into (3.29) and (3.30), we obtain ) for G = SU() and ) for G = SO(), where we use  * =   .These results confirm the analytic estimations of Eqs.(2.15) and (2.20), where   ∼  2  ∼  2 and  ★ =   ∼  2 .

Discussion and conclusions
We have obtained a clear dynamical picture of the thermalization of pure YM plasma during perturbative reheating after inflation by numerically solving the Boltzmann kinetic equation by appropriately considering the LPM effect.Our results confirm the results of previous analytic studies based on quasi-equilibrium ansatz applicable to two regimes of  soft ≪  ≪  max or  max ≪ , but also describe the transient epoch between these regimes, providing a better estimation for the maximal temperature of our Universe.The maximal temperature obtained can be significantly smaller than the estimation based on the instantaneous thermalization approximation, highlighting the importance of understanding thermalization.Moreover, our results demonstrate the robustness of the quasi-equilibrium solution under the consistent change of the initial conditions.Furthermore, the late-time behavior is stable even when the initial conditions are applied far away from the consistent initial conditions.The implications of our results for particle cosmology are broad.In particular, the maximal temperature of the Universe is the key ingredient for understanding possible cosmological phase transitions.Recently, hidden pure YM theories have been gaining attention because they involve the glueballs as a candidate of DM [52][53][54][55][56][57][58][59][60][61][62][63][64][65] and would lead to the first-order phase transitions [66][67][68][69], possibly accompanied by cosmic strings [70][71][72].Our results are essential to understanding their implications, such as the prediction of the GW spectrum as well as the condition to obtain the confinement phase transition after inflation.
As discussed in the main text, we believe our results can be applicable to the SM plasma when the inflaton dominantly decays into the SM gluons of SU(3).This expectation is based on our previous study [39], which showed that the thermalization is dominated by the gluons, although the analysis is restricted to the quasi-equilibrium regime.The complete dynamical analysis of the SM plasma, including the cases in which the inflaton decays into other SM species, is worthwhile to investigate in future studies.the maximal temperature.Still, all results reach the maximal temperature within the time scale of order  max , even if the simulation starts at a later time.All results agree at a later time ( ≫  max ), at which high-energy particles thermalize within the Hubble-time scale.We expect that the cases with  h (  0 )/(10 5 ×  −2 0  −3 0 ) = 1 (blue curves) are consistent initial conditions in which the ambient plasma is generated by the energy injection via the thermalization of high-energy particles.These cases agree with each other except for  ∼  0 , even if we change the value of  0 / 1/2 0 .This supports the fact that it consistently starts within the attractor regime of (2.19).This is not the case for different values of  h (  0 )/(10 5 ×  −2 0  −3 0 ), as shown in the figure.However, we note that the other initial conditions may also be interesting in some of the cases mentioned above.Our numerical simulations can also be used to analyze such cases.

Figure 2 :
Figure 2: Spectra at the end of numerical simulations.(left) The blue/brown solid curves represent the hard/thermal spectra respectively.We also show the initial spectra with the dashed lines in the same color.(right) The red dotted curve represents the stationary solution corresponding to Γ split / → 0 given in Eq. (3.21).The color codings for the blue solid/dashed and brown curves are the same as Fig. 1.

Figure 4 :
Figure 4: Same as Fig. 1 but with different initial conditions.See Eq. (A.1) for their color codings.
.14)Throughout this paper, we restrict ourselves to  ≫  soft .With time,  split grows continuously, and eventually  split () >   , or equivalently Γ split (  ) > , which occurs at Physical temperature as a function of physical time.All curves are rescaled according to(3.15)and(3.16).The red and green dotted lines represent the analytic dependence for  s () in the regimes of  ≪  max ( s () ∝ ) and  ≫  max ( s () ∝  −1/4 ), respectively.The blue solid/dashed curves represent the numerical results for (  0 ,  0