Analytic solutions for neutrino-light curves of core-collapse supernovae

Neutrino is a guaranteed signal from supernova explosions in the Milky Way and is the most valuable messenger that can provide us with information about the deepest part of supernovae. In particular, neutrinos will provide us with physical quantities, such as the radius and mass of protoneutron stars (PNS), which are the central engine of supernovae. It requires a theoretical model that connects observables such as neutrino luminosity and average energy with physical quantities. Here we show analytic solutions for the neutrino-light curve derived from the neutrino radiation transport equation by employing the diffusion approximation and the analytic density solution of the hydrostatic equation for the PNS. The neutrino luminosity and the average energy as functions of time are explicitly presented, with dependence on PNS mass, radius, the total energy of neutrinos, surface density, and opacity. The analytic solutions provide good representations of the numerical models from a few seconds after the explosion and let our rough estimate of these physical quantities to be made from observational data.


Introduction
Neutrino is a guaranteed signal from a core-collapse supernova (SN) occurring in the nearby universe. We will firmly observe many neutrinos when an SN appears inside the Milky Way with currently running facilities and future ones. It is a strong point of neutrino observations comparing with other observable signals: For instance, even inside the Milky Way, the optical emission would be absorbed by the dense dust if an SN happens near the Milky Way center, and the gravitational waves, whose amplitudes are highly uncertain and model dependent, are difficult to predict the detectability at the current moment. For neutrinos, multiple running detectors, such as Super-Kamiokande [1], IceCube [2], NOνa [3], and also future detectors expected within the next decade like Hyper-Kamiokande [4], JUNO [5], and DUNE [6] will be able to detect many neutrinos from the next nearby SNe. Preparations for observation are well underway.
On the other hand, not much progress has been made in theory side. So far, only a few numerical studies for long-term protoneutron star (PNS) cooling associated with supernova explosions have been carried out using advanced technologies [7][8][9][10][11][12] (see Refs. [13,14] for SN1987A and Ref. [15] and references therein for short-term simulations). These studies just picked up a few progenitor models and presented a limited number of long-term evolution. To analyze the real data when an SN takes place, we need a systematic neutrino event rate evolution, i.e., neutrino-light curves. For this purpose, the detailed numerical simulations are not good enough for a parameter study, so that we need simplified analytic formulae of neutrino-light curves. The same analogy applies in modeling the optical light curves of SNe, where Arnett model [16], in which bolometric light curves are calculated by simply assuming photon diffusion in the expanding media, is widely used to extract physical parameters from observational data. This model allows one to extract the mass of radioactive 56 Ni, the mass of the ejecta, and the velocity of the ejecta, with the bolometric light curves. In this paper, we give such analytic solutions for neutrino-light curves, which give simple and useful rules.
Although the neutrino transport is a complex phenomenon dependent on the detailed structure of the ambient matter and the neutrino spectrum, we derive simple analytic formulae. We approximate the PNS's density structure by the Lane-Emden solution (Section 2) and also approximate the neutrino transfer with the diffusion process (Section 3). Combining two components (early-time and late-time solutions) gives a good agreement with detailed numerical solutions. Based on the analytic formulae, we give a simple expression of the neutrino detection rate and positron energies, which are useful for the data analysis (Section 4). The caveats are discussed in Section 5 before summarizing in Section 6.

