Perturbative reheating in Sneutrino-Higgs cosmology

The theory of perturbative reheating in the Sneutrino-Higgs cosmology of the B-L MSSM is presented. It is shown that following an epoch of inflation consistent with all Planck2015 data, the inflaton begins to oscillate around its minimum at zero and to reheat to various species of standard model and supersymmetric matter. The perturbative decay rates to this matter are computed, both analytically and numerically. Using these results, the Hubble parameter and the relative energy densities for each matter species, including that of the inflaton, are calculated numerically. The inflaton energy density is demonstrated to vanish at an energy scale of O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(1013) GeV, signaling the end of the period of reheating. The newly created matter background is shown to be in thermal equilibrium, with a reheating temperature of ≃ 1.13 × 1013 GeV. To allow for a B-L breaking scale sufficiently smaller than the reheating scale, we extend the statistical method of determining the soft supersymmetry breaking parameters developed in previous work. The result is that one can determine a large number of phenomenologically realistic initial conditions for which the B-L breaking scale is an order of magnitude or more smaller than the reheating scale.


JHEP09(2018)001
The introduction of Higgs inflation in the standard model of particle physics [1][2][3][4][5][6], emphasized the important idea that the "inflaton" associated with inflationary cosmology might well be -and, perhaps, should be -a fundamental scalar field in the standard model. In the non-supersymmetric case, this could only be the neutral Higgs scalar. A careful analysis of this possibility led to the result that an acceptable theory of inflation could potentially arise in this context -but only at the cost of assuming several unnatural features. For example, the square of the Higgs magnitude must necessarily be coupled to the curvature scalar in the Lagrangian density with an unnaturally large coupling parameter. Generalizing this idea to the supersymmetric standard model (MSSM), and some variants thereof, potentially extends this idea by introducing a large number of scalar fields into the theory, a subset of which might play the role of the inflaton. Although at first seeming promising, these models [7][8][9][10] continued to exhibit various fundamental problems. Not the least of which is the fact that the plethora of scalar fields now greatly complicates that potential energy. In particular, it now becomes very difficult to obtain a potential for the "inflaton" that is stable with respect to the other directions in field space. That is, the inflating field usually rolls off into another scalar direction, rapidly losing potential energy and ending the inflationary period far short of the required 60 e-foldings. Fortunately, there is one specific supersymmetric theory -called the B-L MSSMthat naturally solves this difficulty. The B-L MSSM is a vacuum solution of heterotic M-theory [11][12][13][14][15][16][17][18][19] whose observable sector has precisely the spectrum of the MSSM, where each of the three lepton families contains a right-handed neutrino chiral supermultiplet. The theory differs from the MSSM in that the gauge group is enhanced from the SU(3) C × SU(2) L × U(1) Y group of the standard model to SU(3) C × SU(2) L × U(1) Y × U(1) B−L . That is, the global B-L symmetry of the MSSM is gauged. This theory is known to remain anomaly free and is the simplest possible extension of the MSSM. This theory was introduced in [20], and its properties studied in a series of papers [21][22][23][24][25][26][27][28][29][30][31]. Importantly, it was shown in [30] that for a very wide and diverse range of initial conditions, the B-L MSSM is completely consistent with all phenomenological data; that is, the B-L symmetry is appropriately radiatively broken, electroweak symmetry is also radiatively broken with the correct values of the Z 0 and W ± masses, all sparticles exceed the present experimental bounds and, remarkably, the prediction for the mass of the neutral Higgs boson is within 2σ of the present ATLAS bounds from CERN [32][33][34]. Thus, if a successful inflationary theory can be developed in this context, it would simultaneously be consistent with all low energy phenomenology.
Remarkably, a completely successful theory of inflationary cosmology can be developed within this context. The reason for this is relatively straightforward. The potential energy of the scalar field has three contributions; the D-term potential V D = 1 2 a D 2 a , the Fterm potential V F and contributions from the soft supersymmetry breaking terms V soft . As will be discussed below, the F-term potential essentially vanishes after the beginning of inflation. What is remarkable about the B-L MSSM is that, for a precise combination -denoted by φ 1 -of the right-handed sneutrino, the left-handed sneutrino and the JHEP09(2018)001 neutral up Higgs field, all D a contributions also vanish, rendering V D = 0. This fact, which is due to the existence the right-handed sneutrino, lowers the entire D-flat potential into a "valley" that is stable against roll-off into any other scalar direction. It follows that the inflaton stably inflates along the V soft potential -completely consistent with all Plank2015 cosmological data. This inflationary theory, which we have named "Sneutrino-Higgs inflation" was developed in detail in [35]. In fact, this theory is very similar to the so-called "no scale" supergravity theory developed in a different context in [36][37][38][39]. What we believe is important about Sneutrino-Higgs inflation in the B-L MSSM is that it is, as stated above, entirely consistent with all low energy particle physics phenomenology. That is, not only is it an accurate inflationary theory, it is within the context of realistic low-energy particle physics as well.
With this in mind, in the present paper we move beyond the purely inflationary epoch and study the precise theory of perturbative reheating of the Sneutrino-Higgs theory. The subject of reheating in inflationary models has been studied extensively; see for example [40][41][42][43][44][45][46][47][48][49][50][51][52][53]. As discussed in those papers, reheating can occur after the inflationary epoch via preheating (that is, parametric resonance), via perturbative decay of the inflaton or as a combination of both of these processes. Which of these two processes dominates reheating is highly model dependent. Of all the theories discussed in the preceding references, only the models in [51][52][53] are closely related to the Sneutrino-Higgs inflation in the B-L MSSM discussed in our work. Specifically, the theories in [51][52][53] involve Higgs and sneutrino inflatons coupled to so-called "no-scale" N = 1 supergravity -very similar to Sneutrino-Higgs inflation in the B-L MSSM. Within this context, the first two papers found that perturbative reheating is the dominant process for both decay to standard model particles [51] as well as to the gravitino [52]. This was due to the fact that the inflaton in these papers was purely a right-handed sneutrino and, hence, the coupling to matter is via small Yukawa and gravitational interactions. In [53], however, the authors analyzed the case of a pure Higgs inflaton, which couples to matter via gauge interactions. They conclude that non-perturbative production may or may not be a significant component of reheating, depending on the specific parameters of the model. However, it is possible that preheating could be the more efficient reheating mechanism. Be that as it may, in the present work we will compute only the perturbative decays of the Sneutrino-Higgs inflaton, since this is expected to play a role in the reheating process. The possible effects of preheating, which may dominate over perturbative decays, will be explored elsewhere. As we will show, in the Sneutrino-Higgs theory, this reheating epoch is completely amenable to exact computation. The end result is that this theory also reheats in a technically determined manner and appears to be completely acceptable physically.
Specifically, we will do the following. In section 2, we give a brief review of the theory of Sneutrino-Higgs inflation presented in [35]. We outline the basic formalism and present the specific parameters that lead to a successful theory of inflation -completely consistent with the Planck2015 data [54]. Using the statistical method for choosing the soft supersymmetry breaking parameters introduced in [30], we show that there is a large and diverse set of phenomenologically valid initial conditions that are consistent with the appropriate cosmological parameters. We conclude this section by presenting a new result JHEP09(2018)001 -specifically, we compute the range of values of the B-L breaking scale in this context and plot this as a histogram against the number of valid points at a given B-L scale. We find that the smallest B-L scale that one can attain in this formalism is 2 × 10 12 GeVand this for only a very small number of valid initial points. Although this is marginally acceptable, we would like to modify our formalism in such a way as to attain smaller values for the B-L scale. The reason is that we expect the reheating temperature to be of O(10 13 GeV) and would, for simplicity, prefer that B-L breaking, and, hence, the breaking of baryon number, lepton number and R-parity, to occur at a much smaller scale. This allows one to separate the reheating epoch from the period of baryo-and lepto-genesis [55]. This is carried out in detail in appendix A, where we introduce an extension of our previous formalism that allows one to lower the B-L scale arbitrarily, while remaining completely phenomenologically realistic. In this appendix, we give an explicit example where the B-L scale is reduced to a range from a low of 10 6 GeV to well above 10 13 GeV -all this occurring for phenomenologically acceptable initial conditions. We also give an explicit formula for the degree of fine-tuning required to set the B-L scale to any given value. Such a wide range is not required for cosmology. Hence, section 3 is devoted to explaining this new formalism, and then using it to alter the B-L scale to a more reduced range of values -namely from 10 10 GeV to 10 12 GeV. We find, in this context, that the greatest number of phenomenologically valid initial points occurs for a B-L scale of 10 11 GeV and, hence, when specificity is required later in our analysis, we will choose this as the value of B-L symmetry breaking.
The remainder of the paper concentrates on the post-inflationary epoch and reheating. In section 4, we give a detailed discussion of the classical behavior of both the inflaton and the Hubble parameter -ignoring for the time being the quantum mechanical decays of the inflaton into matter. We begin by calculating these quantities numerically and show that at a given time after the end of inflation, denoted by t osc , the inflaton begins to oscillate around its minimum at zero. Then, using an iterative method, we analytically compute both the inflaton and Hubble parameter to a high degree of approximation. We show that this analytic solution becomes well-defined after a time, t MD > t osc , which marks the beginning of the period of matter-domination. The numerical and analytic results are plotted simultaneously and shown to very closely approximate each other for t > t MD . In section 5, using these results, we compute the decay of the inflaton into matter and, hence, reheating. We show that there are four major decay channels -up-type standard model particles, charginos, neutralinos and gauge bosons -and analytically calculate the decay rates Γ of the inflaton into each of these species. We present the details of this analysis both in section 5 and in appendix C. In section 6, we analyze the full set of differential equations for the inflaton and Hubble parameter -now sourced by the decay rates and associated energy density of the matter species into which the inflaton decays -as well as the differential equations for the evolution of the matter itself. These equations are solved numerically for t > t MD . The time-dependent behavior of the individual decay rates, as well as the behavior of the Hubble parameter, are plotted. Using these results, the fraction of the individual energy densities, ρ i ρ total , including the inflaton density, are computed and plotted. The time of reheating, t R , is defined to be the time at which the JHEP09(2018)001 inflaton density vanishes -indicating that its energy has then been entirely transferred to matter. Finally, in section 7 we show that by t R all the newly created matter species are in thermal equilibrium, and calculate the associated reheating temperature.

