The Spectrum of Dark Radiation as a Probe of Reheating

After inflation the Universe presumably undergoes a phase of reheating which in effect starts the thermal big bang cosmology. However, so far we have very little direct experimental or observational evidence of this important phase of the Universe. In this letter, we argue that measuring the spectrum of freely propagating relativistic particles, i.e. dark radiation, produced during reheating may provide us with powerful information on the reheating phase. To demonstrate this possibility we consider a situation where the dark radiation is produced in the decays of heavy, non-relativistic particles. We show that the spectrum crucially depends on whether the heavy particle once dominated the Universe or not. Characteristic features caused by the dependence on the number of the relativistic degrees of freedom may even allow to infer the temperature when the decay of the heavy particle occurred.

After inflation the Universe presumably undergoes a phase of reheating which in effect starts the thermal big bang cosmology. However, so far we have very little direct experimental or observational evidence of this important phase of the Universe. In this letter, we argue that measuring the spectrum of freely propagating relativistic particles, i.e. dark radiation, produced during reheating may provide us with powerful information on the reheating phase. To demonstrate this possibility we consider a situation where the dark radiation is produced in the decays of heavy, non-relativistic particles. We show that the spectrum crucially depends on whether the heavy particle once dominated the Universe or not. Characteristic features caused by the dependence on the number of the relativistic degrees of freedom may even allow to infer the temperature when the decay of the heavy particle occurred.
Introduction.-Most models of modern particle cosmology predict a reheating phase. While there is no shortage of possible scenarios (cf., e.g. [1] for a review) remarkably little is known for sure about this crucial event in the history of the Universe. For example, the reheating temperature could have been as low as ∼ few MeV [2][3][4][5][6][7][8][9][10][11] or it could have been much higher ∼ 10 16 GeV (cf., e.g. [12]). It is therefore an interesting task to garner evidence for the existence of reheating and find ways to collect information on its details.
One reason that makes it hard to probe reheating is that during this phase the Standard Model (SM) particles thermalize and therefore most information carried by them is lost. One way to overcome this challenge could be very weakly interacting particles that are created during reheating and freely propagate until today. Suitable "messengers" may be gravitons/gravitational waves [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27], axion-like particles (ALPs) [28][29][30][31], right handed neutrinos [32] or some other very weakly interacting new particle. Preferably, such a messenger should remain relativistic until today, thereby constituting dark radiation. The reason for this preference is that for relativistic particles it is easier to relate the energy/momentum measured in a local experiment on Earth to the one imprinted at production. In particular the energy and direction are less affected by structure formation.
In a general setup, the messengers can be produced from inflaton or modulus decays. In a previous paper by the present authors [32], the transparency condition of the Universe for the messengers, including neutrinos and relativistic dark particles, were clarified. Satisfying these conditions, they can travel over the thermal history until today, and the momentum distributions carry the information of the mother particles and thus can be messengers of the reheating.
The possibility to detect ALPs originating from modulus decays was already discussed in the IAXO white paper [33] (we also note that the authors of Ref. [30] even already calculated the energy spectrum of the ALPs and gave analytic expressions for decays taking place in ei-ther pure matter or radiation domination 1 ). In [32] we discussed further possibilities such as right-handed neutrinos and more general messengers in dark matter, neutrino, and cosmic microwave background experiments.
Beyond the ability to detect the messenger particles the next task is to establish their origin from reheating and to obtain additional information. One feature of the messenger spectrum from reheating is the isotropic angular distribution which can be distinguished from particle spectra from galactic sources [30,32]. However, an isotropic distribution arises as long as the decays of heavy non-relativistic particles happens much before today [34,35]. Thus, one cannot say confidently that an isotropic spectrum must be from reheating.
In this letter, we therefore ask whether we can get clearer evidence of reheating by carefully studying the energy spectrum of freely propagating relativistic particles, which are produced in the decays of heavy nonrelativistic particles. We show that this differential flux contains information on the equation of state of the dominant energy of the Universe during the time when the decays occur (cf. the next section). We find the shape of this flux and point out that by precisely measuring the spectrum, we can get strong evidence for reheating. Moreover, changes in the number of available degrees of freedom with the temperature also leave imprints in the messenger spectrum. Resolving these may give direct information on the reheating temperature (see further below and also Appendix A). Combining this with the peak position of the spectrum [32] may even allow to determine the mass of the heavy particle reheating the Universe.
Our analysis is independent of the specific type of messenger. The results should therefore equally apply in the case of the aforementioned or other messengers, that are sufficiently weakly coupled relativistic particles. 2 The main assumption we make is that the messenger is produced in a two-body decay from the precursor particle. A similar analysis is also possible if the decay of the precursor is more complicated (an example is discussed in Appendix B 4) but we expect that the identification of features will be less clean.
In practice, of course, the detection will depend on the type of messenger realized. While the required measurements certainly are quite challenging we can nevertheless conclude (see also Appendix C) that, given the existence of suitable messengers, future measurements at observatories such as IAXO [33,[56][57][58], IceCube [59][60][61][62][63] or DARWIN [64] may shed at least a little light on reheating.
Messenger flux in the expanding Universe.-In this section we start by considering the simple example, where a heavy non-relativistic particle is responsible for the reheating, i.e. it was the dominant form of energy before it decays. The primary decay into SM particles causes the reheating, but a rare decay will serve as the source of our messenger. This situation will turn out to be clearly distinguishable from the case where the initial heavy particle only contributes sub-dominantly to the energy density of the Universe and therefore cannot be the main source of reheating.
Setup.-Let us consider a non-relativistic real scalar field, φ, with mass, m φ , decaying at the cosmic time, t ∼ t decay . We have in mind the case that φ is a scalar modulus or inflaton, but in general, we can also consider the situation where φ is a fermion. This does not change the main conclusions. For us a relevant feature is that the decay is two-body. This is well motivated in many models, e.g. modulus/inflaton decay into axions/ALPs [28][29][30][31] or decays into right-handed neutrinos [32,[65][66][67][68][69]. For our purposes this has the advantage that the produced messenger particles have a definite energy at the time of the decay instead of being distributed over phase space, thereby their energy today contains direct information on the time of decay.
Let us assume that φ has the decay channel where χ is a particle that freely travels to Earth at a velocity close to the speed of light. The decay rate is 2 It could also apply to gravitational waves originating from particle decays [22] but their frequency may be too high to be detected in near future.
where Γ tot is the total decay width of φ and Br φ→χχ ≤ 1 is the branching fraction of φ → χχ. For the decay time we have, We usually assume Br φ→χχ 1 such that the primary decays are to SM particles 3 .
The expansion rate of the Universe is given as where a is the scale factor and ρ φ (ρ r ) is the energy density of φ (radiation), and M pl ≈ 2.4 × 10 18 GeV is the reduced Planck mass. In line with our assumption Br φ→χχ 1 we take ρ r to be dominated by SM particles and neglect the component of χ. We then have, Here, g (and g s appearing soon) is the number of relativistic degrees of freedom (for entropy) in the SM for which we use the values and behavior taken from Ref. [72].
Numerical results for the flux of messenger particles.-We can obtain the time dependence of the energy density of φ as well as the radiation froṁ Here, s r = 2π 2 g s T 3 /45 is the entropy density and The form of c[t] is obtained froṁ ρ r Γ tot ρ φ and a prompt thermalization in a much shorter period than the expansion of Universe. We have neglected the dark radiation contribution in s r . If Γ tot ρ φ → 0, s r a 3 conserves, as expected. From these equations we can then also obtain the time dependence of the Hubble parameter H.
Using the time dependence of ρ φ , ρ r and H, we can now obtain the differential flux of χ. We start with the momentum distribution of χ. This can be obtained from the Boltzmann equation, where f k is the distribution function of χ of momentum k (we have assumed rotational invariance). This formula is justified when the occupation number of χ is not too high (c.f. Refs. [55,73]), and m χ is the mass of χ. By assuming f k (t t decay ) = 0 and neglecting the mass of χ here and hereafter, we obtain the solution where t satisfies a(t)k = a(t )m φ /2, and θ is the step function. Here, the dependence on t , which in turn depends on a, makes explicit the connection between the spectral shape and the expansion history. The differential flux today (at t = t 0 ) is then given by where E = k 2 + m 2 χ is the energy of χ and we have neglected m χ on the r.h.s. (we will also do so in the following). Importantly one finds that the distribution function depends on the Hubble parameter H(t ) at the time of the decay that contributes to the flux at k. This will make it possible to probe the expansion history in the early Universe by precisely measuring the differential flux. Now, we can calculate the flux by solving Eqs. (6) and (7). Let us define CASE A and B from two initial conditions depending on whether φ once dominated the Universe or not: Here, we have fixed the branching fraction for CASE A to be Br. The shape of the spectrum does not depend on the value of Br as long as it is much smaller than 1. The overall signal strength is simply proportional to it. In both cases we fix Γ tot = g (T φ )π 2 T 4 φ /(90M 2 pl ), with T φ = 200 MeV being the decay temperature, i.e. φ decays at the same time scale in both cases. This is also roughly the reheating temperature if the decays reheat the Universe. We also fix m φ = 100 TeV. The branching ratio in CASE B is taken to fit CASE A in the high energy region. The numerical result for the spectrum Fig.1. 4 Here, we divide the differential flux by Br. In this way our result applies to a wide range of Br unless it is so high that χ becomes the dominant component of the Universe.
We find that at energies below the peak the spectrum clearly depends on whether φ once dominated the Universe or not. In particular, the slopes of the fluxes (in the log-log plane) differ by O(1) before the peak. 4 In our logarithmic plots we use dΩd log 10 E because this facilitates getting an impression of the total flux in an energy interval by simply multiplying the plotted value with the logarithmic width of the interval.

C A S E A ( R e h e a t i n g f l u x ) C A S E B
CA SE C Figure. 1. The differential flux of the messenger particle, d 2 Φ/d log 10 EdΩ. CASE A (φ once dominated the Universe) and CASE B (φ never dominates the Universe and decay in the radiation dominant epoch) are shown in red and black lines, respectively. We also show the flux for CASE C where a subdominant φ decays in the matter dominant era as the blue dashed line.
As we will see from the analytical estimates below, the slope of the flux below the peak depends mainly on the equation of state of the dominant energy form. Nevertheless, the flux in CASE A can be also distinguished from the case of a subdominant φ decaying during a matter dominant epoch. This is shown as the blue dashed line (CASE C) in Fig. 1 where the decay rate of φ satisfies Γ tot = H| z=100 i.e. in the matter dominant era. We use

and the forms of H and t[z]
given in Ref. [32]. (This form can be obtained by solving Eq. (6) and then using Eqs. (9) and (10) while taking H in the matter dominated era and neglecting the energy density of φ.) We take C mat ≈ 6.1 × 10 −42 GeV 4 , m φ = 5 keV. to match the IR flux to the one of CASE A. We can see that with E around or above the peak energy, we have more than O(10)% flux differences. More drastic differences may be observed if φ decays in the dark energy dominant epoch. In particular, if φ is decaying today (and has been nonrelativistic for a sufficiently long time), φ tends to gather around the galactic center due to the gravitational interaction. In contrast to CASE A this flux has an angular dependence that (amongst other things) superimposes a narrow peak at m φ /2.
Analytical approach to the flux.-To better understand the different behaviors and the ways to distinguish the different scenarios let us consider some analytical estimates. Note that f E for a given energy E is proportional to H −1 at the scale factor a = 2E/m φ . H depends on the dominant energy of the Universe at a, which means that the differential fluxes, or spectral intensities, at different E scan the dominant energy of the Universe at different a.
Let us consider the epoch before the typical decay time t < t decay = 1/Γ tot , i.e. a a decay with a decay being the scale factor at Γ tot = H. Since f E ∝ a −3+3/2(1+w) in this epoch we obtain (in agreement with [30] where w is the equation of state of the Universe, e.g. w = 0 (1/3) for matter (radiation) dominated Universe. Therefore, the slopes of the fluxes (in the log-log plane) differ by O(1) for CASE A and B when E < E decay = m φ /(2a decay ). χ with E > E decay , are produced from the few remaining φ after the typical decay time. Accordingly, the flux gets an additional exponential suppression.
If φ reheats the Universe, i.e. it dominates the Universe before the decay and transfers most of its energy into radiation around and after the decay, the differential flux of χ has two typical regimes where κ is an order 1 numerical coefficient. Such a reheating flux can also be distinguished from the decay of a subdominant species during a matter dominated phase of the Universe (CASE C). The reason is that if the decaying particle is responsible for reheating, decays after the typical decay time occur during a radiation dominated epoch. For a subdominant species such a changeover at the decay time usually does not happen. 5 This is then reflected in the spectrum at energies higher than the peak energy, where we have (again in accord with [30] Note that the relative shapes of the differential flux do not depend on parameters of the model such as Br φ→χχ or m φ . Thus the reheating flux shape is a robust prediction of reheating caused by non-relativistic particle decays to relativistic messengers. If we can identify this flux, it should be strong evidence of reheating. Measurement of reheating parameters.-Let us now outline a strategy for an ultimate possibility to measure reheating parameters in the future. A slight change of the equation of state of the Universe will also be induced by the decoupling of the relativistic components in the radiation dominant era. Since soon after the reheating the Universe is radiation dominated, this effect affects the reheating flux shape. We depict a comparison of reheating fluxes in Fig. 2 with and without taking account of this decoupling effects in red-solid and gray-dotted lines, respectively, with T φ = 400 MeV. The CASE A (T φ = 200 MeV) flux is also shown in the blue-dashed line. The decoupling affects the flux via the changes of g , g s , which is especially significant around the QCD phase transition. (See more details in Appendix A) Again these flux shapes do not depend on Br and m φ . Neglecting the decoupling effect (gray-dashed line) the shape even does not depend on T φ . 5 Such a changeover at the same time as the decay would be an unlikely coincidence. Conversely by carefully measuring the spectral shape, the reheating temperature ≈ T φ can be obtained in principle. For example, in Fig. 2 we can infer that the flux indicated by the red solid line originates from around QCD phase transition era. In addition, we can then obtain m φ via a determination of m φ /T φ from the peak position [32]. Then, T φ can be translated into Γ tot as well as ρ r and thus ρ φ from energy conservation. From the intensity of messenger flux we can also derive Br φ→χχ .
Conclusions and discussion.-Reheating is a central part of our current modelling of the early Universe. However, it has not yet been confirmed by direct experiment/observation.
In this letter, we have studied the spectra of freely propagating relativistic particles produced by the decays of heavy non-relativistic particles and shown that the spectra depend significantly on whether the heavy particle once dominated the Universe. We have demonstrated that the energy spectrum of relativistic messengers from the decay of a non-relativistic particle exhibits clear features that can tell us whether the decaying particle was responsible for reheating. These imprints arise via the equation of state. If the spectrum can be measured with sufficient precision we may even be able to tell the reheating temperature and the mass of the decaying particle. One may also wonder what happens if the mother particle is relativistic, featuring its own non-trivial spectrum. In this case (cf. Appendix B where we consider typical spectra arising in preheating scenarios) one cannot easily tell whether the decaying particle dominated the energy density or not. However, it is nevertheless usually distinguishable from the case of a non-relativistic mother particle. Furthermore, even if the messengers arise from a cascade of two subsequent two body decays, some small traces of the phase during which they originated may still be visible (cf. Appendix B 4).
Our discussion here can apply to various reheating sce-narios with mother particles coupled to light-weakly coupled messenger particles, such as gravitons, ALPs, dark matter, neutrinos, etc.. Let us conclude with some comments on the experimental opportunities. The discrimination between CASE A and B is relatively straightforward, since the flux at energies below the peak are quite different in a large range. In Appendix C, we argue that, even if we can just measure the flux in two bins with a relative width of 25% corresponding to a quite moderate energy resolution, the discrimination is possible at the 2σ level with O(1000) events, possibly even significantly less if the data is used more efficiently than in our simplistic estimate. This gives us an optimistic expectation for suitable experiments such as IceCube, IAXO and DARWIN [33,[56][57][58][59][60][61][62]64] (see also [74]).
The discrimination between a reheating flux and subdominant φ decays in a matter dominated epoch as well as a measurement of the reheating temperature is more difficult due to the exponential suppression of the flux at high energy. A good energy resolution (expected for at least some of the above mentioned experiments) will be critical to do this. As an additional check of the possible reheating origin, the angular distribution of the flux can also play a useful role [30,32,34,35]. The last question is whether enough events can be observed. In fact, Ice-Cube [63,75,76], ANITA [77,78] and XENON1T [79] experiments have already observed anomalous events that may be from BSM physics. 6 As an example we show in Fig. 6 (in Appendix B 4) fluxes for the different cases discussed in this letter fitting the events [63] of IceCube. Aside from inviting intriguing speculation this also suggests that a sufficiently good measurement to distinguish the fluxes may be feasible in the not too distant future. If these hints persist and are not explained by other effects, we may be in the fortuitous situation that measuring their energy dependence may allow us to get a glimpse of the beginning of the thermal history. The shape of the flux also depends on the changes in the number of relativistic degrees of freedom. This becomes noticeable for φ decays that happen in or close to the radiation dominated era. This holds because when from Eq. (9). The effect can be clearly seen as the slight kinks near the peak of the spectrum of CASE B in Fig. 1. The changing of the slope just before the peak reflects the changes of g , g s due to the QCD phase transition.
This decoupling effect is also important for the reheating flux as discussed in the main part since soon after reheating radiation is dominant. The time evolution of T and ρ φ in CASE A is shown in Fig.3. When t × Γ tot ∼ 1 reheating ends which also roughly corresponds to the time when the peak of the reheating flux is created. We see that at T ∼ 200 MeV the QCD phase transition slightly slows down the decrease of T . This is different from the flux with constant g and g s . A comparison of the fluxes is shown in Fig. 2. Here the flux divided by a parameter C is shown. The CASE A flux is given in blue dashed line (C = Br) and the case with constant g and g s is indicated by the gray dotted line. For the gray dotted line, T φ = 400 MeV, C = 1.5Br and m φ = 240 TeV as well as assuming constant g (T ) = g (T φ ) ≈ 66.9 and g s (T ) = g s (T φ ) ≈ 66.3. We have also checked that the shape does not change when changing T φ , Br, m φ as long as we take constant g and g s , as can be expected. However, the peak energy and flux intensity may change. The difference can be more significant if T φ is slightly larger so that the QCD phase transition happens mostly in the radiation dominant epoch but with not too suppressed ρ φ . The red solid line represents C = 1.21Br with T φ = 400 MeV, m φ = 221 TeV. As we can see the fluxes can differ noticeably. By carefully measuring the energy scale and the size of the depression, in principle, we can identify T φ .
This discussion applies not only to the reheating flux but also in the more general case where we can measure T φ given a flux from φ decays in the radiation dominated epoch. By carefully measuring the flux at two energies, we can in principle measure g −1/2 g −1

S
at two different temperatures around T φ . Then T φ can be obtained under the assumption that the SM accounts for the D.O.F. Appendix B: Messenger flux from the decay of relativistic particles

Setup -Relativistic particles from preheating
If the reheating is caused by the decay of an oscillating scalar field with a non-parabolic potential, φ may be "preheated" [36][37][38][39][40][41][42] and have a non-trivial spectrum before the decay. Let us consider the φ potential and assume a large enough initial field (or tiny enough mass term). When the amplitude is large, parametric resonance [36][37][38][39][40][41][42] 7 becomes important. Soon afterwards the system enters into a turbulence regime in which a scalar field spectrum follows a power-law [100][101][102][103]. For a sufficiently long period, almost no homogeneous mode of φ remains. The relativistic φ particle spectrum undergoes a self-similar evolution until the decays. A typical expectation for the occupation number is [100][101][102][103] More generally we can consider a power-law spectrum of the form, We will now focus on the phase when such a power-law has been established, the relativistic φ dominates the energy density and its decay into SM particles reheats the Universe. Via couplings to the SM particles and χ, φ can decay both into SM particles and χs. Again we assume a two body decay to χ, φ → χχ. (B3) Moreover, we assume that these couplings are so small that the decay processes can be treated in perturbation theory. The decay rate of the k mode of φ is where the pre-factor is the Lorentz factor, and again Γ tot is the total decay width in the rest frame. We can get the decay rate to a χ pair as Assuming k m φ , we have Γ φ→χχ [k], Γ tot [k] ∝ k −1 .

Numerical result
The equations for the time dependence of ρ φ and ρ r (neglecting the small χ density), arė . H is given in (4). The relevant modes of φ are assumed to behave as radiation. Indeed, hereafter we mostly neglect m φ , m χ .
The χ differential flux can be obtained by assuming a subdominant decay to φ → χχ with Br φ→χχ 1. In addition to Eqs. (B6) and (B7), we can solvė and P [k, k ] = 2(k/k ) 2 θ(k − k) represents the phase space distribution of χ from a relativistic φ decay with m φ m χ . We get the differential flux, These equations can be applied to the χ flux for any initial relativistic spectrum of φ decaying into SM radiation and χ. Now let us come back to the preheating scenario discussed above. The initial conditions for the reheating from the self-resonant φ decays are 1 so that before the decay φ particles dominate the Universe. The exponential in ρ φ,k (t ini ) represents the decay of the power-law set by hand, and k min,max is set for the convenience of calculation. We also set by which we define the reheating temperature. For convenience we use In Fig. 4, we show the resulting cosmological temperature (red solid line) and the energy density of ρ φ (black solid line). The corresponding differential flux, φ→χχ , is plotted in Fig. 5. For comparison, we also show the initial flux of φ. One can see that the χ flux mimics the original φ flux for n = 5/2 but not for n = 3/2.

Analytical understanding and general features of the flux from relativistic particle decays
The behavior of the χ flux for different n can be understood from the equatioṅ which is derived from Eq. (B8), by definingρ χ,k ≡ a 4 ρ φ,k andρ φ,k (t) ≡ a 4 ρ φ,k withk ≡ a × k being the comoving momentum. In the second equality we have usedρ φ,k (t) = exp − t tini dt Γ tot (k ) ρ φ,k (t ini ). Integrating over time we havê where we have assumed t much larger than the typical decay time. Due to the combination of the step-function in P [k,k ] and the power ofk in the integrand, the slope of ρ χ,k at the IR end can be approximated aŝ Here, we have used that thek integral is dominated by values ofk aroundk for n > 2. In contrast for n ≤ 2 the integral is dominated by the UV cutoff (for n = 2 only logarithmically) and thus approximately independent of k. Therefore, in this caseρ χ,k is proportional tok 2 from P [k,k ]. This agrees well with the numerical results. Translating this into the flux, we obtain a maximal slope, d 2 Φ/dΩdE ∝ ρ χ,k /k 2 ∝ E 0 . As a result, the reheating flux (∝ E 1/2 ) cannot be mimicked by the decays of relativistic particles with a simple power-law.
From Eq. (B15) we can also see that the form of ρ χ,k does not depend on H in the early Universe as long as all φ particles decay. This is different from the case of nonrelativistic φ decays. We therefore can neither distinguish if φ dominated the Universe nor infer the behavior of H by the messenger flux from relativistic φ decays.

An example with a cascade decay
Another plausible scenario for the origin of a relativistic mother particle of the messenger is a cascade decay. Let us consider a situation where a non-relativistic particle undergoes a two-body decay which is followed by another two-body decay to our messenger particle. The spectrum of the relativistic intermediary particles is that of the messengers considered in the main text. From this we find that a decay in the matter (radiation) dominated phase corresponds to n = 3/2 (n = 1) for energies far below the peak. From Eq. (B16) we can now see that the resulting spectrum for the final messengers is then independent of the original slope n. The second decay, unfortunately obscures the information on the slope.
However, there are still some subtle traces of the expansion history/slope left. The reason is that the k 2 behavior of Eq. (B16) is obtained in a region where the integral in Eq. (B15) is completely dominated by the UV. However, when k approaches the cutoff, i.e. the end of the power-law behavior, the lower boundary of the integral becomes relevant and differences in the shape of the original spectrum start to matter.
To see this let us consider a concrete example: φ → N N → ννππ with N (ν, π) being the right-handed neutrino (left-handed neutrino, pion). Here, for simplicity, we assume that the two body decay of N is dominant, which happens in a certain parameter range [104]. Moreover, we neglect the masses of the decay products as well as neutrino flavor. As discussed in [32], this process can carry the information from reheating to Earth without interacting with the plasma.
We consider, again, three cases for the flux of N . In CASE A' φ reheats the Universe. In CASE B' and C', a subdominant φ decays in the radiation dominated and matter dominated epoch, respectively. More precisely, we take the parameters for CASE A' corresponding to the CASE A set-up in the main text as In either case the other parameters are left unchanged. Then we obtain the ν flux by using (B15), i.e. d 2 Φν We add a 1/2 because on average half of the energy of N is transferred into ν. The resulting spectra are shown in Fig. 6. (Note that we display E 2 ν d 2 Φ ν /dΩdE ν to match the presentation chosen by IceCube.) The additional red solid line represents the reheating flux (CASE A) assuming a direct decay φ → νν with T φ = 10 MeV, m φ ≈ 9.5 × 10 15 GeV and Br φ→νν = 0.5 × 10 −7 . 8 For an illustration of the experimental situation we also display current experimental data [63] from the IceCube experiment. For now all curves are consistent with the measurement (we have chosen the parameters such that they roughly fit the present data around the bump). A data point around 1 PeV, i.e. just out of the range of the figure, is slightly above the curves. However, it may be easily explained by considering certain SM cosmic-rays e.g. [63]. On the contrary the displayed bump may be difficult to explain within the SM since it is above the Waxman-Bahcall bound [105]. But, envisioning future improvements, it seems conceivable that the curves can be distinguished.  Figure. 6. Neutrino spectra from a cascade decay φ → N N → ννππ. Here CASE A', B', and C' are shown in black solid, purple dotted, and gray dashed lines, respectively, representing the case that φ decays to reheat the Universe (T φ = 10 MeV), φ decays in the radiation dominant era (T φ = 10 MeV), and decays in the matter dominant era (z = 100). For comparison, we also show a reheating flux assuming φ → νν. 8 IceCube 7.5 years' data (blue points) is taken from Ref. [63].
We remark that the cascade decay spectrum of ν is a four body decay spectrum with a special phase space distribution. Therefore, even in the case of cascade decays or perhaps even multi-body decays, reheating may be probed in experiments by precisely measuring the spectra.
Appendix C: Statistics for measuring the reheating flux.
Let us roughly estimate what level of signal size we need in order to identify the reheating flux. For simplicity of discussion, let us suppose that the detector features only two energy bins close to the peak but in the powerlaw region of the flux. We can now try to discriminate the different power-laws of the flux that distinguish CASE A and B.
We expect the number of events in a small energy bin to behave as (assuming that the sensitivity of the detector is energy independent over the relevant region), where we have also assumed a reasonably small bin size ∆E and that the flux scales as a power γ of the energy. Considering two energy bins 1 and 2 covering the energy interval (E 0 + 2∆E, E 0 + ∆E) and (E 0 + ∆E, E 0 ) we can now ask how many events we need in order to distinguish two different powers γ such as 1/2 and 1, corresponding to CASE A and B, respectively.
Treating the counting errors ∼ √ N 1 ∼ √ N 2 in the two bins as statistically independent (adding the uncertainties in quadrature) we find, For an nσ distinction between two values of γ differing by δγ we therefore need, events.
Using ∆E/E 0 = 0.25 and δγ = |1/2 − 1| = 0.5 (CASE A vs. CASE B) this yields roughly ∼ 1000 events for a 2σ detection. Dropping the small (∆E/E 0 ) approximation this slightly increases to about ∼ 1400. That said, the procedure employed here is far from optimal. Choosing more and/or better separated bins this number can probably be decreased by an order of magnitude.
So far we have focused on the distinction between CASE A and CASE B. Distinguishing CASE A from CASE C will be significantly more challenging because the flux drops rapidly in the relevant region. In this case, we will probably need to fit the flux by using a sizeable number of bins. Here, a very good measurement of the energy of the events, allowing for small bins as well high statistics will be required.
To put the above example into context let us just remark that in various future experiments, e.g. IAXO [33,[56][57][58], IceCube [59][60][61][62] or DARWIN [64], the energy resolution for a new physics flux can be as small as d log 10 (E/ GeV) ∼ O(0.1) or even better, for at least some energy range. In addition increased sensitivity, potentially by orders of magnitude, may then allow to collect the required statistics.