Density and temperature structure
We start from an analytic solution of Lane-Emden equation for n = 1 as follows [17], where ρ c is the central density, ξ is dimensionless radius. With PNS mass M PNS and radius R PNS , we can write r = αξ with α = R PNS /π and ρ c = M PNS /4π 2 α 3 . The surface of PNS is given at ξ = π.
Assuming the constant entropy in whole PNS gives the temperature structure. The entropy per nucleon is [18] where k B is the Boltzmann's constant, T is temperature, n N is nucleon's number density, n i is number density of protons (i = p) and neutrons (i = n), and µ i = ( 2 /2m)(3π 2 n i ) 2/3 is chemical potential of nucleons with m being the mass of nucleons. Here, we omit finite temperature correction, since we are now interested in relatively cold PNS. Since n N = ρ/m, we get with Y e being the electron fraction, i.e. n p = n N Y e and n n = n N (1 − Y e ). With Y e = 0.1, f (0.1) = 0.699 and Y e = 0.5, f (0.5) = 0.630. Hereafter, we 2/12 fix f (Y e ) = 0.630 for simplicity because f (Y e ) changes only slightly for the typical value of Y e . With thermal energy per nucleon [18] the total thermal energy is given by Note that the entropy outside the neutrinosphere might be different, but the position of neutrinosphere is almost at the surface of PNS for timescale we are interested in (see Eq. 40). Thus we assume that the entropy is constant in the whole PNS for simplicity. For comparison, the total gravitational binding energy of this density structure is given as which is larger than the case of constant density that gives 3.1 × 10 53 erg.