Supergravitational B-L MSSM
In the flat superspace limit, the particle content of the B-L MSSM is precisely that of the N = 1 supersymmetric MSSM; that is, three families of quark and lepton chiral supermultiplets, each family with a right-handed neutrino chiral multiplet, and a pair of Higgs, Higgs conjugate chiral superfields. The gauge group, however, is extended from the standard model gauge group to (2.1) The additional Abelian gauge group U(1) B−L is anomaly free and arises from the fact that the B-L MSSM is a low energy vacuum of heterotic M-theory [11,15,17,22] compactified to four-dimensions on a specific Calabi-Yau threefold [58] with a given holomorphic vector bundle [20]. As discussed in [30], it is convenient (and physically equivalent) to use the group U(1) 3R , with generator T 3R = Y − 1 2 (B − L) in place of U(1) Y , since this eliminates kinetic mixing of the associated field strength with that of B-L. 1 The formalism for coupling of a generic four-dimensional compactification of M-theory to N = 1 supergravity was presented in [59], and can easily be applied to the B-L MSSM. We note that the Planck mass used in this formalism is the so-called "reduced" Planck mass given by M P = 2.435 × 10 18 GeV. For the bosonic sector of the theory, this was carried out in [15]. Suffice it here to say that, to order κ 2/3 in heterotic M-theory, by appropriately redefining the matter scalars C i and setting the Planck mass M P = 1, the matter scalar Kähler potential and each of the gauge kinetic functions in the Lagrangian become respectively. Recalling that the matter scalar kinetic energy terms in the Lagrangian are given by it follows that for small values of C i these kinetic energy terms are canonically normalized. Importantly, however, this is no longer true when one or more of these fields are of the

JHEP09(2018)001
order of the Planck scale. Note that throughout this paper, we will often, for notational simplicity, set M P = 1. However, there are sections where the physical content is more transparent when expressed in mass dimensions of GeV. Which notation is being used will be explicitly indicated.
In addition to the kinetic terms, the B-L MSSM Lagrangian contains a holomorphic superpotential given by where both generational and gauge indices have been suppressed. As discussed in [30], the Yukawa couplings can be taken to be flavor diagonal and real for the purposes this analysis. The µ parameter can also be chosen to be real without loss of generality. This superpotential leads to the F-term contributions to the potential energy for matter scalars in the Lagrangian, given by the usual expression Similarly, eliminating the auxiliary fields in the vector superfields associated with the SU(3) C , SU(2) L , U(1) 3R and U(1) B−L gauge groups leads to the D-term contributions to the matter scalar potential energy given by and T r (a) , r = 1, . . . , dim G a are the generators of the group G a . Both V F and V D were computed for all matter scalars in [35]. We refer the reader there for details. Finally, it is assumed in [30] that supersymmetry is spontaneously broken in a hidden sector. This leads to "soft" supersymmetry breaking terms in the observable low-energy Lagrangian. These contain a final contribution to the matter scalar potential energy given by We will make the usual assumption that each of the dimensionless couplings for the cubic terms is proportional to the associated Yukawa coupling. Since it will become important in the analysis of reheating given later in this paper, we point out that, in addition to the scalar terms in (2.8), the soft supersymmetry breaking terms also include mass terms for each of the gauginos given by

JHEP09(2018)001
The fermionsg,W ,W R andB in (2.8) are the SU(3) C , SU(2) L , U(1) 3R and U(1) B−L gauginos respectively. Finally, we point out that the coupling of the B-L MSSM to N = 1 supergravity using the M-theory formalism given in [15], leads to a canonical Lagrangian for the metric tensor g µν . That is, we find that the pure gravitational action is simply given by It follows that in the Sneutrino-Higgs inflation theory, we do not require any "non-canonical" coupling of matter to the curvature tensor R.

The inflaton and the inflationary potential
The inflationary potential in [35] was constructed as follows. In order to help guarantee the stability of the potential energy, we first searched for a direction in scalar field space for which all of the above D-terms vanish. We were naturally lead to the field space with all other fields set to zero. We note that only in a model such as the B-L MSSM with right-handed neutrino superfields would such a D-flat direction arise. To proceed, we defined three new fields φ i , i = 1, 2, 3 by The field φ 1 corresponds to the D-flat field direction while φ 2 and φ 3 are two orthogonal directions. We note that (2.13) and the associated quadratic soft mass squared is given by (2.14) Setting all fields to zero with the exception of φ 1 , the V D potential vanishes and the Lagrangian becomes

JHEP09(2018)001
Here Y ν3 is the third-family sneutrino Yukawa coupling and µ is the usual supersymmetric Higgs parameter. In [35] it was verified in detail that, modulo several acceptable caveats, the field φ 1 lies in a valley of the total potential energy. That is, the potential energy gets larger as one moves away from the pure φ 1 direction by making any other scalar field non-zero. We henceforth, for simplicity, considered only the real part of φ 1 , which we continue to denote by the same symbol. Since the value of φ 1 can be of Planck scale during inflation, it follows from (2.2) and (2.3) that its kinetic energy is not canonically normalized in this regime. To canonically normalize the kinetic term, we make a field redefinition to a real scalar ψ given by In terms of the new field ψ, Lagrangian (2.15) now becomes where V F (ψ) is obtained from the first term in (2.16) using (2.17) and We emphasize again that throughout this paper we will, unless otherwise specified, set the reduced Planck mass to unity. Hence, all terms in Lagrangian density (2.18) are dimensionless -including the coordinates t and x i , i = 1, 2, 3.

Inflation and the primordial parameters
To analyze inflation, one chooses the metric to be in the Friedman-Robinson-Walker form and defines the Hubble parameter as H =ȧ/a. Then the complete set of dynamical equations associated with the canonically normalized gravity action (2.10) and the normalized ψ Lagrangian (2.18) is given by where the potential V is given by A careful discussion of the primordial parameters involved in inflationary cosmology, as well as their exact relationship to the potential energy function V (ψ) given in (2.24), was presented in [35]. The results were the following. Choosing the parameters  it was shown that the inflaton produced 60 e-foldings of inflation between the end of inflation at ψ end and a starting value of ψ * where ψ end 1.21 , ψ * 6.25 . (2.26) It followed from this that the spectral index n s , the tensor-to-scalar ratio r and the energy scale of inflation associated with the above parameters are given by respectively. These results are completely consistent with all of the Planck2015 data [54]. For the explicit parameters in (2.25), the graph of the complete potential energy associated with the inflationary epoch is shown in figure 1. Similarly, for the same choice of parameters, the plot of the Hubble parameter H between the time t end 9.89 × 10 6 at the end of inflation and the time t * = 0 when inflation begins is shown in figure 2. It is physically enlightening to restore mass units for both the mass of the inflaton as well as the µ coefficient. Then the results in (2.25) are given in by m = 1.58 × 10 13 GeV , Y ν3 ∼ 10 −12 , µ = 1.20 × 10 10 GeV (2.28) respectively. It follows that the scale of supersymmetry breaking -which is associated with m -is very large, whereas the µ-term is considerably smaller. Both scales play an important role in the analysis of Sneutrino-Higgs inflation and reheating. The parameters presented in (2.25), (2.28) are specific representatives of a wider class of parameters, all of which satisfy the Planck2015 data. Generically, one can choose any parameter m ∼ 10 13 GeV and Y ν3 ∼ 10 −12 . The allowed value of the parameter µ, however, is much wider. As was discussed in [35], its value has an upper bound given by µ = 1.20 × 10 10 GeV presented in (2.25), (2.28). However, any smaller value is allowed as long as it is consistent with the breaking of electroweak symmetry. For these values of µ, it was shown in [35] that V F V soft both during and after the inflationary epoch. That is, V F is only substantial prior to inflation.

