Neutrino Oscillations Within the Induced Gravitational Collapse Paradigm of Long Gamma-Ray Bursts

The induced gravitational collapse (IGC) paradigm of long gamma-ray bursts (GRBs) associated with supernovae (SNe) predicts a copious neutrino-antineutrino ($\nu\bar{\nu}$) emission owing to the hypercritical accretion process of SN ejecta onto a neutron star (NS) binary companion. The neutrino emission can reach luminosities of up to $10^{57}$MeV s$^{-1}$, mean neutrino energies 20 MeV, and neutrino densities $10^{31}$ cm$^{-3}$. Along their path from the vicinity of the NS surface outward, such neutrinos experience flavor transformations dictated by the neutrino to electron density ratio. We determine the neutrino and electron on the accretion zone and use them to compute the neutrino flavor evolution. For normal and inverted neutrino-mass hierarchies and within the two-flavor formalism ($\nu_{e}\nu_{x}$), we estimate the final electronic and non-electronic neutrino content after two oscillation processes: (1) neutrino collective effects due to neutrino self-interactions where the neutrino density dominates and, (2) the Mikheyev-Smirnov-Wolfenstein (MSW) effect, where the electron density dominates. We find that the final neutrino content is composed by $\sim$55% ($\sim$62%) of electronic neutrinos, i.e. $\nu_{e}+\bar{\nu}_{e}$, for the normal (inverted) neutrino-mass hierarchy. The results of this work are the first step toward the characterization of a novel source of astrophysical MeV-neutrinos in addition to core-collapse SNe and, as such, deserve further attention.