Neutrino transfer equation and solutions
The Boltzmann equation in the spherical symmetry within O(v/c) is given by [19] df cdt + µ ∂f ∂r + µ d ln ρ cdt + 3v cr where f = f (t, r, µ, E) is the distribution function of neutrinos, µ is the cosine of the angle between radial direction and neutrino propagation, E is the neutrino energy, v is the fluid velocity with respect to the laboratory frame, c is the speed of light, j is the emissivity, χ is the absorptivity, and R is the isoenergetic scattering kernel. We denote the Lagrangian time derivative in the comoving frame by d/dt. Hereafter, we assume v = 0 because we treat a static PNS as the background matter. For deep inside the PNS, the distribution function can be approximated by f (t, r, µ, E) ≈ f (0) (t, r, E) + µf (1) (t, r, E). With it, the specific energy density and the specific flux are given by where Ω is a solid angle and h is the Planck's constant. By introducing them to Eq. (9) and taking the zeroth and the first angular moments, we get the following two equations: and assume f (0) = f ′(0) because we employ only the elastic scattering for R. By omitting the term ∂F E /∂t and introducing Eq. (13) into Eq. (12), we get Since for the local equilibrium state the RHS of Eq. (14) vanishes, the equilibrium distribution function is given as . We rewrite the Eq. (14) as where we also rewrite opacities with the following variables, κ t := j + χ + φ and κ a := j + χ. It gives the time evolution of neutrino specific energy density.
Integrating over the energy gives where with a = 8π 5 k 4 B 15(hc) 3 = 7.56 × 10 −15 erg cm −3 K −4 being the radiation constant. Here, we employ the Fermi-Dirac function without the chemical potential for the equilibrium spectrum of neutrinos as ε eq E = 4πE 3 (hc) 3 1 1 + e E/kBT . The opacities are expressed by the following mean as, The integrated flux is given by In the following, we assume that 1/ κ t = 1/ κ t eq and κ a = κ a eq , since inside the PNS the spectrum is largely determined by the equilibrium state so that ε E ≈ ε eq E . Then, taking κ t =κ t (E/m e c 2 ) 2 , κ a =κ a (E/m e c 2 ) 2 , with m e being the electron mass, we get Opacities are given by Eqs. (21) and (22) of Ref. [20] as Here, we fix Y = 0.5 of Eq. (25) for simplicity (see [20] for detail). β is a boosting factor of the scattering due to the existence of heavy nuclei, which significantly amplify the scattering cross section by the coherent scattering [21], in the PNS crust. β ≈ 3 for the free nucleons, in which the absorption is also taken into account.
With thermal neutrino spectrum, we rewrite Eq. (16) as, From Eq. (27), we can estimate the following two timescales; the absorption and emission timescale, τ a , and the diffusion timescale, τ dif , as By applying the approximation ε = ε eq , we have the following two equations, as The boundary condition is given by F eq = const. near the surface of a PNS, where the plane-parallel approximation can be applied. On replacing r by the mean optical depth τ , Eq. (31) becomes Integrating Eq. (33) with respect to τ gives ε eq (τ ) = ε eq (τ = 0) 1 + 3 2 τ .
Here, we impose F eq (τ = 0) = c 2 ε eq (τ = 0) by assuming that the neutrinos at the surface is isotropic along the outward directions and vanishing in the inward directions. As ε eq = 7 16 aT 4 , we have where T 0 is the surface temperature, i.e. T 0 := T (τ = 0). It is related to the brightness temperature, which is given by F = 7 16 σT 4 br with σ = ac/4 being the Stefan-Boltzmann constant, by According to Eqs. (35) and (36), the brightness temperature becomes the same as the matter temperature at τ = 2/3, which is the definition of the neutrinosphere in the following. The radius of the neutrinosphere, R ν is given by the following equation, as where ξ ν := R ν /α andκ t =κ t /ρ. We introduce a new parameter g which accounts for the different structure of the PNS surface from the Lane-Emden solution with n = 1. 1 The integral can be approximated as π ξν sin ξ ξ The error of this expression compared to the numerical integration is at most 17.4% at ξ ν ≈ 1.2. Combining them gives It tells that the neutrinosphere locates near the PNS surface. 1 The g factor is typically smaller than unity since the polytropic index near the PNS surface is n ≈ 3 (relativistic electrons are dominant), which leads to the steeper decline of the density than n = 1. It is shown in Table 1, in which g = 0.04 -0.1 is used to fit the numerical solutions.
Next, we derive the time evolution of the neutrino luminosity. Since the PNS thermal energy is decreased by neutrino emission, where we take into account all six types of neutrinos. By combining Eqs. (5) and (44) where t is time and t 0 is the time origin, which gives the initial condition of the entropy. The corresponding neutrino luminosity is given as The average energy of neutrinos is given by In the model, there are five parameters, i.e. the PNS mass M PNS , the PNS radius R PNS , the density correction factor g, the opacity boosting factor by coherent scattering β, and 7/12 the total energy emitted by all flavors of neutrinos E tot . Note that the boosting factor β is time-dependent because the heavy nuclei in the crust are absent for the early phase and appear later once the temperature decreases below the Coulomb energy of the lattice structure [11]. Therefore, we propose a two-component model to reproduce numerical models of neutrino-light curves. The first component represents the early time without coherent scattering (β = 3) and the second component represents the late time with the opacity boost by the coherent scattering (β ≫ 1). The neutrino luminosity is given by the total luminosity of two components, L 1 + L 2 , and the average energy is estimated by the harmonic mean,     Table 1 Model parameters reproducing numerical solutions presented in Ref. [12].
converts to additional neutrino emission. In the very late phase ( 30 s), the approximations (e.g., the thermal spectrum, a constant g factor) may break down. Table 1 summarizes the model parameters reproducing numerical solutions presented in Ref. [12]. These parameter sets can be used for the mock data production to study the detector responses for the next nearby supernovae. They are also useful to complete the neutrino-light curve with the detailed (multi-dimensional) hydrodynamics simulation, which are typically calculated up to ∼ 1 s, by connecting the numerical data with the analytic formula for the PNS cooling phase. Here, we compare the PNS mass and radius obtained from the analytic solutions with those of the numerical models. The PNS mass is consistent with the numerical model (147S has 1.47 M ⊙ , M1L and M1H have 1.29 M ⊙ , and M2L and M2H have 2.35 M ⊙ for the baryonic mass) due to degeneracy with other parameters, but the PNS radius differs from the numerical model (147S has 14.3 km, M1L and M1H have 14.4km, and M2L and M2H have 13.4 km). It may be because the analytic solutions apply the Lane-Emden solution to the entire region, and is based on Newtonian gravity. The resolution of these deviations is beyond the scope of this paper, but we will improve them in future work.
Although we demonstrated fitting by the two components works quite well, we limit our discussion within the single component model fitting the experimental data to extract physical quantities in the next section. Fitting will be elaborated on by the two components in the future study.