The search for valid low energy points
In [35], we used the formalism presented in [30] to statistically search the space of initial soft supersymmetry breaking parameters for those points which 1) satisfy the Planck2015 data by using the parameters presented in (2.28) while 2) simultaneously being consistent with all present low energy phenomenological data -that is, appropriate B-L and EW breaking, all lower bounds on SUSY sparticles and the experimentally measured lightest neutral Higgs mass. We refer the reader to [30] for details of this formalism. Since the relevant phenomenological data is usually presented in GeV, we will work in this unit for this entire analysis, including the discussion of lowering the B-L scale presented in the following section. Suffice it here to say that initial dimensional soft SUSY breaking parameters are analyzed by randomly scattering all of them, with the exception of m Hu , in the interval [m/f, f m], where m = 1.58 × 10 13 GeV is the cosmologically consistent mass presented in (2.28) and f = 3.3. The parameter m Hu is then fixed, for each initial throw, by demanding that it satisfy the cosmological constraint given by (2.14) and (2.28). In this paper, to be consistent with the analysis below, we extend the results given in [31] by statistically throwing the initial parameters 50 million times instead of the 10 million times presented previously. The results satisfying both requirements 1) and 2) are shown as the "valid" black points in figure 3.
It is of interest to use the formalism presented in [30] to compute the B-L breaking scale associated with each of the 1406 valid black points presented in figure 3. This quantity JHEP09(2018)001 Hu is fixed by the cosmological constraint. We find that 4,209,300 points break B-L but not electroweak symmetry, and 860,084 points appropriately break both B-L and electroweak symmetry. Of the latter, 545,753 points are consistent with current LHC bounds on sparticle searches. Finally, we have 1406 points which satisfy all these conditions and are within the 2σ window of the measured Higgs mass. The black points are enlarged for legibility.The axes are two dominant parameters of the renormalization group equations and are defined in [30].
was not analyzed in any of the previous papers on the B-L MSSM and Sneutrino-Higgs inflation. However, as we will see below, knowledge of the B-L breaking scale is important in discussing reheating. Therefore, we have computed the B-L scales for all valid black points and present a statistical graph of the results in figure 4. For completeness, we have indicated the percentage of points for which the B-L scale M BL exceeds or is smaller than the supersymmetry breaking scale M SUSY defined in [30]. This is referred to as "rightside-up" and "upside-down" B-L breaking respectively. It is important to note that the smallest B-L scale associated with the valid black points is approximately 2 × 10 12 GeV. As we will see in the following, appropriate reheating will occur most naturally for values of M BL 10 12 GeV. It is important, therefore, to see if one can modify the initial statistical input of the soft breaking parameters so that one obtains physically realistic valid black points for which M BL is substantially smaller than 10 12 GeV. The answer is affirmative, as will be shown in the following section. discussed in [30] and used in the previous section. Rather, we will input the B-L scale as an arbitrary parameter as part of the initial data.

JHEP09(2018)001
To do this, we recall from [30] that the U(1) B−L symmetry breaks when the righthanded sneutrinoν c 3 obtains a non-vanishing vacuum expectation value. This occurs when the parameter m 2 JHEP09(2018)001 B-L breaking scale and use (3.3) and the relevant RGEs to determine the required soft sneutrino mass parameter at M U . Once this process is accomplished, we then have a complete set of initial data against which the phenomenological constraints 1), 2) and 3) above can be verified. To carry this out in detail relies heavily on a generalization of the formalism for the renormalization group equations of the B-L MSSM previously given in [30]. This is technically non-trivial and, hence, we present the mathematical details in appendix A of this paper. In this section, we will simply use the results obtained in that appendix.
As discussed in subsection A.2, one can choose a range over which one wants to input the B-L scale and then, using the formalism presented there, determine the phenomenologically acceptable valid black points whose B-L scales lie in that range. In the appendix, we carried this out over a very wide range of B-L scales, specifically from 10 6 GeV to 10 14 GeV. However, for the discussion of reheating in this paper, such a wide range for M BL is not required. Instead, we will limit our discussion in the text to B-L scales in the range 10 10 GeV ≤ M BL ≤ 10 12 GeV. Let us implement the procedure outlined in the appendix, now, however, for this restricted range of B-L scales. We will generate 50 million initial throws of the soft masses with the inputted scale of U(1) B−L breaking randomly generated from a log-uniform distribution between 10 10 GeV and 10 12 GeV. Carrying out our checks, we find that this ultimately leads to 215 sets of initial data which satisfy all phenomenological constraints. These physically valid black points are shown in figure 5. The distribution of the B-L breaking scale for the black points in figure 5 is given in figure 6.
Finally, using the formalism discussed in subsection A.3 of the appendix, we compute the amount of fine-tuning required to lower the B-L scale into the 10 10 GeV to 10 12 GeV range. The results for the 215 valid black points are shown in figure 7. Note that the JHEP09(2018)001

Post-inflationary epoch: classical behavior of ψ and H
In this section, we will begin our discussion of the post-inflationary epoch, assuming, for the time being, that the inflaton does not decay to normal matter. Within this context, we will calculate the classical behavior of the inflaton ψ and the Hubble parameter H numerically, and then present analytic solutions for both quantities which closely approximate the numerical results. Having achieved this, we will, in the next section, begin our discussion JHEP09(2018)001 of the details of the inflaton decay to normal matter and reheating. For this analysis, it is far more convenient to work in units where M P = 1.
In the inflationary and post-inflation epoch, the equations of motion for ψ and the Hubble parameter H are specified by equations (2.21), (2.22) and (2.23). The values for the parameters m, Y ν3 and µ will be chosen to be those given in (2.28) to ensure that the inflating epoch is consistent with all phenomenological data. As discussed above, for this choice of the Y ν3 and µ the F-term potential satisfies V F V soft in both the inflationary and post-inflation regimes. Therefore, in this section, we can, to a high degree of accuracy, simply take the potential energy to be V = V soft . Then the relevant equations of motion are given by These equations can be solved numerically for both ψ and H as functions of time. The results for ψ(t) and H(t) starting 1) at the beginning of inflation at t * = 0, 2) running through the inflationary epoch to t end 9.89 × 10 6 , and then 3) continuing into the postinflation epoch with for t > t end , are shown in figure 8(a) and (b) respectively. As is clear from figure 8(a), shortly after the end of the inflationary period, the inflaton will begin to oscillate around the minimum of its potential at ψ = 0. Taylor expanding V (ψ) in (4.4) around the origin, one obtains When ψ 3, V ≈ 1 2 m 2 ψ 2 and the mass of inflaton is m ψ = V ψψ ≈ m. Noting that m = 6.5 × 10 −6 , it follows from figure 8(b) that m H everywhere in the post-inflationary period. Hence, ψ will oscillate coherently around the minimum of V with a frequency ω = m ψ , although with decreasing amplitude. From figure 8(a), we can numerically show that the height of the first oscillatory peak corresponds to ψ 3 2 4.2 × 10 −3 , which easily satisfies the above criterion that ψ 3. Henceforth, for specificity, we consider this first peak as the beginning of the oscillatory phase and will denote the corresponding time, which we numerically compute, to be t osc 1.096 × 10 7 . This time is indicated by a dashed red line in figure 8(a). For all t > t osc we will, henceforth, take V (ψ) = 1 2 m 2 ψ 2 . During the oscillatory period, that is, when, t > t osc , equations (4.1), (4.2) and (4.3) then become ψ + 3Hψ + m 2 ψ = 0 .  Figure 8. The numerical solutions for ψ(t) and H(t), where we have set M P = 1. Note that t * = 0 and t end 9.89 × 10 6 mark the beginning and end of the inflationary period. The times t > t end correspond to the post inflationary epoch. As defined in the text, t osc 1.096 × 10 7 marks the time at which the potential energy is well approximated by V = 1 2 mψ 2 and t MD 1.387 × 10 7 is the time at which our analytic solutions for ψ and H become valid.

