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 ${\cal{O}}(10^{13})~{\rm 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 $\simeq 1.13 \times 10^{13}~{\rm 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.


Introduction
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 loosing 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 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 F-term potential V F and contributions from the soft supersymmetry breaking terms V sof t . 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 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 sof t 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 reheating of the Sneutrino-Higgs theory. The subject of reheating in inflationary models has been studied extensively in [40,41,42,43,44,45,46,47,48,49,50]. 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 [51]. 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-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 GeV-and 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 [52]. 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 M D > 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 M D . 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 M D . 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 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 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 [55] 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 . 5 The formalism for coupling of a generic four-dimensional compactification of M-theory to N = 1 supergravity was presented in [56], and can easily be applied to the B-L MSSM . We note that the Planck mass used in this formalism is the socalled "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 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 where the D a , a = 3, 2, 3R, B − L functions are and T r (a) , r = 1, . . . , dim G a are the 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 (8), the soft supersymmetry breaking terms also include mass terms for each of the gauginos given by

−L
Setting all fields to zero with the exception of φ 1 , the V D potential vanishes and the Lagrangian becomes where 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) and (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 (15) now becomes where V F (ψ) is obtained from the first term in (16) using (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 (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 (10) and the normalized ψ Lagrangian (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 (24), was presented in [35]. The results were the following. Choosing the parameters m = 6.49 × 10 −6 , Y ν3 ∼ 10 −12 , µ = 4.93 × 10 −4 (25) 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 .
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 [51]. For the explicit parameters in (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 (25) are given in by m = 1.58 × 10 13 GeV , Y ν3 ∼ 10 −12 , µ = 1.20 × 10 10 GeV (28) Figure 1: The black line is a graph of the potential V sof t + VF for the parameters m = 6.49 × 10 −6 , Yν3 ∼ 10 −12 and µ = 4.93 × 10 −4 . For these values of the parameters, the vertical red dashed lines mark ψ end 1.21 and ψ * 6.25 respectively. We have set MP = 1. 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 (25), (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 (25), (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 sof t 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 the parameters presented in (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 (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 (14) and (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 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 SU SY defined in [30]. This is referred to as "right-side-up" and 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].
"upside-down" B-L breaking respectively. It is important to note that the smallest  Figure 3. The B-L and supersymmetry breaking scales are computed using the formalism presented in [30]. We indicate what fraction of each bin consists of those cases in which MBL > MSUSY (right-side-up), and in which MBL < MSUSY (upside-down). For example, between 1.0 and 1.1 ×10 13 GeV the number of right-side-up valid points is ≈ 120 whereas the number of upside-down valid points is ≈ 210 − 120 = 90.
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.

Lowering the B-L scale
We would like to arbitrarily lower the scale of B-L breaking for a given set of initial data so that the RGE evolution of this data is consistent with all phenomenological constraints; that is, 1) the electroweak scale is radiatively broken with the correct Z and W boson masses, 2) all sparticle masses exceed their present experimental lower bounds and 3) the Higgs mass is given to within 2σ of its measured value of 125 GeV. To accomplish this, we will no longer compute the scale of B-L breaking as an "output" of the initial class of data 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.
To do this, we recall from [30] that the U (1) B−L symmetry breaks when the right-handed sneutrinoν c 3 obtains a non-vanishing vacuum expectation value.
This occurs when the parameter m 2 ν c 3 turns negative as one runs down in energymomentum from the unification scale M U defined in [31]. Since natural reheating will require a lower value of M BL , we will assume that our B-L scale will be less than the scale of SUSY breaking; that is, we will only consider the "upsidedown" hierarchy where M BL < M SU SY . The B-L scale M BL is defined in [30] via the recursive relation where M Z R is the mass of the Z boson, which receives a mass when the righthanded sneutrino develops a vacuum expectation. We also have the relation It follows that where m 2 (M BL ) also be fixed. If we continue to randomly generate mν c 3 at the unification scale M U , as was done in all previous analyses, it is exceedingly improbable that, upon running the sneutrino mass down, we will arrive at our desired value of m 2 ν c 3 (M BL ) and, hence, of M BL . This problem is concretely expressed by the results shown in Figure 4.
We therefore change our approach from previous work, and no longer randomly generate the sneutrino mass m 2 ν c 3 at the unification scale. Instead, we specify the desired B-L breaking scale and use (31) 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 degree of fine-tuning is of O(10 4 − 10 6 ) for B-L scale of order 10 10 GeV and of O(10 2 − 10 3 ) for B-L scale of order 10 12 GeV. We note in passing that all the black points in Figure 5 are in the so-called "upside-down" hierarchy, with M BL < M SU SY .

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 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 (21), (22) and (23). The values for the parameters m, Y ν3 and µ will be chosen to be those given in (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 sof t 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 sof t . Then the relevant equations of motion are given by against the B-L scale, for the valid black points shown in Figure   5 from the scan of 50 million sets of initial conditions. The quantity where 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 post-inflation 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 (35) 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 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, tosc 1.096 × 10 7 marks the time at which the potential energy is well approximated by V = 1 2 mψ 2 and tMD 1.387 × 10 7 is the time at which our analytic solutions for ψ and H become valid. 8(a), we can numerically show that 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 (32), (33) and (34) then become 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 (39). 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 (40), we have Plug these expressions into (37) and (38), yields where H 1 is the first order approximation to H. It then follows that d dt 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 (40). 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 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 (47) and (49) into equation (39). This gives The solution is where A 1 is a constant. Putting (47), (49) and (53) into (37), we find that (54) 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 (56) 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 M D and indicated it by a red line in Figure 8(a). To summarize, when t > t M D 1.387 × 10 7 where c 2 = t osc − 2 3H(tosc) 9.67 × 10 6 and c 1 9.78 × 10 6 .
The numerical values of H and (ψ(t)/3) 2 at t * , t end , t osc and t M D 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 M D −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 (57) and for ψ(t) in (58), we compare them in Figure 10 (a) and (b) respectively against the exact numerical solutions for H and ψ in the region t > t M D . It is clear that our analytic solution is a very accurate approximation.

Post Inflationary Epoch: Decay of the ψ 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 (37), (38) and (39) 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 process commences; that is, when Γ d,i becomes nonzero. For convenience, we define the fraction of energy density of the i-th species as where, as follows from (59), 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 (13) to be with the associated quadratic soft mass squared given in (14) by ) .
The value of m was fixed as m = 6.49 × 10 −6 (66) to be consistent with the results of Planck2015 [51]. The relationship of H 0 u , ν L,3 andν c R,3 to φ 1 was presented in (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 (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 (68) can be simplified in the post-inflationary regime. Taylor expanding (68) around ψ = 0, we find As can be seen from Figure 9(b), in the matter dominated period t > t M D 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.
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 . 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.

Up-Type Standard Model Fermions
Note from (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 where y HF is the usual Yukawa parameter for coupling to the Higgs and, using (67) and (70), 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 Fig. 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 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. 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 (79) that Γ d will approach a constant. That is, which is the maximal value of Γ d . Using (66) and (76) we find that Γ max d,u = 1.801 × 10 −18 , Γ max d,c = 6.116 × 10 −13 , Γ max d,t = 1.585 × 10 −8 (82) 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 nonrelativistic 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 (59)- (62). 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 (59) and (60), as well as Γ d in (61). It follows that H can still be approximated by as in (57), where c 2 9.67×10 6 . Putting this back into (61) and taking Γ d = Γ max d , one can solve this equation for ψ. We find that where c 1 9.78 × 10 6 and Note that in the limit Γ max d → 0, this expression reproduces the result in (58). Putting expressions (81), (83), (84) and (85) into (62), 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 It follows from (63) and (87) 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 .
with t M D 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 ) 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 (92) 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 [57,30].
The mixing terms can arise from two sources; 1) the superpotential and 2) the"super-covariant 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 arẽ 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.
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 (193), 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 Note that (95) is a mass term for a Dirac mass with 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 where 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 intoC ± 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 .
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 intõ C ± 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 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 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 7, but note that to simplify our expressions we again make the assumptions given in (107) and (110). Additionally, we define Mass 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 whereÑ 2a andÑ 2b have the same mass presented in Table 2 and 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.

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 (73). Inserting this into the final expression in (116) and, as previously discussed, denoting δψ simply as ψ, we find that where the masses are given by (118) 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 respectively. For decays into two different W-bosons, one multiplies expression (119) by 2. It then follows from (117) that the decay rates of the inflaton to the four gauge bosons, and the associated equation of state parameters, are where the gauge boson masses are given in (118). 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 (110).

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 Since ψ 2 will decrease with time, ψ → ψ − d τ c R will be the first decay process to become non-zero, whereas the decay ψ →Ñ2Ñ2 will be turned on last. We have set MP = 1.
decay, one must examine the mass matrix M whose elements are where i, j run over all scalars other than the inflaton and To simplify our calculations, we assume that all scalars have identical soft masses given by 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.

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 (59)- (62) to account for all the relevant decay processes simultaneously. Therefore, in this section, we will numerically solve (59)-(62) 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.  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.

Initial conditions
To find the solutions for H, ψ and ρ i by numerically solving (59)-(62), 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 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 (96). This expression was first presented in (91) 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 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 (58) and (57), 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 Fig. 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 .
(135) 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 (59)-(62), 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 ψ 2dependence becomes complicate after t I . Hence, we will solve them by iteration. We accomplish this by using Eq. (58), 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 (59)- (62) 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 (59)- (62) 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 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 7, 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. (136) 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 Fig. 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 Fig. 13, while A n is the black curves in Fig. 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. 6 We plot the decay rates Γ d,i and the Hubble parameter H in Fig. 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 Fig. 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 Fig. 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 Fig. 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 Fig. 15. Ω N 2 N 2 includes the decay products for all three processes ψ →Ñ 2aÑ2a ,Ñ 2bÑ2b ,Ñ 2aÑ2b . We did not specify them individually because even their sum is very small. Note that i Ω i = 1 − Ω ψ by definition. Eventually, i Ω i → 1, since Ω ψ → 0. We can define the end of the reheating epoch as the time, t R , when Ω ψ → 0, which means that all of the energy of the inflaton has been converted to relativistic species of matter. It is clear from Fig. 15 that However, due to numerical imprecision, we may find that i Ω i ≈ 0.9999 at t R , which is more than sufficient for our purposes. Note that the values of the Ω i in Fig. 15(b) have been rounded to three decimal places. When added together, we find that they sum to 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 as C ± 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 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 5000 magenta points are randomly chosen from the curve of the output ψ(t), in the last (nth) round of iteration, during the time intervals tI ≤ t ≤ tII and t ≥ tII respectively. Their upper boundaries are identical to the input A(t) in the last round of iteration, that is, An, which is the black curves in (a) and (b). In (c), we plot the upper boundaries of the magenta points, that is, An+1, from tI 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 , tW W * = 6.93 × 10 9 , t tt * = 6.95 × 10 9 , tBB * = 7.06 × 10 9 and tNN * = 7.08 × 10 9 . In (d), we plot the solution for H in the last round of iteration, that is, Hnum, and also (57) as a reference, from tI to sometime after reheating. We always set tI = 8 × 10 8 and tII = 5 × 10 9 in our numerical calculations for technical convenience and everywhere set MP = 1.

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 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 and, hence, 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 [58]. 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 (146) 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 (146) 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, see (132), 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 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 [59]. 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 (147), 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 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 The rate ) is plotted as a function of time along with the Hubble parameter in Figure 18. It is clear that condition (146) 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 (146) 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 [60,61,58]. 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 scatteringalthough 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.

Appendix A
In this Appendix, we present the details of the renormalization group equations that allow us to specify the desired B-L breaking scale, and then run in energy-momentum up to the unification scale M U to determine the associated initial values for both m 2 ν c R,3 (t U ) 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.
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 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 α BL (t BL ) = 2 5 The beta functions for the various regimes are b BL,S = 6, b R,S = 7 and b BL = at a given scale. We find that for t SU SY < 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 (166) and (167), we find that However, we must now use the fact that S B−L (t U ) and S R (t U ) contain m 2 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 (169) and (170) are evaluated at t = t U . Doing this, and re-arranging terms to extract , we arrive at the expression

A.2. Cosmological Constraint
Recall from (14) and (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 (172) 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 (171) via the S-terms, we can see that equations (171) and (172) are intertwined and must be solved simultaneously. To do this, we can re-express these equations in the form where The equations in (173) are then soluble by inverting the matrix expression to give and inserting A, B and C in (174). In particular, for m 2 The solutions to m 2 Hu and m 2 ν c R,3 given by equation (177)  (M U ), a condition which defines the third family right-handed sneutrino. If this is not satisfied, we discard this set of initial data and then return to the first step. 7. The complete set of soft masses, the final value of M SU SY and the inputted value of the B-L breaking scale are then used to carry out the remaining physical checks on the scale of electroweak breaking, the Higgs mass and the sparticle mass bounds in the manner described in [30].
In the text of this paper, we will use this formalism to generate physically acceptable valid black points whose B-L scale is in the range 10 10 GeV ≤ M BL ≤ 10 12 GeV. This range comfortably accommodates our theory of reheating. Be that as it may, it is of interest to see whether the B-L scale can be substantially reduced to much lower values. With that in mind, in this Appendix we will implement the above procedure and 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 6 GeV and 10 14 GeV. Carrying out our checks, we find that this ultimately leads to 305 sets of initial data which satisfy all phenomenological constraints 1), 2) and 3) presented earlier. These physically valid black points are shown in Figure 20. The distribution of the B-L breaking scale for the black points in Figure 20 is given in Figure 21and verifies that our approach has indeed allowed us to extend the range of the B-L scale to lower values.

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 reexamine equation (171). Re-arranging this equation as an expression for M BL , we find 1 2 This can be schematically expressed as Equation (181) 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 against the M BL for those sets of initial data which survive all the phenomenological checks previously outlined. The quantity (183) is the Barbieri-Giudice (B-G) sensitivity [62] for the B-L breaking scale. The results of doing this for the 305 "black points" in Figure 20 are shown in Figure 22.  Figure   20 from the scan of 50 million sets of initial conditions. The quantity 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 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 (193) as where 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 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.

The Small µ Limit
When terms involving µ are dropped, the system to be diagonalized (presented in (197) ) simplifies tremendously. We are able to decouple the τ c R , ψ − d states as there is no longer any mixing between ψ − d and ψ + u . Indeed, examining the effective mass Lagrangian in equation (193), we see that which looks like a Dirac mass for the fermion with mass y τ ν L,3 . We can then simplify the system given in equations (196), (197) to Expressing the matrix X schematically as where x 1 = M 2 , x 2 = g 2 ψ / √ 6, it follows that The eigenvalues of (206) are then given by The matrices U and V become It follows that the associated mass eigenstates are with masses mC 1 , mC 2 . 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 where 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 where ψ T = B ,W R ,W 0 , ψ 0 d , ψ 0 u , ν L,3 , ν c R,3 and For simplicity, let us schematically re-express this to give x 8 x 9 0 0 0 0 0 x 10 0 0 0 x 5 x 8 x 10 0 0 0 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 the following table.
Without the VEVs, the last six terms in (220) are 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 since N iu = N iL = N iν = 1 √ 3 . 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 ψÑ † 2x ( p 1 , s 1 ) Figure 24: Feynman diagrams which contribute to the processes ψ →Ñ2aÑ2a and ψ →Ñ 2bÑ2b .
This allows the the inflaton decays ψ →Ñ 2aÑ2a , ψ →Ñ 2bÑ2b and ψ →Ñ 2aÑ2b . 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 with α = β = iO c . The above expressions simplify to give equations (113) and (114).