Observables and parameter extraction
With the simple analytic formula, we can estimate the event rate evolution with the water-Cherenkov detector like Super-Kamiokande. Here, we focus only on the second component of the previous section, which dominates the late-time properties and is useful to extract physical parameters. The event rate is approximately given by the total number of protons in the detector, the number of anti-electron type neutrinos coming into the detector, and the cross section of inverse beta decay that is the main interaction capturing neutrinos. The event rate is given by ≈ 210 s −1 M det 32.5 kton where M det is the detector mass (32.5 kton corresponds to the entire volume of the inner tank of Super-Kamiokande), m is the nucleon mass, D is the distance between an SN and the 9/12 earth, σ is the cross section of inverse beta decay (p +ν e → n + e + ). For the cross section, we use σ(E ν ) = σ 0 (E ν /MeV) 2 with σ 0 = 9.4 × 10 −44 cm 2 [22] and by assuming the thermal spectrum σ = ( The average energy of positrons is given by Note that the positron energy given above is valid only when the typical energy of positrons is higher enough than the threshold energy of data analysis (see Section 4.2 of Ref. [12]). Once we would detect neutrinos from the next nearby SN, these formulae can be applied to narrow down the parameter space, which would give a starting point of more detailed calculations. For instance, by dividing Eq. (54) by the fifth power of Eq. (55), one finds R PNS = 10 km R 720 s −1 which is time independent. Next, by taking maximum of tR one finds which gives t 0 . Since M det is given by experiments and D would be measured by the optical or infrared observations, unknowns are M PNS , g and β. Unfortunately, at the current moment they are degenerate and the combination M PNS (gβ) 2/3 is only measurable. The degeneracy would be resolved by combining all the available data, for instance, the total energy emitted by neutrinos gives information of GM 2 PNS /R PNS (see Eq. 8). Also, a consistency relation of the analytic formulae is given by where the dot denotes the time derivative. By measuring it the model consistency can be tested.

Discussion
We discuss here caveats of our drastic assumptions to make the simple formulae. The assumptions of the analytic solutions are the following. We employ the Lane-Emden solution with n = 1 for the density and temperature profiles of neutron stars and assume the constant entropy profiles. For the neutrino transfer equation, we employ the thermal equilibrium (i.e., the Fermi-Dirac distribution). Also, the spherical symmetry is applied. These assumptions lead to caveats as follows: First, the density profiles are dependent on the nuclear equation of state so that we need a more realistic equation of state for the detailed evolution, which makes the analytic treatment difficult. Second, the entropy profile is not constant just after the explosion, but in the late phase, the neutrino diffusion 10/12 produces almost the constant entropy profiles [10,11]. Therefore, at least for the late time, a constant entropy profile is not a wrong assumption. Third, PNS convection may change the neutrino luminosity in the early time, but it becomes weak in the late time so that our analytic formula would not be affected. Forth, neutrino's spectrum is not purely thermal when the temperature gets low, and the diffusion timescale becomes shorter than the emission/absorption timescales, in the very late time. To give more realistic luminosity and spectrum, we need to solve the multi-energy neutrino transfer equation, which is the next step. Lastly, neutrino oscillation is not included in this study. But it is known that, in the late time, all flavors of neutrinos have similar luminosity and spectrum so that the neutrino oscillation does not change the analytic solutions [12].

Summary
In this paper, we derive analytic solutions of the neutrino radiation transfer equation in protoneutron star cooling after the core-collapse supernovae. The luminosity evolution is given by Eq. (47) and average energy evolution is given by Eq. (49). The evolution of the detection rate and the positron energies, for the water-Cherenkov detectors, are also given by Eqs. (54) and (55). The timescale in these equations is Eq. (48). With these equations, the neutrino emissions for the late time, in particular 1 s after the explosion, are given.
The analytic solutions presented in this paper are the very first step of analytic expression of detectable neutrino light curves so that formulae will be updated to give more useful expressions for the next galactic supernova in the forthcoming papers.