JHEP09(2018)001
We now want to find an approximate analytic solution to these equations for both H and ψ. To do this, we first neglect the effect of 3Hψ in (4.8). That is, we will take the lowest order approximation to H to be H 0 = 0. Then the solution of ψ to this order is nothing but an harmonic oscillator. That is where A 0 and c 1 are constants. Using (4.9), we have 1 2ψ Plug these expressions into (4.6) and (4.7), yields where H 1 is the first order approximation to H. It then follows that d dt (4.14) Integrating this from t osc to t(> t osc ) yields Clearly, for times t where t − t osc 1/m, we have where .
By numerically evaluating H(t osc ) using the results displayed in figure 8(b), we find that c 2 9.67 × 10 6 .
Having found an approximate analytic expression for H(t) beyond leading order, we now want to find the next order correction to ψ 0 given in (4.9). Due to the non-vanishing expansion of the Universe given by H ≈ H 1 , the amplitude of the oscillations of ψ will necessarily be damped. However, the frequency of the ψ oscillations will hardly be effected JHEP09(2018)001 as long as m H which, as mentioned previously, will be true for all t > t osc and, hence, for t − t osc 1/m. Thus we can set where A 1 (t) can be obtained by inserting expressions (4.16) and (4.18) into equation (4.8). This gives and, hence,Ä The solution is where A 1 is a constant. Putting (4.16), (4.18) and (4.22) into (4.6), we find that When t − c 2 1/m, which is automatically satisfied when t − t osc 1/m, we simply have Therefore, the next order analytic solution for ψ(t), valid in the region where t−t osc 1/m, is given by where c 2 9.67 × 10 6 was evaluated above. By matching (4.25) with the oscillations in the numerical solution of ψ, see figure 8(a), we can find that c 1 9.78 × 10 6 . For specificity, we note from the numerical calculation that the time associated with the fourth oscillatory peak in figure 8(a) is given by t 1.387 × 10 7 and satisfies t − t osc > 6π/m 1/m. Hence, to a high degree of approximation, the analytic solutions for H and ψ are both valid for any time larger than the time of the fourth oscillation peak. Since, as we will show below, this corresponds to the period of matter domination, we henceforth denote this time as t MD and indicated it by a red line in figure 8(a). To summarize, when t > t MD 1.387 × 10 7 where c 2 = t osc − 2 3H(tosc) 9.67 × 10 6 and c 1 9.78 × 10 6 . at the beginning and end of inflation, and at the beginning of both the oscillatory and matter dominated regimes respectively. We have set M p = 1. The numerical values of H and (ψ(t)/3) 2 at t * , t end , t osc and t MD are displayed in table 1. The regimes of inflation and matter domination are shown as the yellow and blue regions of figure 9 respectively. The duration of the intermediate phase, that is, the gray area in figure 9, is given by ∆t t MD − t end 3.97 × 10 6 . As will be shown below, this is negligible compared with the duration of the reheating period. For that reason, this "transition" regime will, henceforth, be ignored. Finally, as a check on our approximate analytic solution for H(t) in (4.26) and for ψ(t) in (4.27), we compare them in figure 10 (a) and (b) respectively against the exact numerical solutions for H and ψ in the region t > t MD . It is clear that our analytic solution is a very accurate approximation.

Post inflationary epoch: perturbative decay of ψ to matter
In the previous section, we ignored the quantum mechanical decay of the inflaton into various species of matter, focussing instead on its purely classical behavior and the associated classical behavior of the Hubble parameter. However, the ψ does decay into various species of matter, thus reheating the Universe. In this section we commence our discussion of these decays.

Dynamics of ψ and H during particle decay
Different decay processes can be occurring simultaneously, although they may have started at different times. In general, taking account of the decay of the inflaton, the ψ and H equations (4.6), (4.7) and (4.8) are modified to become where Γ d,i is the decay rate of ψ into the i-th matter species, and ρ i and p i are the energy density and pressure respectively of the i-th species in the decay products. The quantities ρ i and p i are related by the relation p i = ω i (t)ρ i , where ω i = 0 and 1/3 respectively for matter and radiation. The initial conditions for ψ and H are set by their classical values at the end of the inflationary epoch, and can be determined from the results in the previous section. Additionally, we have ρ i = 0 until the time at which the i-th decay JHEP09(2018)001 process commences; that is, when Γ d,i becomes non-zero. For convenience, we define the fraction of energy density of the i-th species as where, as follows from (5.1), the total energy density of the Universe is given by ρ total = 3H 2 . The fraction of energy density of the inflaton can be defined by In our theory, the inflaton is defined in (2.13) to be with the associated quadratic soft mass squared given in (2.14) by The value of m was fixed as m = 6.49 × 10 −6 (5.8) to be consistent with the results of Planck2015 [54]. The relationship of H 0 u ,ν L,3 andν c R,3 to φ 1 was presented in (2.12). Setting φ 2 and φ 3 to zero in those expressions -since their values vanish in the D-flat potential energy valley of the inflaton -gives Hence, each of these three fields can each be replaced by φ 1 in the Lagrangian density. However, as discussed above, φ 1 has a non-trivial Kähler potential and, hence, non-canonical kinetic energy. By performing the field redefinition in (2.17), that is we find that the ψ field is canonically normalized. We therefore used ψ in all our previous analysis. To analyze inflaton decay it is, therefore, essential that we re-express φ 1 in terms of ψ in the Lagrangian. Happily, expression (5.10) can be simplified in the post-inflationary regime. Taylor expanding (5.10) around ψ = 0, we find As can be seen from figure 9(b), in the matter dominated period t > t MD we find ψ 1. Hence, to a high degree of approximation, one can simply set in the Lagrangian. We do this in the following analysis. It follows that the decay of H 0 u , ν c R,3 andν L,3 can be viewed as the decay of the canonical scalar field ψ. As will be shown in detail, ψ is coupled with different classes of particles; including the standard model particles, charginos, nuetralinos, gauge bosons and scalar particles.

JHEP09(2018)001
A specific decay process can occur only when the total mass of the decay products is smaller than the mass of ψ. Since shortly after inflation, namely, when t > t osc , ψ will oscillate around the minimum of its potential V (ψ) = 1 2 m 2 ψ 2 , the mass of the inflaton is The mass of potential decay products can have two origins; namely, from soft mass terms in the Lagrangian or from the non-zero expectation value of the inflaton ψ. In fact, as well as inducing mass terms, the expectation value of the inflaton will also give rise to mixing terms between different fields. For example, the coupling 1 √ 2 g 2 ψW −ψ+ u , which arises from a super-covariant derivative term, gives rise to mixing between theW − gaugino and the Higgsino ψ + u . This will be discussed in more detail below. Given the essential role of the inflaton "expectation value" in determining the masses and couplings of its decay products, it is essential that this be carefully defined. Since ψ oscillates around 0 when t > t osc , its "naive" expectation value will vanish. However, it is clear that this is not the physical expectation value effecting the inflaton decay. Rather, we will use the root mean squared value of ψ, that is, ψ 2 , where ψ 2 can be defined by with δ being the period of the oscillations of ψ. As long as the total decay rate i Γ d,i m, which as shown below will always be satisfied, then the frequency of the oscillations of ψ can accurately be taken to be δ = 2π/m ψ . Having defined this root mean squared VEV for the inflaton, we will henceforth expand ψ as where δψ is a small fluctuation. As we will see below, using this expansion in the Lagrangian density will have two important ramifications; first, it will produce time dependent mass terms proportional to ψ 2 for each particle species and second, it will lead to a coupling of the inflaton fluctuation to matter -thus inducing quantum mechanical reheating of the universe. We will, for simplicity, often abuse notation and denote the fluctuation δψ simply as ψ. The correct meaning of the symbol will always be clear from the context.

Decay classes
In this subsection, we present a detailed analysis of the different types of matter into which the inflaton can decay. These are 1. up type standard model particles (ψ → tt, cc, uū); . These classes are distinguished by whether the decay products are fermions or bosons, and whether their masses must be determined by diagonalizing a mass matrix which arises due to mixing terms from the inflaton VEV. Figure 11. ψ → tt.