Introduction
The emergent picture of gamma-ray bursts (GRBs) is that both short-duration and long-duration GRBs originate from binary systems (Ruffini et al. 2016b).
Short bursts originate from neutron star-neutron star (NS-NS) or neutron star-black hole (NS-BH) mergers (see, e.g., Goodman 1986;Paczynski 1986;Eichler et al. 1989;Narayan et al. 1991). For this case Narayan et al. (1992) introduced the role of neutrino-antineutrino (nn ) annihilation leading to the formation of an electron-positron plasma ( -+ e e ) in NS-NS and NS-BH mergers. Such a result triggered many theoretical works, including the general relativistic treatment by Salmonson & Wilson (2002) of the nn annihilation process giving rise to the -+ e e plasma in an NS-NS system. For long bursts we stand on the induced gravitational collapse (IGC) paradigm (Ruffini et al. 2006(Ruffini et al. , 2008(Ruffini et al. , 2015bIzzo et al. 2012;Fryer et al. 2014), based on the hypercritical accretion process of the supernova (SN) ejecta of the explosion of a carbon-oxygen core (CO core ) onto an NS binary companion. In the above processes, the emission of neutrinos is a key ingredient.
We focus hereafter on the neutrino emission of long bursts within the IGC scenario. The role of neutrinos in this paradigm has been recently addressed by Fryer et al. (2014Fryer et al. ( , 2015 and Becerra et al. (2015Becerra et al. ( , 2016. The hypercritical accretion of the SN ejecta onto the NS companion can reach very high rates of up to 10 −2 M e s −1 , and its duration can be of the order of 10-10 4 s depending on the binary parameters. The photons become trapped within the accretion flow and thus do not serve as an energy sink. The high temperature developed on the NS surface leads to e − e + pairs that, via weak interactions, annihilate into nn pairs with neutrino luminosities of up to 10 52 erg s −1 for the highest accretion rates. Thus, this process dominates the cooling and gives rise to a very efficient conversion of the gravitational energy gained by accretion into radiation. We refer to Becerra et al. (2016) for further details on this process.
The above hypercritical accretion process can lead the NS to two alternative fates, leading to the existence of two long GRB subclasses (Fryer et al. 2014Becerra et al. 2015Becerra et al. , 2016Ruffini et al. 2016b): I. The hypercritical accretion leads to a more massive NS companion but not to a BH. These binaries explain the X-ray flashes (XRFs): long bursts with isotropic energy E iso 10 52 erg and rest-frame spectral peak energy E p,i 200 keV (see Ruffini et al. 2016b, for further details). The local observed number density rate of this GRB subclass is (Ruffini et  . II. The hypercritical accretion is high enough to make the NS reach its critical mass, triggering its gravitational collapse with consequent BH formation. These binaries explain the binary-driven hypernovae (BdHNe): long bursts with E iso 10 52 erg and E p,i 200 keV (see Ruffini et al. 2016b, for further details). The local observed number density rate of this GRB subclass is (Ruffini et al. 2016b . Simulations of the hypercritical accretion process in the above binaries have been presented in Fryer et al. (2014Fryer et al. ( , 2015 and Becerra et al. (2015Becerra et al. ( , 2016. It has been shown how, thanks to the development of a copious neutrino emission near the NS surface, the NS is allowed to accrete matter from the SN at very high rates. The specific conditions leading to XRFs and BdHNe, as well as a detailed analysis of the neutrino production in these systems, have been presented in Becerra et al. (2016). Neutrino emission can reach luminosities of 10 52 ergs −1 and the mean neutrino energy of the order of 20MeV. Under these conditions, XRFs and BdHNe become astrophysical laboratories for MeV-neutrino physics additional to core-collapse SNe.
On the other hand, the emission of TeV-PeV neutrinos is relevant for the observations of detectors such as the IceCube (Aartsen et al. 2013). High-energy neutrino emission mechanisms have been proposed within the context of the traditional model of long GRBs. In the traditional "collapsar" scenario (Woosley 1993;Paczyński 1998;MacFadyen & Woosley 1999) the gravitational collapse of a single, fast-rotating, massive star originates a BH surrounded by a massive accretion disk (see, e.g., Piran 2004, for a review), and the GRB dynamics follows the "fireball" model that assumes the existence of an ultrarelativistic collimated jet with Lorentz factor Γ ∼ 10 2 -10 3 (see, e.g., Shemi & Piran 1990;Meszaros et al. 1993;Piran et al. 1993;Mao & Yi 1994). This scenario has been adopted for the explanation of the prompt emission, as well as both the afterglow and the GeV emission of long GRBs. The GRB light-curve structures are there described by (internal or external) shocks (see, e.g., Rees & Meszaros 1992. The high-energy neutrinos in this context are produced from the interaction of shock-accelerated cosmic rays (e.g., protons) with the interstellar medium (see, e.g., Agostini et al. 2017;Kumar & Zhang 2015, and references therein). A recent analysis of the thermal emission of the X-ray flares observed in the early afterglow of long GRBs (at source rest-frame times t ∼ 10 2 s) shows that it occurs at radii ∼10 12 cm and expands with a mildly relativistic Γ4 (see Ruffini et al. 2018, for further details). This rules out the ultrarelativistic expansion in the GRB afterglow traditionally adopted in the literature. Interestingly, the aforementioned mechanisms of high-energy neutrino production conceived in the collapsar-fireball model can still be relevant in the context of BdHNe and authentic short GRBs (S-GRBs, NS-NS mergers with E iso 10 52 erg leading to BH formation; see Ruffini et al. 2016b, for the classification of long and short bursts in seven different subclasses). The emission in the 0.1-100 GeV energy band observed in these two GRB subclasses has been shown to be well explained by a subsequent accretion process onto the newly born BH (Ruffini et al. 2015a(Ruffini et al. , 2015b(Ruffini et al. , 2016a(Ruffini et al. , 2016bAimuratov et al. 2017; see also Aimuratov et al., in preparation). Such GeV emission is not causally connected either with the prompt emission or with the afterglow emission composing the flaring activity (Ruffini et al. 2018). An ultrarelativistic expanding component is therefore expected to occur in BdHNe and S-GRBs, which deserves to be explored in forthcoming studies as a possible source of high-energy neutrinos. Specifically, this motivates the present article to identify the possible additional channels to be explored in the hypercritical accretion not around an NS but around a BH.
The aim of this article is to extend the analysis of the MeVneutrino emission in the hypercritical accretion process around an NS in the XRFs and BdHNe to assess the possible occurrence of neutrino flavor oscillations.
We shall show in this work that, before escaping to the outer space, i.e., outside the Bondi-Hoyle accretion region, the neutrinos experience an interesting phenomenology. The neutrino density near the NS surface is so high that the neutrino self-interaction potential, usually negligible in other very well known scenarios such as the Sun, the upper layers of Earth's atmosphere, and terrestrial reactor and accelerator experiments, becomes more relevant than the matter potential responsible for the famous Mikheyev-Smirnov-Wolfenstein (MSW) effect (Wolfenstein 1978;Mikheev & Smirnov 1986). A number of papers have been dedicated to the consequences of the neutrino self-interaction dominance (Notzold & Raffelt 1988;Pantaleone 1992;Qian & Fuller 1995;Pastor & Raffelt 2002;Sawyer 2005;Duan et al. 2006bDuan et al. , 2007Duan et al. , 2008aDuan et al. , 2008bDuan et al. , 2010Fuller & Qian 2006;Esteban-Pretel et al. 2007, 2008Fogli et al. 2007;Chakraborty et al. 2008;Dasgupta & Dighe 2008;Dasgupta et al. 2008a;Sawyer 2009;Wu & Qian 2011), most of them focused on SN neutrinos. In these cases, the SN induces the appearance of collective effects such as synchronized and bipolar oscillations leading to an entirely new flavor content of emitted neutrinos when compared with the spectrum created deep inside the star. The density of neutrinos produced in the hypercritical accretion process of XRFs and BdHNe is such that the neutrino self-interactions, as in the case of SNe, dominate the neutrino flavor evolution, giving rise to the aforementioned collective effects. The main neutrino source, in this case, is the nn pair production from e − e + annihilation (Becerra et al. 2016), which leads to an equal number of neutrinos and antineutrinos of each type. This equality does not happen in the SN standard scenario. We will show that bipolar oscillations, inducing very quick flavor pair conversions n n n n n n « « m m t t¯ē e , can occur with oscillation length as small as O(0.05-1) kilometers. However, the n n -¯symmetry characterizing our system leads to the occurrence of kinematic decoherence, making the neutrino flavor content reach equipartition deep inside the accretion zone. In the regions far from the NS surface where the neutrino density is not so high, the matter potential begins to dominate and MSW resonances can take place. As a result, an entirely different neutrino flavor content emerges from the Bondi-Hoyle surface when compared with what was originally created in the bottom of the accretion zone.
This article is organized as follows. In Section 2 we outline the general features of the accretion process onto the NS within the IGC paradigm and present the processes responsible for the neutrino creation. From these features, we obtain the distribution functions that describe the neutrino spectrum near the NS surface. Section 3 shows a derivation of the equations that drive the evolution of neutrino oscillations closely related to the geometrical and physical characteristics of our system. We discuss some details on the neutrino oscillation phenomenology. Since we have to face a nonlinear integro-differential system of equations of motion (EoM), we introduce the singleangle approximation to later recover the full realistic phenomenology after generalizing our results to the multiangle approach and, consequently, de-coherent picture. In Section 5 the final neutrino emission spectra are presented and compared with those in which neutrinos are created in the accretion zone. Finally, we present in Section 6 the conclusions and some perspectives for future research on this subject.

Neutrino Creation during Hypercritical Accretion
The SN material first reaches the gravitational capture region of the NS companion, namely, the Bondi-Hoyle region. The infalling material shocks as it piles up onto the NS surface, forming an accretion zone where it compresses and eventually becomes sufficiently hot to trigger a highly efficient neutrino emission process. Neutrinos take away most of the infalling matter's gravitational energy gain, letting it reduce its entropy and be incorporated into the NS. Figure 1 shows a sketch of this entire hypercritical accretion process.
It was shown in Becerra et al. (2016) that the matter in the accretion zone near the NS surface develops conditions of temperature and density such that it is in a nondegenerate, relativistic, hot plasma state. The most efficient neutrino emission channel under those conditions becomes the electron-positron pair annihilation process: The neutrino emissivity produced by this process is proportional to the accretion rate to the 9/4 power (see below). This implies that the higher the accretion rate, the higher the neutrino flux; hence, the largest neutrino flux occurs at the largest accretion rate.
We turn now to estimating the accretion rate and thus the neutrino emissivity we expect in our systems.

Accretion Rate in XRFs and BdHNe
We first discuss the amount of SN matter per unit time reaching the gravitational capture region of the NS companion, namely, the Bondi-Hoyle accretion rate. It has been shown in Bayless et al. (2015) and Becerra et al. (2016) that the shorter (smaller) the orbital period (separation), the higher the peak accretion rateṀ peak and the shorter the time at which it peaks, t peak .
The Bondi-Hoyle accretion rate is proportional to the density of the accreted matter and inversely proportional to its velocity. Thus, we expect the accretion rate to increase as the denser and slower inner layers of the SN reach the accretion region. Based on these arguments, Becerra et al. (2016) derived simple, analytic formulae forṀ peak and t peak as a function of the orbital period (given all the other binary parameters) that catch both the qualitatively and quantitatively behaviors of these two quantities obtained from full numerical integration. We refer the reader to Appendix A of that article for further details. For the scope of this work these analytic expressions are sufficient to give us an estimate of the hypercritical accretion rates and related timescale developed in these systems: where P is the orbital period, m is the index of the power-law density profile of the pre-SN envelope, v star,0 is the velocity of the outermost layer of the SN ejecta, M = M CO + M NS is the total binary mass, and M CO = M env + M νNS is the total mass of the CO core given by the envelope mass and the mass of the central remnant, i.e., the new NS (hereafter νNS) formed from the region of the CO core that undergoes core collapse (i.e., roughly speaking the iron core of density ρ core and radius R core ). We here adopt = where R star 0 is the total radius of the pre-SN CO core , r core and R core are parameters of the pre-SN density profile introduced to account of the finite size of the envelope, and m is the powerlaw index followed by the density profile at radii r > R core (see Becerra et al. 2016, for further details). Figure 2 shows the peak accretion rate in Equation (2) as a function of the orbital period. In this example, we consider the following binary parameters (see Becerra et al. 2016, for details): a CO core produced by a zero-age main-sequence (ZAMS) progenitor with M ZAMS = 20 M e , i.e., M CO = 5.4 M e , an initial NS mass of 2.0 M e , and a velocity of the outermost ejecta layer v star,0 = 2 × 10 9 cm s −1 . For these parameters, η ≈ 0.41. Becerra et al. (2015Becerra et al. ( , 2016 showed the existence of a maximum orbital period, P max , over which the accretion onto an NS companion is not high enough to bring it to the critical mass for gravitational collapse to a BH. As mentioned in the Introduction, CO core -NS binaries with P > P max lead to XRFs, while the ones with PP max lead to BdHNe. For the binary parameters of the example in Figure 2, P max ≈ 127 minutes (vertical dashed line). We can therefore conclude that BdHNe can have peak accretion rates in the range~-

Neutrino Emission at Maximum Accretion
For the accretion rate conditions characteristic of our models at peak ∼10 −4 to 10 −2 M e s −1 , pair annihilation dominates the neutrino emission and electron neutrinos remove the bulk of the energy (Becerra et al. 2016). The e + e − pairs producing the neutrinos are thermalized at the matter temperature. This temperature is approximately given by where P shock is the pressure of the shock developed on the accretion zone above the NS surface,Ṁ acc is the accretion rate, v acc is the velocity of the infalling material, σ is the Stefan-Boltzmann constant, and c is the speed of light. It can be checked that, for the above accretion rates, the system develops temperatures and densities (T10 10 K and ρ10 6 g cm −3 ; see, e.g., Figure16 in Becerra et al. 2016) for which the neutrino emissivity of the e + e − annihilation process can be estimated by the simple formula (Yakovlev et al. 2001) where k B is the Boltzmann constant. The accretion zone is characterized by a temperature gradient with a typical scale height Δr ER = T/∇T ≈ 0.7 R NS . Owing to the strong dependence of the neutrino emission on temperature, most of the neutrinos are emitted from a spherical shell around the NS of thickness (see Figure 1) 0.08 . 6 e e e e

ER NS
Equations (4) and (5) imply that the neutrino emissivity satisfies  µ -+Ṁ e e acc 9 4 as we had anticipated. These conditions lead the neutrinos to be efficient in balancing the gravitational potential energy gain, allowing the hypercritical accretion rates.
The effective accretion onto the NS can be estimated as where ΔM ν and L ν are, respectively, the mass and neutrino luminosity in the emission region, and = E g D +D , 8 e e NS 2 with  -+ e e being the neutrino emissivity in Equation (5). For M NS = 2 M e and temperatures 1-10 MeV, Equations (7) and (8) give » -M 10 eff 10 to 10 −1 M e s −1 and L ν ≈ 10 48 -10 57 MeV s −1 .

Neutrino Spectrum at the NS Surface
After discussing the general features of neutrino emission during the accretion process, it is necessary for our analysis of the neutrino oscillations to determine the neutrino spectrum at the NS surface. Specifically, we need to determine the ratios at which the neutrinos of each flavor are created and their average energy so that we can find a fitting distribution function f ν with these characteristics.
Since the main source of neutrinos is the e − e + pair annihilation process, we can conclude that neutrinos and antineutrinos are created in equal number. Furthermore, the information about the neutrino and antineutrino emission of a given flavor i can be calculated from the integral (Yakovlev et al. 2001 where h  e is the electron (positron) degeneracy parameter, including its rest mass. The Dicus cross section σ i is written in terms of the electron and positron four-momenta = as (Dicus 1972)  The factors  C i , 2 are written in terms of the weak interaction vector and axial-vector constants:  . Peak accretion rate,Ṁ peak , as a function of the binary orbital period, as given by Equation (2). This example corresponds to the following binary parameters: a CO core formed by an M ZAMS =20M e progenitor, i.e., M CO =5.4M e , an initial NS mass of 2.0M e , v star,0 =2×10 9 cms −1 , η≈0.41, and index m=2.946 (see Becerra et al. 2016, for further details). For these parameters the largest orbital period for the induced collapse of the NS to a BH by accretion is P max ≈127 minutes, which is represented by the vertical dashed line.
For m=0 and m=1 Equation (9) gives the neutrino and antineutrino number emissivity (neutrino production rate) and the neutrino and antineutrino energy emissivity (energy per unit volume per unit time) for a certain flavor i, respectively. Hence, not only are we able to calculate the total number emissivity with but we can also calculate the neutrino or antineutrino energy moments with We wish to construct a Fermi-Dirac-like fitting formula for the neutrino spectrum as is usually done in supernova neutrino emission (Janka & Hillebrandt 1989a, 1989b, that is, a function like Equation (10) in terms of two parameters: the effective neutrino temperature nn T and the effective neutrino degeneracy parameter h nn , otherwise known as the pinching parameter (Raffelt 1996;Keil et al. 2003). To that end, it is enough to calculate the first two moments. In particular, for a relativistic nondegenerate plasma Table 1) Equation (9) can be approximated with a very good accuracy by (Yakovlev et al. 2001 , and adding over every flavor, this expression reduces to Equation (5). With Equations (13) and (14) we find , . However, as in Equation (16), we will simplify our description using only two flavors: the electronic neutrinos and antineutrinos n n , e e , and a superposition of the other flavors n n , . This can be understood as follows. Since the matter in the accretion zone is composed by protons, neutrons, electrons, and positrons, ν e and n e interact with matter by both charged and neutral currents, while n m , n t , n m , and n t interact only by neutral currents. Therefore, the behavior of these states can be clearly divided into electronic and nonelectronic. This distinction will come in handy when studying neutrino oscillations. F c i i , respectively, and using Equation (16), we can recollect two important facts: , , 17a  (15), from here on out we will use the approximation where p is the neutrino momentum. 4. From Equation (13) we obtain the same energy moments for both neutrinos and antineutrinos, but, as Misiaszek et al. (2006) points out, these energies should be different since, in reality, this expression returns the arithmetic mean of the particle and antiparticle energy moments, that is, . However, if we calculate the differences between the energy moments with Equations (41) and (46)  2 . These differences are small enough that we can use the same effective temperature and pinching parameter for both neutrinos and antineutrinos.

Solving the equations
for any value of T in Table 1 , and then we multiply by p e = á ñá ñ = D n n n n n n n n n ( ) It can be checked that these distributions obey and with these conditions satisfied we can conclude that Equations (21) are precisely the ones that emulate the neutrino spectrum at the NS surface. In Table 1 we have collected the values of every important quantity used in the calculations within this section for the range of accretion rates in which we are interested.
Considering that the problem we attacked in this section reduces to finding a normalized distribution whose first two moments are fixed, the choice we have made with Equation (21) is not unique. The solution depends on how many moments are used to fit the distribution and what kind of function is used as an ansatz. A different solution based on a Maxwell-Boltzmann distribution can be found in Keil et al. (2003), Fogli et al. (2005), and Misiaszek et al. (2006).
At this stage, we can identify two main differences between neutrino emission in SNe and in the IGC process of XRFs and BdHNe, within the context of neutrino oscillations. The significance of these differences will become clearer in the following sections, but we mention them here to establish a point of comparison between the two systems since SN neutrino oscillations have been extensively studied. x . It has been argued that this difference between neutrinos and antineutrinos is enough to dampen kinematical decoherence, so that bipolar oscillations are a feature present in SN neutrinos (see, e.g., Esteban-Pretel et al. 2007).
In the next section, we will use the results presented here to determine the neutrino flavor evolution in the accretion zone.
As we discussed in Section 2, in our physical system of interest neutrinos are mainly created by electron-positron pair annihilation, and so the number of neutrinos is equal to the number of antineutrinos. Such a fact creates an interesting and unique physical situation, different from, for example, SN neutrinos, for which traditional models predict a predominance of electron neutrinos mainly due to the deleptonization caused by the URCA process (see, e.g., Esteban-Pretel et al. 2007).
The neutrino self-interaction potential decays with the radial distance from the NS faster than the matter potential. This is a direct consequence of the usual 1/r 2 flux dilution and the collinearity effects due to the neutrino velocity dependence of the potential. Consequently, we identify three different regions along the neutrino trajectory in which the oscillations are dominated by intrinsically different neutrino phenomenology. Figure 3 illustrates the typical situation of the physical system we are analyzing. Just after the neutrino creation in the regions of the accretion zone very close to the surface of the NS, neutrinos undergo kinematic decoherence along the same length scale of a single cycle of the so-called bipolar oscillations. Bipolar oscillations imply very fast flavor conversion between neutrino pairs n n n n n n « « m m t t¯ē e , and, amazingly, the oscillation length in this region can be as small as of the order of tens of meters. Note that kinematic decoherence is just the averaging over flavor neutrino state processes resulting from quick flavor conversion, for which oscillation length depends on the neutrino energy. It does not imply quantum decoherence, and thus, neutrinos are yet able to quantum oscillate if appropriate conditions are satisfied. In fact, as can be observed from Figures 4 and 5, bipolar oscillations preserve the characteristic oscillation pattern, differently from quantum decoherence, which would lead to a monotonous dumping figure.
Kinematic decoherence is relevant when three conditions are met: (i) the self-interaction potential dominates over the vacuum potential, (ii) the matter potential does not fulfill the MSW condition, and (iii) there is a low asymmetry between the neutrino and antineutrino fluxes. We will see that our system satisfies all three conditions.
As the self-interaction potential becomes small and the matter potential becomes important, oscillations are suppressed and we do not expect significant changes in the neutrino flavor content along this region. This situation changes radically when the matter potential is so small that it is comparable to neutrino vacuum frequencies Δm 2 /2p, where Δm 2 is the neutrino squared mass difference and p is the norm of the neutrino momentum p. In this region, the neutrino self-interaction potential is negligible and the usual MSW resonances can occur. Therefore, we can expect a change in the neutrino spectrum.
We dedicate this section to a detailed derivation of the EoM of flavor evolution. In later sections, we will analyze the neutrino oscillation phenomenology to build the neutrino emission spectrum from a binary hyperaccretion system.

Equations of Motion
The EoM that govern the evolution of an ensemble of mixed neutrinos are the quantum Liouville equations Taking into account the current-current nature of the weak interaction in the standard model, the Hamiltonian for each equation is (Dolgov 1981;Sigl & Raffelt 1993;Hannestad et al. 2006) ò ò p r r p where Ω p is the matrix of vacuum oscillation frequencies, l p andl p are matrices of occupation numbers for charged leptons built in a similar way to the neutrino matrices, and = v p p p is the velocity of a particle with momentum p (either neutrino or charged lepton).
As in Section 2, we will only consider two neutrino flavors: e and x=μ+τ. Three-flavor oscillations can be approximated to two-flavor oscillations as a result of the strong hierarchy of the squared mass differences D 2 12 2 (see Table 2). In this case, only the smallest mixing angle θ 13 is considered. We will drop the suffix for the rest of the discussion. Consequently, the relevant oscillations are n n  e x and n n ē x , and each term in the Hamiltonian governing oscillations becomes a 2×2 Hermitian matrix.
Let us first present the relevant equations for neutrinos. Due to the similarity between H p andH p , the corresponding equations for antineutrinos can be obtained in an analogous manner. In the two-flavor approximation, ρ in Equation (23) can be written in terms of Pauli matrices and the polarization vector P p as Figure 3. Interaction potentials as functions of the radial distance from the NS center for selected accretion ratesṀ (see Table 1). Each plot runs from the NS surface to the Bondi-Hoyle surface. m r stands for the self-interaction neutrino potential, λ r is the matter potential, and ω H and ω L are the higher and lower resonances corresponding to the atmospheric and solar neutrino scales, respectively, defined in Equation (59). Outside the Bondi-Hoyle region the neutrino and electron densities depend on the direction of their path relative to the SN and the particular ejecta density profile.
sin 2 , 0, cos 2 , and θ is the smallest neutrino mixing angle in vacuum.
The other two terms in Equation (24) are special since they make the evolution equations nonlinear. Even though they are very similar, we are considering that the electrons during the accretion form an isotropic gas; hence, the vector v q in the first integral is distributed uniformly on the unit sphere and the factor · v v q p averages to zero. After integrating, the matter Hamiltonian is given by is the charged current matter potential and = ( ) L 0, 0, 1 . Such simplification cannot be made with the final term. Since neutrinos are responsible for the energy loss of the infalling material during accretion, they must be escaping the accretion zone and the net neutrino and antineutrino flux is nonzero. In this case the factor · v v q p cannot be averaged to zero. At any rate, we can still use Equation (25) and obtain (Pantaleone 1992;Malkus et al. 2016;Zhu et al. 2016) (23), and using the commutation relations of the Pauli matrices, we find the EoM for neutrinos and antineutrinos for each momentum

Introducing every Hamiltonian term in Equation
Solving the above equations would yield the polarization vectors as a function of time. However, in our specific physical system, both the matter potential λ and the neutrino potential vary with the radial distance from the NS surface, as well as the instant t of the physical process, which can be characterized by the accretion rateṀ. As we will see later, the time dependence can be ignored. This means that Equation (32) must be written in a way that makes explicit the spatial dependence, i.e., in terms of the vector r. For an isotropic and homogeneous neutrino gas or a collimated ray of neutrinos the expression dt=dr would be good enough, but for radiating extended sources the situation is more complicated. In Equation (23) we must replace the matrices of occupation numbers by the spacedependent Wigner functions r p r , (and r p r , ) and the total time derivative by the Liouville operator (Cardall 2008 We will ignore the third term of the Liouville operator since we would not consider the gravitational deflection of neutrinos. The distances traveled by a neutrino in these times are r≈3×10 12 -3×10 18 cm. These distances are much larger than the typical binary separation a. As a consequence,  we can consider the neutrino evolution to be a stationary process. This fact allows us to neglect the first term in Equation (33). Putting together these results, the EoM become where H p r , andH p r , are the same as in Equation (24) but the matrices of densities (as well as the polarization vectors) depend on the position r. Note, however, that the electrons in the accretion zone still form an isotropic gas, Equation (30)  . The first two terms in the Hamiltonian remain virtually unchanged. On the other hand, projecting the EoM onto the radial distance from the NS and using the axial symmetry of the system, the integral in the neutrino-neutrino interaction term can be written as Since farther from the NS the interacting neutrinos approach a perfect collinearity, the projected velocities J v r become decreasing functions of the position. In this particular geometry the diagonal elements of the matrix of densities are written as a product of independent distributions over each variable J f p, , , where the f dependence has been integrated out. The one over p is the normalized Fermi-Dirac distribution, and the one over ϑ is assumed uniform owing to symmetry. The r dependence is obtained through the geometrical flux dilution. Knowing this, the diagonal elements of matrices of densities at the NS surface are where the functions n f i are given by Equation (21).

Single-angle Approximations
The integro-differential Equations (32) and (34) are usually numerically solved for the momentum p and the scalar · v v q p . Such simulations are quite time-consuming, and the result is frequently too complicated to allow for a clear interpretation of the underlying physics. For this reason, the analytic approximation called the single-angle limit is made. Such an approximation consists of imposing a self-maintained coherence in the neutrino system, i.e., it is assumed that the flavor evolution of all neutrinos emitted from an extended source is the same as the flavor evolution of the neutrinos emitted from the source along a particular path. Under this premise, the propagation angle between the test neutrino and the background neutrinos is fixed. In expression (35) this is equivalent to dropping the J¢ dependence of ρ and replacing the projected velocity J v r either by an appropriate average at each r (as in Dasgupta & Dighe 2008) or by a representative angle (usually 0 or π/4). We will follow the former approach and apply the bulb model described in Duan et al. (2006a). Within this model it is shown that the projected velocity at a distance r from the neutrino emission zone is where v R NS is the projected velocity at the NS surface. By redefining the matrices of density with a change of variable , , and using Equation (25), we can write the full EoM as where we have replaced v r by its average value All the interaction potentials now depend on r, and each effective potential strength is parameterized as follows (Dasgupta & Dighe 2008): It is worth mentioning that all the effective potential strengths are affected by the geometry of the extended source through the projected velocity on the right-hand side of Equation (34). For the neutrino-neutrino interaction potential, we have chosen the total neutrino number density as parameterization. This factor comes from the freedom to renormalize the polarization vectors in the EoM. A different choice has been made in Esteban-Pretel et al. (2007). Of the other two r-dependent factors, one comes from the geometrical flux dilution and the other accounts for collinearity in the single-angle approximation. Overall μ r decays as r 1 4 . In Figure 3 the behavior of the effective potentials within the single-angle formalism is shown for = -M 10 2 , 10 −4 , 10 −6 , and 10 −8 M e s −1 . In all cases, the neutrino energy is the corresponding average reported in Table 1. Since the oscillatory dynamics of the neutrino flavors are determined by the value of the potentials, and the value of the potentials depends on the data in Table 1, it is important to establish how sensible this information is to the model we have adopted, in particular, to the pre-SN envelope density profile index m. The reported accretion rates can be seen as different states in the evolution of a binary system or as peak accretion rates of different binary systems. For a given accretion rate, the temperature and density conditions on the NS surface are fixed. This, in turn, fixes the potentials involved in the equations of flavor evolution and the initial neutrino and antineutrino flavor content. To see the consequences of changing the index m, we can estimate the peak accretion rates for new values using Equation (2). Since we are only interested in SNe Ic, we shall restrict these values to the ones reported in Table1 of Becerra et al. (2016) (that is, m = 2.771, 2.946, and 2.801), and in each case, we consider the smallest binary separation such that there is no Roche lobe overflow. For these parameters, we find peak accretion ratesM peak 1 with peak times at t peak ≈7-35 minutes. Because these accretion rates are still within the range in Table 1, the results contained in Section 4 apply also to these cases with different values of the m-index.
The profiles for the electron and positron number densities were adopted from the simulations presented in Becerra et al. (2016). Due to the dynamics of the infalling matter, close to the NS, the behavior of --+ ( ) ( ) n r n r e e is similar to μ r . At the shock radius, the electron density's derivative presents a discontinuity, and its behavior changes, allowing for three distinct regions inside the Bondi-Hoyle radius. The matter potential is always higher than the neutrino potential, yet, in most cases, both are higher than the vacuum potential, so we expect neutrino collective effects (neutrino oscillations) and MSW resonances to play a role in the neutrino flavor evolution inside the Bondi-Hoyle radius. Outside the capture region, as long as the neutrinos are not directed toward the SN, they will be subjected to vacuum oscillations.

Single-angle Solutions and Multi-angle Effects
The full dynamics of neutrino oscillations is a rather complex interplay between the three potentials discussed in Section 3, yet the neutrino-antineutrino symmetry allows us to generalize our single-angle calculations for certain accretion rates using some numerical and algebraic results obtained in Hannestad et al. (2006), Fogli et al. (2007, Esteban-Pretel et al. (2007), and references therein. Specifically, we know that if μ r ?ω r , as long as the MSW condition λ r ;ω r is not met, collective effects should dominate the neutrino evolution even if λ r ?μ r . On the other hand, if μ r ω r , the neutrino evolution is driven by the relative values between the matter and vacuum potentials. With this in mind, we identify two different ranges of values for the accretion rate:

High Accretion Rates
For accretion rates ´- hence, we expect strong effects of neutrino self-interactions. In order to appreciate the interesting physical processes that happen with the neutrinos along their trajectory in the accretion zone, we begin this analysis with a simplified approach to the EoM for a monochromatic spectrum with the same energy for both neutrinos and antineutrinos. Let us introduce the following definitions: The role of the matter potential is to logarithmically extend the period of the bipolar oscillations, so we can ignore it for now. Also, we will restrict our analysis to a small enough region at + D n R r NS so that we can consider w m » ( ) 0 d dr r r (adiabatic approximation). Then, by summing and subtracting Equation (39) and using definitions (45) and (46), we obtain We are now able to build a very useful analogy. The equations above are analogous to the EoM of a simple mechanical pendulum with a vector position given by Q, precessing around an angular momentum D, subjected to a force wmB with a moment of inertia proportional to the inverse of μ. With Equations (17) and (26) the initial conditions for the polarization vectors are

NS NS
We can easily show that , it can be checked that this value is conserved.
The analogous angular momentum is . Thus, the pendulum moves initially in a plane defined by B and the z-axis, i.e., the plane xz. Then, it is possible to define an angle j between Q and the z-axis such that j j = | |( ) ( ) Q Q sin , 0, cos . 50 Note that the only nonzero component of D is the y-component, and from Equations (47) and (48) we find j m = | | ( ) D d dr 51 and The above equations can be equivalently written as where we have introduced the inverse characteristic distance k by which is related to the anharmonic oscillations described by the nonlinear EoM (51) and (52). The logarithmic correction to the oscillation length due to matter effects is (Hannestad et al. 2006) The initial conditions (49) imply To investigate the physical meaning of the above equation, let us assume for a moment that q 2 is a small angle. In this case j ( ) R NS is also a small angle. If > k 0 2 , which is true for the normal hierarchy D > m 0 2 , we expect small oscillations around the initial position since the system begins in a stable position of the potential associated with Equations (51) and (52). No strong flavor oscillations are expected. On the contrary, for the inverted hierarchy Δm 2 <0, k 2 <0 and the initial j ( ) R NS indicates that the system begins in an unstable position, and we expect very large anharmonic oscillations. P z (as well as P z ) oscillates between two different maxima passing through a minimum P z ( P -¯z) several times. This behavior implies total flavor conversion: all electronic neutrinos (antineutrinos) are converted into nonelectronic neutrinos (antineutrinos) and vice versa. This has been called bipolar oscillations in the literature (Duan et al. 2010).
We solved numerically Equation (39) for both normal and inverted hierarchies using a monochromatic spectrum dominated by the average neutrino energy for = -M 10 ,  Table 1 with the initial conditions given by Equations (17) and (36). The behavior of the electronic neutrino survival probability inside the accretion zone is shown in Figures 4 and 5 for inverted hierarchy and normal hierarchy, respectively. For the inverted hierarchy, there is no difference between the neutrino and antineutrino survival probabilities. This should be expected since for these values of r the matter and self-interaction potentials are much larger than the vacuum potential, and there is virtually no difference from Equation (39). Also, note that the antineutrino flavor proportions discussed in Section 2.3 remain virtually unchanged for normal hierarchy, while the neutrino flavor proportions change drastically around the point λ r ∼ω r . The characteristic oscillation length of the survival probability found on these plots is t » -( ) 0.05 1 km, 57 which agrees with the ones given by Equation (55) calculated at the NS surface up to a factor of order 1. Such a small value of τ suggests extremely quick n n n n «ē e x x oscillations. Clearly, the full EoM are highly nonlinear, so the solution may not reflect the real neutrino flavor evolution. Concerning the single-angle approximation, it is discussed in Hannestad et al. (2006), , and Fogli et al. (2007) that in the more realistic multi-angle approach, kinematic decoherence happens. And in Esteban-Pretel et al. (2007) the conditions for decoherence as a function of the neutrino flavor asymmetry have been discussed. It is concluded that if the symmetry of neutrinos and antineutrinos is broken beyond the limit of O(25%), i.e., if the difference between emitted neutrinos and antineutrinos is roughly larger than 25% of the total number of neutrinos in the medium, decoherence becomes a subdominant effect.
As a direct consequence of the peculiar symmetric situation we are dealing with, in which neutrinos and antineutrinos are produced in similar numbers, bipolar oscillations happen, and, as we have already discussed, they present very small oscillation length as shown in Equation (57). Note also that the bipolar oscillation length depends on the neutrino energy. Therefore, the resulting process is equivalent to an averaging over the neutrino energy spectrum and an equipartition among different neutrino flavors is expected . Although, for simplicity, we are dealing with the two neutrino hypothesis, this behavior is easily extended to the more realistic three-neutrino situation. We assume, therefore, that at a few kilometers from the emission region neutrino flavor equipartition is a reality: Note that the multi-angle approach keeps the order of the characteristic length τ of Equation (55) unchanged and kinematics decoherence happens within a few oscillation cycles (Sawyer 2005;Hannestad et al. 2006;. Therefore, we expect that neutrinos created in regions close to the emission zone will be equally distributed among different flavors in less than a few kilometers after their creation. Once the neutrinos reach this maximally mixed state, no further changes are expected up until the matter potential enters the MSW resonance region. We emphasize that kinematics decoherence does not mean quantum decoherence. Figures 4 and 5 clearly show the typical oscillation pattern, which happens only if quantum coherence is still acting on the neutrino system. Differently from quantum decoherence, which would reveal itself by a monotonous dumping in the oscillation pattern, kinematics decoherence is just the result of averaging over the neutrino energy spectrum resulting from quick flavor conversion, for which oscillation length depends on the neutrino energy. Therefore, neutrinos are yet able to quantum oscillate if appropriate conditions are satisfied.
We discuss now the consequences of the matter potential.

Matter Effects
After leaving the emission region, beyond r≈R NS +Δr ν , where Δr ν is the width defined in Equation (6), the effective neutrino density quickly falls in an asymptotic behavior μ r ≈1/r 4 . The decay of λ r is slower. Hence, very soon the neutrino flavor evolution is determined by the matter potential. Matter suppresses neutrino oscillations, and we do not expect significant changes in the neutrino flavor content along a large region. Nevertheless, the matter potential can be so small that there will be a region along the neutrino trajectory in which it can be compared with the neutrino vacuum frequencies and the higher and lower resonant density conditions will be satisfied, i.e., where Δm 2 and Dm 21 2 are, respectively, the squared mass differences found in atmospheric and solar neutrino observations. Table 2  for m 3 < m 1 <m 2 . When the above resonance conditions are satisfied, the MSW effects happen and the flavor content of the flux of electronic neutrinos and antineutrinos will again be modified.
The final fluxes can be written as In order to evaluate n ( ) F E e and n ( ) F E e after matter effects, we have to estimate the survival probability at the resonant regions. There are several articles devoted to this issue; for instance, we can adopt the result in Fogli et al. (2003) The factor X, the conversion probability between neutrino physical eigenstates, is given by Petcov (1987), Fogli et al. (2003), and Kneller & McLaughlin (2006): The factor X is related to how fast physical environment features relevant for neutrino oscillations change, such as neutrino and matter densities.
For slow and adiabatic changes  X 0, while for fast and nonadiabatic changes  X 1. In our specific cases, the MSW resonances occur very far from the accretion zone, where the matter density varies very slow and therefore  X 0, as can be explicitly calculated from Equation (63). Consequently, it is straightforward to estimate the final fluxes of electronic and nonelectronic neutrinos and antineutrinos.

Low Accretion Rates
For accretion rates <´-- M M 5 10 s 5 1 , either the matter potential is close enough to the vacuum potential and the MSW condition is satisfied, or both the self-interaction and matter potentials are so low that the flavor oscillations are only due to the vacuum potential. In both cases, bipolar oscillations are not present. In Figure 6 we show the survival probability for = -- M M 10 s 6 1 as an example. We can see that neutrinos and antineutrinos follow different dynamics. In particular, for antineutrinos there are two decreases. The first one, around r≈(1-2) R NS , is due to bipolar oscillations that are rapidly damped by the matter potential as discussed in Section 4.1.1. The second one happens around r≈(10-20) R NS . It can be seen from the bottom left panel of Figure 3 (that one for = -- M M 10 s 6 1 ) that around r≈(1-2)× 10 7 cm (or, equivalently, r ≈ (10-20) R NS ) the higher MSW resonance occurs (l w r r H ). For inverted hierarchy, such resonance will affect antineutrinos depleting its number, as can be seen from Equation (60). Without bipolar oscillations, it is not possible to guarantee that decoherence will be complete and Equation (58) is no longer valid. The only way to know the exact flavor proportions is to solve the full Equation (32).

Neutrino Emission Spectra
Using the calculations of the previous section, we can draw a comparison between the creation spectra of neutrinos and antineutrinos at the NS surface ( n n F n , c c ), initial spectra after kinematic decoherence ( n n F n , 0 0 ), and emission spectra after the MSW resonances ( n n F n , ). Table 3 contains a summary of the flavor content inside the Bondi-Hoyle radius. With these fractions and Equation (21) Figure 7. In such figures, the left column corresponds to normal hierarchy and the right column corresponds to inverted hierarchy. The first two rows show the number fluxes after each process studied. The last row shows the relative fluxes Each column corresponds to a neutrino mass hierarchy: normal hierarchy on the left and inverted hierarchy on the right. The first two rows show the number fluxes after each process studied. n F C , n F 0 , and n F are the creation flux at the bottom accretion zone due to e + e − pair annihilation, the flux after the region with dominant neutrino-neutrino potential, and the final emission flux after the region with dominant neutrino-matter potential, respectively. The last row shows the relative fluxes n n F F C between the creation and emission fluxes. so that each one is a normalized Fermi-Dirac distribution multiplied by the appropriate flavor content fraction. To reproduce any other case, it is enough to use Equation (21) with the appropriate temperature.
At this point two comments have to be made about our results: 1. As we mentioned before, the fractions in Table 3 were obtained by assuming a monochromatic spectrum and using the single-angle approximation. This would imply that the spectrum-dependent phenomenon called the spectral stepwise swap of flavors is not present in our analysis even though it has been shown that it can also appear in multi-angle simulations (Fogli et al. 2007). Nevertheless, we know from our calculations in Section 2.3 that neutrinos and antineutrinos of all flavors are created with the exact same spectrum up to a multiplicative constant. Hence, following Raffelt & Smirnov (2007a, 2007b, by solving the equation we find that the critical (split) energy is E c =0. This means that the resulting spectrum should still be unimodal and the spectral swap in our system could be approximated by a multiplicative constant that is taken into account in the decoherence analysis of Section 4. 2. The fluxes of electronic neutrinos and antineutrinos shown in these figures and in Equation (60) represent fluxes at different positions up to a geometrical 1/r 2 factor, r being the distance from the NS radius. Also, since we are considering the fluxes before and after each oscillatory process, the values of r are restricted to r=R NS for n F C , t < < r r M H for n F 0 , and r>r L for F ν . To calculate the number flux at a detector, for example, much higher values of r have to be considered, and it is necessary to study vacuum oscillations in more detail. Such calculations will be presented elsewhere.
From Figure 7 one can observe that the dominance of electronic neutrinos and antineutrinos found at their creation at the bottom of the accretion zone is promptly erased by kinematic decoherence in such a way that the content of the neutrinos and antineutrinos entering the MSW resonant region is dominated by nonelectronic flavors. After the adiabatic transitions provoked by MSW transitions, electronic neutrinos and antineutrinos dominate again the emission spectrum except for nonelectronic antineutrinos in the normal hierarchy. Although no energy spectrum distortion is expected, the flavor content of neutrinos and antineutrinos produced near the NS surface escape to the outer space in completely different spectra when compared with the ones in which they were created, as shown in the last row of Figure 7.

Concluding Remarks
We can now proceed to draw the conclusions and some astrophysical consequences of this work: 1. The main neutrino production channel in XRFs and BdHNe in the hypercritical accretion process is pair annihilation: nn  -+ē e . This mechanism produces an initial equal number of neutrinos and antineutrinos and an initial 7/3 relative fraction between electronic and other flavors. These features lead to a different neutrino phenomenology with respect to the typical core-collapse SN neutrinos produced via the URCA process. 2. The neutrino density is higher than both the electron density and the vacuum oscillation frequencies for the inner layers of the accretion zone, and the self-interaction potential dictates the flavor evolution along this region, as illustrated by Figure 3. This particular system leads to very fast pair conversions n n n n « m t m tē e , , induced by bipolar oscillations with oscillation length as small as O (0.05-1) km. However, due to the characteristics of the main neutrino production process, neutrinos and antineutrinos have very similar fluxes inside the neutrino emission zone and kinematic decoherence dominates the evolution of the polarization vectors. 3. The kinematic decoherence induces a fast flux equipartition among the different flavors that then enters the matter-dominated regions in which MSW resonances take place. 4. Therefore, the neutrino flavor content emerging from the Bondi-Hoyle surface to the outer space is different from the original one at the bottom of the accretion zone. As shown in Table 3, The initial 70% and 30% distribution of electronic and nonelectronic neutrinos becomes 55% and 45% or 62% and 38% for normal or inverted hierarchy, respectively. Since the n n «¯oscillations are negligible (Pontecorvo 1957(Pontecorvo , 1968Xing 2013), the total neutrino-to-antineutrino ratio is kept constant.
We have shown that such a rich neutrino phenomenology is uniquely present in the hypercritical accretion process in XRFs and BdHNe. This deserves the appropriate attention since it paves the way for a new arena of neutrino astrophysics besides SN neutrinos. There are a number of issues that still have to be investigated: 1. We have made some assumptions that, albeit being a first approximation to a more detailed picture, have allowed us to set the main framework to analyze the neutrino oscillation phenomenology in these systems. We have shown in Becerra et al. (2015) that the SN ejecta carry enough angular momentum to form a disk-like structure around the NS before being accreted. However, the knowledge of the specific properties of such possible disk-like structure surrounding the NS is still pending more accurate numerical simulations at such distance scales. For instance, it is not clear yet whether such a structure could be modeled via thin-disk or thick-disk models. We have adopted a simplified model assuming isotropic accretion and the structure of the NS accretion region used in Becerra et al. (2016), which accounts for the general physical properties of the system. In order to solve the hydrodynamics equations, the neutrino emission region features, and the neutrino flavor oscillation equations, we have assumed spherically symmetric accretion onto a nonrotating NS, a quasi-steady-state evolution parameterized by the mass accretion rate, a polytropic equation of state, and subsonic velocities inside the shock radius. The matter is described by a perfect gas made of ions, electrons, positrons, and radiation with electrons and positrons obeying a Fermi-Dirac distribution. The electron fraction was fixed and equal to 0.5. We considered pair annihilation, that such an accretion process onto the BH can explain the observed 0.1-100GeV emission in BdHNe (Ruffini et al. 2015a(Ruffini et al. , 2015b(Ruffini et al. , 2016a(Ruffini et al. , 2016bAimuratov et al. 2017; see also Aimuratov et al. in preparation). The interaction of such an ultrarelativistic expanding emitter with the interstellar medium could be a possible source of highenergy (e.g., TeV-PeV) neutrinos, following a mechanism similar to the one introduced in the traditional collapsar-fireball model of long GRBs (see, e.g., Kumar & Zhang 2015;Agostini et al. 2017, and references therein). 6. Although the symmetry between the neutrino and antineutrino number densities has allowed us to generalize the results obtained within the single-angle and monochromatic spectrum approximations, to successfully answer the question of detectability, full-scale numerical solutions will be considered in the future to obtain a precise picture of the neutrino emission spectrum. In particular, it would be possible to obtain an r-dependent neutrino spectrum without the restrictions discussed in Section 5. 7. For low accretion rates ( ´-- M M 5 10 s 5 1 ) the matter and self-interaction potentials in Equation (39) decrease and the general picture described in Figure 3 changes. The resonance region could be located around closer to the NS surface, anticipating the MSW condition λ r ∼ω r and interfering with the kinematic decoherence. This changes the neutrino flavor evolution and, of course, the emission spectrum. Hence, the signature neutrino emission spectrum associated with the least luminous XRFs might be different from the ones reported here.