Up-type standard model fermions
Note from (2.13) that the inflaton contains H 0 u as a component field. It follows that the inflaton is able to decay via the Yukawa interactions directly into up-type standard model fermions. Since, as we will discuss below, the up-type leptons, that is, the neutrinos, can also mix with Higgsinos, we will treat these separately. Here, we consider only up-type quarks, since they cannot mix with other fermions. If we denote by F any of the u, c and t quarks then L ⊃ y HF HFF = y ψF ψFF , (5.16) where y HF is the usual Yukawa parameter for coupling to the Higgs and, using (5.9) and (5.12), the Yukawa parameter for coupling to ψ is The values for y HF depend on the energy scale at which they are evaluated, and can be determined at any given scale using the renormalization group analysis presented in [30]. As discussed in subsection 6.4, an appropriate scale in the interior of the reheating interval is 5.8 × 10 13 GeV. The values for y HF at this energy are found to be y Hu = 6.47 × 10 −6 , y Hc = 3.77 × 10 −3 , y Ht = 6.07 × 10 −1 .
That is, the inflaton can decay, in the order of the coupling strength, as ψ → tt, ψ → cc and ψ → uū. Consider the process ψ → tt as an example. Since ψ can decay into tt (see figure 11) where we define the four component Dirac spinors The decay rates of ψ to tt, cc and uū all have the following form. Noting that m F = mF , we find where the mass of the fermion is given by and m ψ = m = 6.49 × 10 −6 . It is important to note that the decay can only occur once 2m F < m ψ . Since m F is determined by ψ 2 , m F can initially be larger than m ψ /2. 2 In this case, Γ d = 0. With the expansion of the Universe, the amplitude of the oscillations of ψ will decrease. When ψ 2 becomes sufficiently small, the decay ψ → FF will become non-zero at some specific time, which we denote by t F * . For t > t F * , Γ d will increase as ψ 2 continues to get smaller. Eventually, when m F m ψ , it follows from (5.21) that Γ d will approach a constant. That is, Therefore, a species with a smaller Yukawa coupling constant will be produced earlier, but with a relatively smaller maximal decay rate than a species with a larger Yukawa coupling. The equation of state for FF is given by for 2m F ≤ m ψ . When 2m F m ψ , the decay products F andF are highly non-relativistic and ω FF ≈ 0. However, when 2m F m ψ , F andF are relativistic and, hence, ω FF ≈ 1/3. In the regime where Γ d Γ max d H, it is possible to give an approximate analytic solution for ρ FF using equations (5.1)-(5.4). However, this condition can only be satisfied for the up and charm quarks since their Yukawa parameters are relatively small. For the top quark, its Yukawa parameter is sufficiently large that Γ max d > H. Therefore, the results in the remainder of this subsection apply to u and c quark decays only. The density function ρ tt for the top quark can only be computed numerically. This will be carried out in section 6.
To lowest order, one can ignore ρ FF and, hence, p FF in (5.1) and (5.2), as well as Γ d in (5.3). It follows that H can still be approximated by as in (4.26), where c 2 9.67 × 10 6 . Putting this back into (5.3) and taking Γ d = Γ max d , one can solve this equation for ψ. We find that Note that as t → t F * , this expression for ρ FF → 0. That is, although derived in the in the regime where Γ d = Γ max d , we find that it remains a good approximation to the density for any t ≥ t F * since physically one knows that After t F * , that is, when 2m F < m ψ , ρ FF will initially increase with time and reach a maximum value of . Then, ρ FF will decrease with time as ρ FF ∼ (t − c 2 ) −1 . It follows from (5.5) and (5.29) that the fraction of energy density of species FF is given by It is interesting to note that for Γ max d H, by using the approximation in appendix B, we find that . (5.33) As we will determine below, the time at which reheating is finalized is given by t R 8×10 9 . It then follows from (5.26) and (5.32) that Ω uū (t R ) < 8.636 × 10 −9 , Ω cc (t R ) < 2.932 × 10 −3 . with t MD 1.387 × 10 7 . We conclude that although up-type fermions with small Yukawa coupling constants, that is, u and c, can be produced relatively early, their contribution to the background evolution of H and the final decay products of ψ are actually negligible. Physically, this is true because if Γ d H, the decay products will be diluted by the expansion of the Universe, thus barely effecting the evolution of H and ψ. As a proof of this, one can compare, for example, Ω cc (t R ) < 2.932 × 10 −3 against the smallest Ω(t R )

JHEP09(2018)001
computed numerically in section 6. This is found to be Ω BB (t R ) = 4 × 10 −3 . Noting that the value of Ω cc (t R ) is actually dramatically reduced relative to its value in (5.34) by the decay of the inflaton into the other species discussed below, we conclude that reheating into charm quarks, and therefore, into up quarks is negligibly small. They will, henceforth, be ignored. Only when a Yukawa parameter is large enough that the decay rate becomes comparable and then larger than H, will that species play an important role in reheating. As we will see below, this will be the case for the top quark.

Charginos
As mentioned previously, the non-zero expectation value of the inflaton -more precisely, the RMS value ψ 2 -gives rise to effective mass terms for fields, as well as to mixing between different particle species. By diagonalizing the mass matrix for such species, one can determine the correct mass eigenstates into which the inflaton decays. We now examine the first class of such mass eigenstates, which we will label "charginos", in direct analogy with the mass eigenstates associated with dynamical electroweak symmetry breaking in the B-L MSSM [30,62].
The mixing terms can arise from two sources; 1) the superpotential and 2) the"supercovariant derivative" of the H u Higgs doublet superfield. The former set are parameterized by y Hi ψ 2 , while the latter have the parameters g a ψ 2 , where y Hi , g a denote Higgs coupled Yukawa parameters and gauge couplings respectively. We give an explicit description of where these terms arise from in appendix C. Since the third family Yukawa coupling parameters are the largest, we will, for simplicity, assume that 1. All Yukawa coupling matrices are diagonal.
2. Only the third family quark and lepton Yukawa coupling parameters need be considered.
3. Since the third family neutrino Yukawa coupling parameter is also negligible, it can be dropped as well.
Dropping all terms which have a neutrino Yukawa coupling y ν and examining the effective mass Lagrangian for the "charginos", we find that one set of fields which are mixed due to the non-zero value of ψ 2 areW In order to construct the inflaton potential given in section 2, we have previously taken the supersymmetric µ parameter to be of O(10 10 GeV). This value is much smaller than the soft masses of the W -gauginos, as well as the initial values of the mixing terms y Hτ ψ 2 and g 2 ψ 2 . We can, therefore, simplify this system further by working in the "small µ" limit, and drop terms involving µ. Of course, as the value of ψ 2 decreases, the value of µ will eventually exceed that of other terms we have not dropped. However, this effect is not significant since it will only occur very near the end of the reheating period. Hence we can, to a good approximation, take µ to be negligible.

JHEP09(2018)001
In this limit, we are able to decouple the τ c R , ψ − d states since there is no longer any mixing between ψ − d and ψ + u . Examining the effective mass Lagrangian in equation (C.8), we see that where, using the formalism presented in [30], we find that the value of y Hτ at 5.8×10 13 GeV is given by y Hτ = 3.88 × 10 −2 . (5.38) Note that (5.37) is a mass term for a Dirac mass The decay rate to the ψ − d and τ c R states is, therefore, analogous to the decay of the inflaton to top quarks, and takes the form The equation of state parameter for ψ d τ is The remaining fermions,W remain mixed and form the new mass eigenstatesC ± 1 andC ± 2 , where Determining the matrices U , V is straightforward. The explicit form of both are given in appendix C.1. The statesC ± 1 andC ± 2 form Dirac fermions with masses mC 1 and mC 2 respectively. These are also presented in appendix C.1. We find that the large value of mC 2 makes the decay of the inflaton toC ± 2 kinematically impossible. Hence, onlyC ± 1 are produced. The decay rate of the inflaton to the mass eigenstatesC ± 1 is then given by The elements of U, V are given in appendix C.1. The equation of state parameter forC ± 1 is It is important to note that the rate and equation of state for the inflaton decay intõ C ± 1 depend on the soft SU(2) L gaugino mass M 2 . However, this will vary statistically over the interval [m/f, f m], where m = 1.58 × 10 13 GeV and f = 3.3. Generically, it will be different for each of the 215 valid black points discussed in section 3. To avoid having to do a separate analysis for each of these 215 black points, we will, instead, note that one expects their average value, denoted by M , to be near the center of the interval. Furthermore, for concreteness, we will henceforth assume that M = m = 1.58 × 10 13 GeV . (5.48) Looking ahead, we note that inflaton decays into different species will depend on the gaugino soft masses M R and M B−L , as well as on M 2 . Therefore, for concreteness, we will henceforth make the generic assumption that Secondly, we note that the rate and equation of state for the inflaton decay intoC ± 1 also depend on the SU(2) L gauge parameter g 2 . This quantity is evaluated, as are all the other gauge couplings, by running it from its measured value at the electroweak scale up to the scale of reheating at ∼ 5.8 × 10 13 GeV. Hence, its value will essentially be the same for all 215 valid black points. Using the formalism developed in [30], we find that at 5.8×10 13 GeV g 2 = 0.57 . (5.50) Again, looking ahead we find that inflaton decays into different species will depend on the gauge couplings g R and g BL as well as on g 2 , all evaluated at the reheating scale of 5.8 × 10 13 GeV. Using the formalism developed in [30], we find that at 5.8 × 10 13 GeV g 3 = 0.60 , g 2 = 0.57 , g R = 0.56 , g BL = 0.56 (5.51) with an average value of g = 0.57. Therefore, for simplicity of calculation, we will henceforth make the generic assumption that As with the top quark, the expressions for the energy densities ρ ψτ and the ρCC charginos can only be computed numerically. We will carry this out in section 6.

Neutralinos
We now turn to the second set of particles which mix due to the non-zero inflaton VEV. We refer to these as "neutralinos", in analogy with states described by the B-L MSSM. Again, we will ignore all terms multiplied by a neutrino Yukawa coupling parameter y Hν . Making the same assumptions about Yukawa couplings as before, the effective mass Lagrangian for the "neutralinos" is such that the states that mix arẽ Once again, we take µ to be negligible compared to other terms in the effective mass Lagrangian. In this limit, the state ψ 0 d decouples. It follows that the effective mixed state mass matrix, MÑ , is six-by-six. To proceed, this must be diagonalized. We leave the details of this to appendix C, but note that to simplify our expressions we again make the assumptions given in (5.49) and (5.52). Additionally, we define Diagonalizing the mass matrix MÑ , we find the six mass eigenstates given in table 2. Of the six eigenstates, onlyÑ 1 andÑ 2a ,Ñ 2b are kinematically accessible to the decay of the inflaton. The decay processes and rates are SinceÑ 2a andÑ 2b have the same mass, their equations of states parameters are given by the same form As with the top quark and the charginos, the expressions for the energy densities ρÑÑ for the neutralinos can only be computed numerically. We will carry this out in section 6.

JHEP09(2018)001
Mass Degeneracy State Table 2. Mass eigenstates of the neutralino mass matrix MÑ . The masses M and u are defined in (5.48) and (5.54) respectively.

Gauge bosons
The covariant derivatives of H 0 u ,ν 3L andν c 3R couple the inflaton ψ to the associated gauge bosons and, furthermore, give an effective mass to these bosons. This occurs via the following terms in the Lagrangian To find the mass for each species of vector boson, as well as to determine their coupling parameter to the inflaton, we expand the inflaton around its root mean squared VEV as in (5.15). Inserting this into the final expression in (5.58) and, as previously discussed, denoting δψ simply as ψ, we find that where the masses are given by For a generic coupling G i ψW µ i W iµ with identical W-bosons, the decay rate and the equation of state parameter are given by ) where the gauge boson masses are given in (5.60).
As with the top quark, charginos and neutralinos, the expressions for the energy densities ρ W W and ρ BB for the gauge fields can only be computed numerically. We will carry this out in section 6. To simplicity the calculations, we will again use the approximation for the gauge couplings presented in (5.52).

Scalars
The inflaton can couple to other scalar fields via the supersymmetric F-term and D-term potentials, as well as the soft supersymmetry breaking terms. These couplings give rise to a potential three-body decay vertex, as well as mass terms and mixing terms. To find out the mass eigenstates into which the inflaton can decay, one must examine the mass matrix M whose elements are where i, j run over all scalars other than the inflaton and

JHEP09(2018)001
To simplify our calculations, we assume that all scalars have identical soft masses given by m = 1.58 × 10 13 GeV (5.73) and take the gauge couplings to have their average value at order 5.8 × 10 13 GeV, as discussed above. That is, Diagonalizing M, we find that the eigenvalues of M are either m 2 or larger. It follows that there are no decays of the inflaton to scalars. Additionally, we note that as we are describing a supersymmetric theory of inflation, it is also natural to consider the spin-3 2 superpartner of the graviton, the gravitino. The over-production of gravitinos during reheating can potentially be hazardous to nucleosynthesis [63][64][65]. However, in our present discussion, we expect that the mass of the gravitino is of order m ∼ 10 13 GeV (that is, the same order as the soft supersymmetry breaking scale) and hence its production from the decay of the inflaton is heavily suppressed.

Numerical calculation
As discussed in previous sections, the inflaton can decay into different species with different time dependent decay rates. The root mean squared value of ψ, that is, ψ 2 , decreases with the decrease of the oscillatory amplitude of ψ due to the expansion of the Universe and the decay of the inflaton. Thus, different decay processes will begin at different times since the masses of the decay products depend on ψ 2 . We plot the decay rates for different processes with respect to ψ 2 in figure 12. We did not plot the decay rates for uū and cc because, as discussed above, they are too small to have any substantial effect. The values of ψ 2 at which each relevant process is turned on can be found in table 3. Even though we have simplified the computations by ignoring the u and c quark decays, it remains impossible to find analytical solutions for (5.1)-(5.4) to account for all the relevant decay processes simultaneously. Therefore, in this section, we will numerically solve (5.1)-(5.4) to find the solutions for H, ψ and ρ i . From this, one can determine the relative energy densities of the different species at the end of the reheating epoch, as well as the reheating temperature.

Initial conditions
To find the solutions for H, ψ and ρ i by numerically solving (5.1)-(5.4), one needs the initial conditions for H, ψ and ρ i . In principle, one can solve these equations starting from the beginning of inflation. Thus, the initial conditions would be set by inflation. However, such an approach would take a great deal of computing time due to the severe oscillations of ψ after t osc . Furthermore, ignoring the u and c quark decays, it follows from table 3 that the first decay process to turn on is ψ → ψ − d τ c R . The time at which this decay commences can be computed from the expression , (6.1)

Decay processes
Value of ψ 2 Value of ψ 2 at turn on (M p = 1) at turn on (GeV)  Table 3. Values of ψ 2 at which each decay process is turned on. We use the Yukawa couplings evaluated at the reheating scale of order 5.8 × 10 13 GeV.
where t osc 1.096 × 10 7 from figure 8, H(t osc ) 5.16 × 10 −7 from table 1 and y Hτ was given in (5.38). This expression was first presented in (5.33) for top-quark decays, but can be used here since both u and c are being neglected. The result is which is much later than t osc , thus exacerbating the computing time even more. Therefore, to save computing time, we will start our calculation from an arbitrarily chosen time JHEP09(2018)001 t = t I , where t I is close to, but smaller than, t ψ − d τ c R * . The exact value of t I is for technical convenience only. We will set t I = 8 × 10 8 . When t ≤ t ψ − d τ c R * , we can neglect all decay effects. Thus, we can approximately set ρ i = 0 at t I for all decay products. Using (4.27) and (4.26), one can obtain the initial conditions for ψ and H at t I , respectively. Of course, the decay rates for each process are zero until the corresponding process is turned on.
As can be ascertained from figure 12 and table 3, the second decay to turn on is ψ →C + 1C − 1 . It is clear from this data that the associated time, tC 1 * , will be much later than t ψ − d τ c R * . As we will see below, tC 1 * 6.45 × 10 9 .
Therefore, to further reduce the time for computation, we will divide the numerical calculations into two parts. First, we will numerically compute from t = t I to some time t = t II < tC 1 * by using the iterative method described in next subsection. Again, the choice of the value of t II is for technical convenience only, which will not make any physical difference as long as t II < tC 1 * . We will choose t II = 5 × 10 9 . Second, we numerically compute, also using the iterative method, from t = t II to the time t R where reheating has been completed. The initial conditions for ψ, H and ρ i at t = t II will be set by the numerical solutions of the first part; that is, for t I < t < t II .

Iterative method
In (5.1)-(5.4), the equation of state parameters ω i and the decay rates Γ d,i depend on the root mean square value of ψ, that is, ψ 2 . This makes these background equations very difficult to solve -even numerically -especially when the ψ 2 -dependence becomes complicate after t I . Hence, we will solve them by iteration. We accomplish this by using eq. (4.27), which is the solution for ψ without considering any decay, as the first input ψ for ψ 2 in ω i and Γ d,i . Then, we treat ω i and Γ d,i simply as some time-dependent functions so that we can find solutions for (5.1)-(5.4) numerically. This gives the first output ψ(t). The first output ψ(t) will be identical with the first input ψ(t) until some observable decay processes are turned on. Once i Γ d,i becomes effectively nonzero, the oscillations of ψ will be damped more quickly. From then on, the output ψ 2 will be smaller than that of the input.
Next, we use this first output ψ(t) as a second input ψ for ψ 2 in ω i and Γ d,i and numerically solve (5.1)-(5.4) again. This then leads to the second output ψ(t) that is closer to the final solution. By repeating this method, the output ψ(t) will become closer and closer to the real solution of ψ. We repeat this method iteratively multiple times until the output ψ(t) is almost identical with the input ψ(t), which means we have found the correct solution for ψ. Using this method also leads to solutions for H and ρ i to some reasonable accuracy.
Additionally, since the inflaton will oscillate rapidly around the minimum of its potential during the reheating phase, the oscillations of ψ are too dense to be easily handled in the numerical calculation. The amplitude of the oscillations of ψ will decrease with time. However, the frequency of the oscillations is almost a constant as long as i Γ d,i m ψ , as we will demonstrate in the numerical results. Hence, as a good approximation we can set

JHEP09(2018)001
where c 1 9.78 × 10 6 . This is very useful in getting rid of the obstacles described above. Then, by using the discussion in appendix B, we can simply replace ψ 2 in ω i and Γ d,i with A(t)/ √ 2. Therefore, for every input ψ(t) in the iterative calculation, we can focus on A(t) instead of ψ(t). We apply this iterative method to both the first part of the calculation, where t I < t < t II , and to the second part, where t > t II , which we respectively denote as part I and part II. Eventually, when the input A(t) and the output oscillatory amplitude are almost same, we can conclude that we have found the correct solution for A(t) (or equivalently ψ(t)), H(t) (or equivalently a(t)) and ρ i (t) (or equivalently Ω i (t)).
The possible corrections to i Ω i from the accuracy of the solution of ψ can be defined as where ρ ψin = 1/2ψ 2 in + V (ψ in ), ρ ψout = 1/2ψ 2 out + V (ψ out ) and ψ in , ψ out are the input and the output ψ, respectively. Note that eq. (6.4) should be used when we transform between ψ and A.

Numerical results
After several iterations, we obtained the final solution for A(t) (or, equivalently, for ψ(t)). For the last round of n such iterations, we plot the input A(t) (let us denote it by A n (t)), the corresponding output ψ(t) or A(t) (let us denote it by ψ n+1 (t) or A n+1 (t), respectively) and the solution of H in figure 13. A n (t) is actually the numerical output A(t) of the (n − 1)th round of iteration. Note that the oscillations of ψ are too dense to be plotted completely. Instead, we simply plot 5000 random points from the curve of ψ n+1 (t) for both part I and part II. Thus A n+1 corresponds to the upper boundaries of the magenta points in figure 13, while A n is the black curves in figure 13. We can see that A n and A n+1 are almost same. We can, therefore, treat ψ n+1 as the actual solution of ψ. We have verified that their deviation is sufficiently small for our interests in this paper. 3 We plot the decay rates Γ d,i and the Hubble parameter H in figure 14. Obviously, i Γ d,i m ψ = 6.49 × 10 −6 (= 1.58 × 10 13 GeV) throughout. Thus the decay of ψ cannot significantly effect its oscillatory frequency. As can be seen in figure 14, the decay process ψ → ψ − d τ c R turns on much earlier than the other species, as was quantified above. Its decay rate reaches its maximal value and then becomes a constant. This is comparable to, but smaller than, H prior to the other species in figure 14 turning on. Thus the backreaction from ψ → ψ − d τ c R cannot be neglected. The decay rate of ψ → tt is similar to that of ψ → ψ − d τ c R , but turns on later and with a much larger maximal value. The decay rates of other species in figure 14 first increase with time after they are turned on. However, since they are proportional to ψ 2 , when ψ 2 is small enough they achieve a maximum and then decrease with time. Note that Γ(ψ →Ñ 2Ñ2 ) is the total decay rate for the three processes ψ →Ñ 2aÑ2a ,Ñ 2bÑ2b ,Ñ 2aÑ2b .
The evolutions of Ω i and Ω ψ are displayed in figure 15. Ω N 2 N 2 includes the decay products for all three processes ψ →Ñ 2aÑ2a ,Ñ 2bÑ2b ,Ñ 2aÑ2b . We did not specify them   Figure 13. In both (a) and (b), the 5000 magenta points are randomly chosen from the curve of the output ψ(t), in the last (nth) round of iteration, during the time intervals t I ≤ t ≤ t II and t ≥ t II respectively. Their upper boundaries are identical to the input A(t) in the last round of iteration, that is, A n , which is the black curves in (a) and (b). In (c), we plot the upper boundaries of the magenta points, that is, A n+1 , from t I to sometime after reheating. In addition, the time at which each specie is turned on is marked with a vertical line, where t τ ψ * = 8.78 × 10 8 , tCC * = 6.45 × 10 9 , t W W * = 6.93 × 10 9 , t tt * = 6.95 × 10 9 , t BB * = 7.06 × 10 9 and t N N * = 7.08 × 10 9 . In (d), we plot the solution for H in the last round of iteration, that is, H num , and also (4.26) as a reference, from t I to sometime after reheating. We always set t I = 8 × 10 8 and t II = 5 × 10 9 in our numerical calculations for technical convenience and everywhere set M P = 1.

-15
10 -12 Figure 14. In (a), we plot the evolutions of H and the decay rates Γ d,i from t I to sometime after the end of reheating. In (b), we plot H and Γ d,i from t II to sometime after the end of reheating. Note that Γ d (ψ → ψ − d τ c R ) is very close to H when t > t II . We set t I = 8 × 10 8 and t II = 5 × 10 9 . In both (a) and (b), Γ(ψ →Ñ 2Ñ2 ) is the total decay rate for three processes ψ →Ñ 2aÑ2a ,Ñ 2bÑ2b , N 2aÑ2b and we have set M P = 1.

JHEP09(2018)001
When t t R , we find H(t R ) 7.8 × 10 −11 ( 1.9 × 10 8 GeV). Since the Universe is now dominated by relativistic particle species, that is, 3M 2 p H 2 = ρ rel , one can, assuming the Universe is in thermal equilibrium, find the reheating temperature from the expression where g * is the (effectively) massless degrees of freedom and T R is the temperature. It follows that the reheating temperature for the Sneutrino-Higgs theory is Since reheating takes place to a mixture of standard model particles (such as the top quark and various W -bosons) and lighter supersymmetric sparticles (such asC ± 1 ), the counting of the degrees of freedom is complicated. However, all of these species will eventually decay to the standard model particles with right-handed neutrinos, which has g * = 118. Hence, it is sufficient for our purposes to make a crude approximation and take this as the number of degrees of freedom at the reheating temperature. It follows that T R 1.13 × 10 13 GeV . (6.10) In this paper, we are requiring that the B-L breaking scale be well separated from the scale at which reheating takes place; that is B-L breaking occurs at a scale 10 13 GeV. This simplifies the reheating calculations and, more importantly, allows reheating to occur prior to the breaking of baryon and lepton number. As discussed in section 3 and appendix A, the B-L scale can be made arbitrarily small, albeit at the expense of some fine-tuning. Clearly, the above requirement will be fulfilled for the B-L scales between 10 10 GeV and 10 12 GeV discussed in section 3. As a concrete example, we see from figure 6 that of the 215 phenomenologically valid black points, the maximal number (≈ 20) occur at a B-L scale of 10 11 GeV. Henceforth, as an example, let us choose this to be the B-L scale. It follows that the associated energy density is ρ BL = 10 44 GeV 4 . Hence Thus H(t R ) H(t BL ), which indicates that B-L breaking will occur much later than the end of the reheating epoch.

The reheating interval
The entire analysis of reheating discussed above depends on inputting the numerical values of specific Yukawa parameters and all of the gauge couplings. However, these quantities all "run" with energy scale, changing their values to satisfy the renormalization group equations. In this paper, we have made two important assumptions that drastically simplify, and clarify, the calculations of reheating. 1) As we will see below, the interval of active reheating to matter is less than one order of magnitude in GeV units. Hence, these parameters vary only minimally over this small energy range. It is therefore a good approximation

JHEP09(2018)001
to choose a point in the interior of the reheating interval and to evaluate all Yukawa and gauge couplings at this fixed scale. 2) We find a point in the interior of the reheating interval as follows. We choose an arbitrary energy scale within an order of magnitude of where we expect to find the end of reheating. Using all Yukawa and gauge couplings evaluated at this energy, we numerically compute t ψ − d τ c R * and t R and the Hubble parameters associated with them. We then take the average of the Hubble parameters, H avg , and convert it to a "matter" energy density using 3H 2 M 2 P = ρ avg . The interior energy scale is then chosen to be the fourth root of ρ avg . Further iteration shows that this interior point remains a good approximation for characterizing the reheating interval.
Specifically, in this paper we do the following. First consider 2). We begin by choosing the arbitrary scale to be 10 12 GeV and use the RGE's to compute all Yukawa and gauge parameters at this energy. The numerical calculation of Γ i , Ω i and H, as well as t ψ − d τ c R * and t R , are carried out as described previously in this section. We find that the associated Hubble parameters are GeV , H(t R ) 1.6 × 10 8 GeV It then follows from assumption 1) that all Yukawa and coupling parameters used in our reheating calculations will be evaluated at 5.8 × 10 13 GeV.

Attaining equilibrium
In order to define a reheating temperature for the plasma of decay products, one must determine that they have attained equilibrium [66]. In this section, we will show that this is indeed the case prior to t R 8 × 10 9 . Thermal equilibrium occurs when the interaction rate of the i-th decay product, which we denote by Γ i int , is sufficiently large for all species i. Specifically, one requires that This implies that the mean interaction length, 1/Γ i int , is within the causal horizon 1/H. To demonstrate that (7.1) is indeed satisfied by the end of reheating, let us consider the elastic scattering of the charginos,C ± 1 , mediated by gauge bosons. As we have shown in figure 15(b), these comprise the largest fraction of inflaton decay products by the end of reheating. We argue that all other interaction processes involving different species present in the plasma will have similar interaction rates -since these will also involve gauge boson mediated scattering, all with similarly large values of gauge couplings. Therefore, if condition (7.1) is satisfied for one species, it is safe to assume that it is satisfied for all of them by the time that reheating is completed.
For simplicity, let us determine the rate for the processC + where, for simplicity, we take the mediating gauge boson to be W 0 µ . This process is also mediated by W Rµ and B µ , but since at this energy scale the gauge couplings are all of similar value, JHEP09(2018)001 see (5.74), the results will be very similar. Feynman diagrams which contribute to this process are shown in figure 16. The differential cross-section for this interaction is where the coupling O a is given by 3) The matrix elements V 1u , V 1W are given in appendix C.1. We have used the fact that in this calculation. Here, as in the rest of this paper, we follow the conventions and notation outlined in [67]. In the center-of-mass frame, the Mandelstam variables s, t and u are where we have assumed that the incoming states carry energy E cm /2, with In deriving equation (7.2), we have ignored the mass of the gauge boson, which quickly becomes negligible as the inflaton VEV decreases. The total cross section is then given by integrating

JHEP09(2018)001
where the cutoff δ must be introduced to remove the collinear divergence. We use the standard result that In figure 17, we plot the resulting cross-section as a function of the parameter x 2 = g 2 ψ 2 / √ 6, which is proportional to the inflaton expectation value. The interaction rate is given by where n is the number density ofC ± 1 and v is the average particle velocity, which we take to be v ∼ c = 1. For a given particle species with average energy E , the number density can be approximated by the expression ) is plotted as a function of time along with the Hubble parameter in figure 18. It is clear that condition (7.1) is satisfied well before the completion of reheating. Since the self-scattering interactions of other decay products involve similarly sized gauge couplings, we expect analogous rates for these species -for example, top quarks and W bosons -to also satisfy condtion (7.1) by the end of reheating, despite their smaller number density.
Strictly speaking, we have determined that the particles above have attained kinetic equilibrium, one of two necessary conditions to ensure that the decay products of the inflaton have thermalized [66,68,69]. The second condition, the achievement of chemical equilibrium, requires the analysis of number-violating 2 → 3 interactions. An example of such a process is given in figure 19. Naively, these interactions may be suppressed by an extra factors of perturbative couplings and hence their rates may be smaller than the number conserving 2 → 2 scattering -although this is not always the case. Without going into the full details of such processes, we will simply argue that their rates are still sufficiently high. That is, a 2 → 3 scattering rate involving the charginos could at worse be Under this assumption, from figure 18, we can see that the condition is also attained before the end of reheating. We expect that this holds for all other decay products as well.  Figure 19. Example of a Feynman diagram for 2 → 3 inelastic scattering involving the charginos.

Acknowledgments
We and m 2 Hu (t U ). To do this, we rely heavily on the description of the RGE's of the B-L MSSM presented in [30], as well as on the phenomenological analyses of the low-energy physics discussed in those papers.

A.1 RG running of the sneutrino mass
In the upside-down hierarchy, the RGEs describing the running of m 2 where we define t = ln( M M Z ). Here, M a , g a are the associated gaugino masses and gauge couplings for a = BL, 3R, and the S a -terms are defined as

JHEP09(2018)001
The RGEs for the M a and S a are given in [30], and can be solved to yield the following analytic expressions.
The running of the gauge couplings g a can be found by solving the RGEs with the boundary conditions The at a given scale. We find that for t SUSY < t < t U , and Given that we know m 2 ν c R,3 (t) at t = t BL , we would like to determine the corresponding value of m 2 ν c R,3 (t) at t = t U . Rearranging equations (A.9) and (A.10), we find that (A.11) -44 -

JHEP09(2018)001
However, we must now use the fact that S B−L (t U ) and S R (t U ) contain m 2 ν c R,3 (t U ). To accomplish this, we redefine the above expressions in terms of two new objects, S B−L (t U ) and S R (t U ) , given by respectively. For simplicity, we have suppressed the fact that all terms in (A.12) and (A.13) are evaluated at t = t U . Doing this, and re-arranging terms to extract m 2 ν c R,3 (t U ), we arrive at the expression (A.14)

A.2 Cosmological constraint
Recall from (2.14) and (2.25) that, in order to construct our inflaton and match the necessary cosmological data, it is required that We have previously satisfied this condition in our computational framework by not randomly generating m 2 Hu (t U ) but, instead, randomly generating the other two masses and using (A.15) to calculate the required value of m 2 Hu (t U ). In the context of our present discussion, we must take this expression into account when we enforce a specific B-L scale. Since this relation involves mν c 3 , and m Hu enters (A.14) via the S-terms, we can see that equations (A.14) and (A.15) are intertwined and must be solved simultaneously. To do this, we can re-express these equations in the form The equations in (A.16) are then soluble by inverting the matrix expression to give and inserting A, B and C in (A.17). In particular, for m 2

JHEP09(2018)001
A.3 Fine-tuning As we have seen, in order to ensure that the scale of B-L breaking have a desired value, the mass of the third family right-handed sneutrino has to be adjusted against the other masses. To quantify the degree of fine-tuning necessary, we re-examine equation (A.14).
Re-arranging this equation as an expression for M BL , we find This can be schematically expressed as .

Equation (A.24)
shows us precisely where the delicate cancellation necessary to produce a low B-L breaking scale arises. To quantify the degree of fine-tuning, we can plot the ratio F X

B Useful approximation
Consider the quantity where δ = 2π/m ψ is the period of the oscillations of ψ. In order to remove of the integral in this expression, we use the approximation that A(t) does not change much in the time interval t − δ to t + δ. This is true since δ H −1 as long as m ψ H. It follows that

C Diagonalization of mass matrices
In the B-L MSSM we have three families of quark and lepton chiral superfields where we assume that the Yukawa parameters are family diagonal and real, and do not display the family index. This gives rise to the Lagrangian as well as the purely scalar F -term potental We also have the soft-mass terms for the SU(2), U(1) 3R and U(1) B−L gauginos given by Additionally, the following terms arise from the gauge super-covariant derivative: The inflaton, φ = (H 0 u +ν L,3 +ν c R,3 )/ √ 3, attains a time-dependent expectation value during both the inflationary and post-inflationary periods. In doing so, it gives rise to fermion mixing terms in the Lagrangian. We now determine the mass eigenstates and eigenvalues due to this mixing.

C.1 Chargino mixing
Dropping terms which have any neutrino Yukawa coupling y ν , the effective mass Lagrangian for the "charginos" is given by

JHEP09(2018)001
Of course, the expectation values need to be expressed in terms of the RMS value of the inflaton field using However, for clarity, we will continue to express the expectation values in terms of the original fields until it becomes necessary to evaluate them. We can re-express (C.8) as and The mass eigenstates of this system can be found by diagonalizing the 3-by-3 matrix X, using the two 3-by-3 unitary matrices U and V defined by (C.14) It is easier to find the matrices U and V from the expressions As it stands, the expressions one gets for the mass eigenvalues and mixing matrices U, V involve solving for the roots of a cubic equation and give very cumbersome expressions. We can simplify the situation by working in the limit in which the µ-term is negligible and can be dropped -a reasonable approximation during the period of reheating, where ψ 2 is initally around 10 14 GeV and the gaugino mass M 2 is always of O(10 13 GeV). As ψ 2 decreases, the µ-term does become comparable to (and indeed larger than) terms such as g 2 ψ 2 that we have kept. However, this will occur after the bulk of reheating has occurred, and thus will have an insignificant effect, so we will simply drop µ in our subsequent calculations.

JHEP09(2018)001
Decays of the inflaton into the charginos arise from the vertices Assuming that the states with mass mC 2 are too heavy to be produced, we will be interested in the decays toC + 1 ,C − 1 only. Rotating the vertices above, we find that (C.28) For later use, we note that in the limit that x 2 x 1 , the above matrix elements have the form x 2 x 1 2 + 11 8 (C.29) The decay rate for the process ψ →C + |U 1W V 1u + U 1τ V 1W | 2 . (C.31) Since α, β are real, we can simplify this further to get Here, the coupling

C.2 Neutralino mixing
Again, ignoring terms which come with a factor of a neutrino Yukawa coupling y ν , the effective mass Lagrangian for the "neutralinos" is given by Here, q Gx denotes the U(1) G charge of the field x = H u , ν c R , L, where G = R, B − L. We can express this as L mass,N = − 1 2 ψ T M N ψ + h.c. , (C. 36) where ψ T = B ,W R ,W 0 , ψ 0 d , ψ 0 u , ν L,3 , ν c R,3 (C. 37) and C.2.1 The small µ limit In this limit, the ψ 0 d state decouples from the mixing. Hence, we only have a six-dimensional system to analyze. This is given by Let us make some further simplifying assumptions in order to find the eigenvalues and eigenvectors of this system.
• Take −x 2 = x 3 = x 5 = −x 6 = −x 8 = x 9 = u ∼ g ψ 2 , where g ∼ 0.57 The mass matrix N then takes the form The eigenvalues and eigenvectors of this system can now be evaluated and are given in table 2 in the main body of this paper.

JHEP09(2018)001
We find that the lightest eigenstates, that is, those into which the inflaton can decay, are given in terms of the states in (C.40) bỹ

JHEP09(2018)001
Rotating the Lagrangian to the lightest mass eigenstatesÑ 1 ,Ñ 2a andÑ 2b , we find that the term ψÑ 1Ñ1 is absent, and hence the decay ψ →Ñ 1Ñ1 is not allowed. Furthermore, we find that L decay,N ⊃ 1 √ 6 It follows that there is no decay of ψ toÑ 1Ñ2a or toÑ 1Ñ2b . Searching the Lagrangian for the remaining kinematically allowed decay terms, we find L decay,N ⊃ 1 √ 6 This allows the inflaton decays ψ →Ñ 2aÑ2a , ψ →Ñ 2bÑ2b and ψ →Ñ 2aÑ2b . Feynman diagrams for these processes are shown in figures 24 and 25. The first two decay rates take the same form where α = β = iO a and iO b for x = a, b respectively. The decay rate toÑ 2aÑ2b takes the form Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.