Neutrino oscillations in a neutrino-dominated accretion disk around a Kerr BH

In the binary-driven hypernova (BdHN) model of long gamma-ray bursts, a carbon-oxygen star explodes as a supernova (SN) in presence of a neutron star binary companion in close orbit. Hypercritical (i.e. highly super-Eddington) accretion of the ejecta matter onto the neutron star sets in, making it reach the critical mass with consequent formation of a Kerr black hole (BH). We have recently shown that, during the accretion process onto the neutron star, fast neutrino flavour oscillations occur. Numerical simulations of the above system show that a part of the ejecta keeps bound to the newborn Kerr BH, leading to a new process of hypercritical accretion. We here address, also for this phase of the BdHN, the occurrence of neutrino flavour oscillations given the extreme conditions of high density (up to $10^{12}$ g cm$^{-3}$) and temperatures (up to tens of MeV) inside this disk. We estimate the evolution of the electronic and non-electronic neutrino content within the two-flavour formalism ($\nu_{e}\nu_{x}$) under the action of neutrino collective effects by neutrino self-interactions. We find that neutrino oscillations inside the disk have frequencies between $\sim (10^{5}$-$10^{9})$ s$^{-1}$, leading the disk to achieve flavour equipartition. This implies that the energy deposition rate by neutrino annihilation ($\nu + \bar{\nu} \to e^{-} + e^{+}$) in the vicinity of the Kerr BH, is smaller than previous estimates in the literature not accounting by flavour oscillations inside the disk. The exact value of the reduction factor depends on the $\nu_{e}$ and $\nu_{x}$ optical depths but it can be as high as $\sim 5$. The results of this work are a first step toward the analysis of neutrino oscillations in a novel astrophysical context and, as such, deserve further attention.


I. INTRODUCTION
Neutrino flavour oscillations are now an experimental fact [1] and, in recent years, its study based only on Mikheyev-Smirnov-Wolfenstein (MSW) effects [2,3] has been transformed by the insight that refractive effects of neutrinos on themselves due to the neutrino self-interaction potential are essential. Their behaviour in vacuum, matter or by neutrino self-interactions have been studied in the context of early universe evolution [4][5][6][7][8][9][10][11][12][13][14][15], solar and atmospheric neutrino anomalies [16][17][18][19][20][21][22][23][24], and core-collapse supernovae (SN)  and references therein. We are here interested in astrophysical situations when neutrino self-interactions becomes more relevant than the matter potential. This implies systems in which a high density of neutrinos is present and in fact most of the literature on neutrino self-interaction dominance are concentrated on supernova neutrinos. It has been there shown how collective effects, such as synchronized and bipolar oscillations, change the flavor content of the emitted neutrinos when compared with the original content deep inside the exploding star.
This article aims to explore the problem of neutrino flavour oscillations in the case of long gamma-ray bursts (GRBs) within the binary-driven hypernova (BdHN) scenario. The GRB progenitor is a binary system composed of a carbon-oxygen star (CO core ) and a companion neutron star (NS) [52][53][54][55][56][57]. The CO core explodes as SN ejecting matter that produces a hypercritical accretion (i.e. highly super-Eddington) process onto the NS companion. The NS reaches the critical mass for gravitational collapse, hence forming a rotating black hole (BH). The emission of neutrinos is a crucial ingredient since they act as the main cooling process that allows the accretion onto the NS to proceed at very high rates of up to 1 M s −1 [56,[58][59][60][61].
In [62], we studied the neutrino flavour oscillations in the aforementioned hypercritical accretion process onto the NS, all the way to BH formation. We showed that, the density of neutrinos on top the NS, in the accreting "atmosphere", is such that neutrino self-interactions dominate the flavour evolution leading to collective effects. The latter induce in this system quick flavour conversions with a short oscillation length as small as (0.05-1) km. Far from the NS surface the neutrino density decrease and so the matter potential and MSW resonances dominate the flavour oscillations. The main result has been that the neutrino flavour content emerging on top of the accretion zone was completely different compared to the one created at the bottom of it.
In the BdHN scenario, part of the SN ejecta keeps bound to the newborn Kerr BH, forming an accretion disk onto it. In this context, the study of accretion disks and their nuances related to neutrinos is of paramount importance to shed light on this aspect of the GRB central engine. In most cases, the mass that is exchanged in close binaries has enough angular momentum so that it cannot fall radially. As a consequence, the gas will start rotating around the star or BH forming a disk. However, the magneto-hydrodynamics that describe the behaviour of accretion disks are too complex to be solved analytically and full numerical analysis are time-consuming and costly. To bypass this difficulty, different models make approximations that allow casting the physics of an accretion disk as a two-or even one-dimensional problem. These approximations can be can be pigeonholed into four categories: symmetry, temporal evolution, viscosity and dynamics. Almost all analytic models are axially symmetric. This is a sensible assumption for any physical systems that rotates. Similarly, most models are timeindependent although this is a more complicated matter. A disk can evolve in time in several ways. For example, the accretion rateṀ depends on the external source of material which need not be constant and, at the same time, the infalling material increases the mass and angular momentum of the central object, constantly changing the gravitational potential. Additionally, strong winds and outflows can continually change the mass of the disk. Nonetheless,Ṁ (x, t) =Ṁ = constant is assumed. Viscosity is another problematic approximation. For the gas to spiral down, its angular momentum needs to be reduced by shear stresses. These come from the turbulence driven by differential rotation and the electromagnetic properties of the disk [63][64][65][66] but, again, to avoid magneto-hydrodynamical calculations, the turbulence accounted for using a phenomenological viscosity α = constant, such that the kinematical viscosity takes the form ν ≈ αHc s , where c s is the local isothermal sound speed of the gas and H is the height of the disk measured from the plane of rotation (or halfthickness). This idea was first put forward by [67] and even though there is disagreement about the value and behaviour of the viscosity constant, and it has been criticized as inadequate [68][69][70][71], several thriving models use this prescription. Finally, the assumptions concerning the dynamics of the disk are related to what terms are dominant in the energy conservation equation and the Navier-Stokes equation that describe the fluid (apart from the ones related to symmetry and time independence). In particular, it amounts to deciding what cooling mechanisms are important, what external potentials should be considered and what are the characteristics of the internal forces in the fluid. The specific tuning of these terms breeds one of the known models: thin disks, slim disks, advection-dominated accretion flows (ADAFs), thick disks, neutrino-dominated accretion flows (NDAFs), convection-dominated accretion flows (CDAFs), luminous hot accretion flows (LHAFs), advection-dominated inflow-outflow solutions (ADIOS) and magnetized tori. The options are numerous and each model is full of subtleties making accretion flows around a given object an extremely rich area of research. For useful reviews and important articles with a wide range of subjects related to accretion disks see [72][73][74][75][76][77][78][79][80][81][82][83][84][85][86] and references therein.
NDAFs are of special interest for GRBs. They are hyperaccreting slim disks, optically thick to radiation that can reach high densities ρ ≈ 10 10 -10 13 g cm −3 and high temperatures T ≈ 10 10 -10 11 K around the inner edge. Under these conditions, the main cooling mechanism is neutrino emission since copious amounts of (mainly electron) neutrinos and anti-neutrinos are created by electron-positron pair annihilation, URCA and nucleon-nucleon bremsstrahlung processes, and later emitted from the disk surface. These νν pairs might then annihilate above the disk producing an e − e + dominated outflow. NDAFs were proposed as a feasible central engine for GRBs in [87] and have been studied extensively since [88][89][90][91][92][93][94][95][96][97][98][99]. In [90] and later in [94], it was found that the inner regions of the disk can be optically thick to ν eνe trapping them inside the disk, hinting that NDAFs may be unable to power GRBs. Yet, the system involves neutrinos propagating through dense media and, consequently, an analysis of neutrino oscillations, missing in the above literature, must be performed. Fig. 1 represents the standard situation of the physical system of interest. The dominance of the self-interaction potential induces collective effects or decoherence. In either case, the neutrino flavour content of the disk changes. Some recent articles are starting to recognize their role in accretion disks and spherical accretion [62,[100][101][102][103]. The energy deposition rate above and accretion disk by neutrino-pair annihilation as a powering mechanism of GRBs in NDAFs can be affected by neutrino oscillation in two ways. The neutrino spectrum emitted at the disk surface depends not only on the disk temperature and density but also the neutrino flavour transformations inside the disk. Also, once the neutrinos are emitted they undergo flavour transformations before being annihilated. Our main objective is to analyze the consequences of neutrino flavour oscillations in NDAFs.
The article is organized as follows. In Sec. II we outline the features of NDAFs making emphasis on the assumptions needed to derive the equations. In Sec. III we discuss the general details of the equation that drive the evolution of neutrino oscillations and use the information of the previous section to build a simple model that adds this dynamic to NDAFs. In Sec. IV we give some details on the initial conditions needed to solve the equations of accretion disks and neutrino oscillations. In Sec. V we discuss the main results of our calculations and analyze in detail the neutrino oscillation phenomenology in accretion disks. Finally, we present in Sec. VI the conclusions of this work. Additional technical details are presented in a series of appendices at the end.

A. Units, velocities and Averaging
Throughout this article we will use Planck units c = G = = k B = k e = 1. To describe the spacetime around a Kerr BH of mass M we use the metric g µν in Boyer- Figure 1. Schematic representation of the physical system. Due to conditions of high temperature and density, neutrinos are produced in copious amounts inside the disk. Since they have a very low cross-section, neutrinos are free to escape but not before experiencing collective effects due to the several oscillation potentials. The energy deposition rate of the process ν +ν → e − + e + depends on the local distribution of electronic and non-electronic (anti)-neutrinos which is affected by the flavour oscillation dynamics.
Lindquist coordinates, with spacelike signature, and with a dimensionless spin parameter a = J/M 2 so that the line element is in coordinates (t, r, θ, φ). The covariant components (g) µν of the metric are and its determinant is with the well known functions Σ = r 2 + M 2 a 2 cos 2 θ and ∆ = r 2 − 2M r + M 2 a 2 . We denote the coordinate frame by CF. Note that these coordinates can be used by an observer on an asymptotic rest frame. The angular velocity of the locally non-rotating frame (LNRF) is and in Eq. (2) it can be seen explicitly that if an observer has a an angular velocity dφ/dt = ω it would not measure any differences between the ±φ directions. The LNRF is defined by orthonormality and the coordinate change φ LNRF =φ = φ − ω t [104,105]. We assume that the disk lies on the equatorial plane of the BH (θ = π/2) . This way we represent the average movement of the fluid by geodesic circular orbits with angular velocity Ω = dφ/dt = u φ /u t plus a radial velocity so that the local rest frame (LRF) of the fluid is obtained by performing, first, an azimuthal Lorentz boost with velocity βφ to a corotating frame (CRF) [106], and then a radial Lorentz boost with velocity βr. Clearly, the metric on the LNRF, CRF and LRF is diag(−1, 1, 1, 1). The expression for the angular velocity of circular orbits is obtained by settinġ r =r = 0 in the r-component of the geodesic equation where (+) is for prograde orbits and (−) is for retrograde orbits. We will limit our calculations to prograde movement with 0 ≤ a ≤ 1 but extension to retrograde orbits is straightforward. Finally, we can get the components of the 4-velocity of the fluid by transforming u LRF = (1, 0, 0, 0) back to the CF leaving βr to be determined by the conservation laws. In Eq. (6) we have replaced βφ with Eq. (A3). A discussion on the explicit form of the transformations and some miscellaneous results are given in Appendix A. We will also assume that the disk is in a steady-state. This statement requires some analysis. There are two main ways in which it can be false: I. As matter falls into the BH, its values M and a change [107,108], effectively changing the spacetime around it. For the spacetime to remain the same we require Ω −1 t acc = ∆M 0 /Ṁ acc , where ∆M 0 is the total mass of the disk andṀ acc is the accretion rate. The characteristic accretion time must be bigger than the dynamical time of the disk so that flow changes due to flow dynamics are more important than flow changes due to spacetime changes. Equivalent versions of this condition that appear throughout disk accretion articles are t dym t visc and where it is understood that the accretion rate obeyṡ M acc ≈ ∆M 0 /t acc . To put this numbers into perspective, consider a solar mass BH (M = 1M ) and a disk with mass between ∆M 0 = (1 − 10)M . For accretion rates up toṀ acc = 1M /s the characteristic accretion time is t acc (1 − 10) s, while Ω −1 ∼ (10 −5 − 10 −1 ) s between r = r ISCO and r = 2000M . Consequently, a wide range of astrophysical system satisfy this condition and it is equivalent to claiming that both ∂ t and ∂ φ are Killing fields.
II. At any point inside the disk, any field Ψ(t, r, θ, φ) that reports a property of the gas may variate in time due to the turbulent motion of the flow. So, to assume that any field is time-independent and smooth enough in r for its flow to be described by Eq. (6) means replacing such field by its average over an appropriate spacetime volume. The same process allows to choose a natural set of variables that split the hydrodynamics into r-component equations and θ-component equations. The averaging process has been explained in [106,109,110]. We include the analysis here and try to explain it in a self-consistent manner. The turbulent motion is characterized by the ed-dies. The azimuthal extension of the largest eddies can be 2π, like waves crashing around an island, but their linear measure cannot be larger than the thickness of the disk, and, as measured by an observer on the CRF, their velocity is of the order of βr so that their period along the r component is ∆t ≈ (Thickness)/βr [e.g. §33 111]. If we denote by H the average half-thickness of the disk as measured by this observer at r over the time ∆t, then the appropriate volume V is composed by the points (t, r, θ, φ) such that t ∈ [t * − ∆t/2, t * + ∆t/2], θ ∈ [θ min , θ max ] and φ ∈ [0, 2π), where we have transformed ∆t and ∆r back to the CF using Eqs. (A4) as approximations. The values θ min and θ max correspond to the upper and lower faces of the disk, respectively. Then, the average takes the form The steady-state condition is achieved by requiring that the Lie derivative of the averaged quantity along the Killing field ∂ t vanishes: L ∂t ψ = 0. Note that the thickness measurement performed by the observer already has an error ∼ M 2 a 2 H 3 /6r 4 since it extends the Lorentz frame beyond the local neighbourhood but, if we assume that the disk is thin (H/r 1), and we do, this error remains small. At the same time, we can take all metric components evaluated at the equator and use Eq. (6) as the representative average velocity. Under these conditions, we have θ max − θ min ≈ 2H/r and the term −g/g rr in Eq. (8) cancels out. It becomes clear that an extra θ integral is what separates the radial and polar variables. In other words, the r-component variables are the vertically integrated fields The vertical equations of motion can be obtained by setting up Newtonian (with relativistic corrections) equations for the field ψ (r, θ) at each value of r [see e.g. 86,109,112,113].

B. Conservation Laws
The equations of evolution of the fluid are contained in the conservation laws ∇ µ T µν = 0 and ∇ µ (ρu µ ) = 0. The most general stress-energy tensor for a Navier-Stokes viscous fluid with heat transfer is [114,115] T = Ideal Fluid (10) where ρ, P , U , ζ, η, q, P and σ are the rest-mass energy density, pressure, internal energy density, dynamic viscosity, bulk viscosity, heat-flux 4-vector, projection tensor and shear tensor, respectively, and thermodynamic quantities are measured on the LRF. We do not consider electromagnetic contributions and ignore the causality problems associated with the equations derived from this stress-energy tensor since we are not interested in phenomena close to the horizon [106]. Before deriving the equations of motion and to add a simple model of neutrino oscillations to the dynamics of disk accretion we must make some extra assumption. We will assume that the θ integral in Eq. (8) can be approximated by for any field ψ. Also, we use Stokes' hypothesis (ζ = 0). Since we are treating the disk as a thin differentially rotating fluid, we will assume that, on average, the only non-zero component of the shearing stress on the CRF is σrφ (there are torques only on the φ direction), and qθ is the only non-zero component of the energy flux (on average the flux is vertical). By u µ σ µν = 0 and Eq. (A7) we have Finally, the turbulent viscosity is estimated to be ∼ l∆u where l is the size of the turbulent eddies and ∆u is the average velocity difference between points in the disk separated by a distance l. By the same arguments in [ §33 111] and in Sec. II B, l can be at most equal to 2H and ∆u can be at most equal to the isothermal sound speed c s = ∂P/∂ρ or else the flow would develop shocks [76]. The particular form of c s can be calculated from Eq. (16). This way we get with α ≤ 1 and Π = ρ + U + P . In a nutshell, this is the popular α-prescription put forward by [67]. As we mentioned at the end of Sec. II A, on the CRF for a fixed value of r, the polar equation takes the form of Euler's equation for a fluid at rest where the acceleration is given by the tidal gravitational acceleration, that is, the θ component of the fluid's path-lines relative acceleration in the θ direction.
with R the Riemann curvature tensor. With uμ ≈ (1, 0, 0, 0), Eq. (11), Eq. (A8) and assuming that there is no significant compression of the fluid under the action of the tidal force, integration of this equation yields the relation up to second order in π/2 − θ where we used the condition P = 0 at the disk's surface. Hence, the average pressure inside the disk is [cf. 86, 94, 113] The equation of mass conservation is obtained by directly inserting into Eq. (B5) the averaged density and integrating vertically where the term 2Hrρu r is identified as the average inward mass flux through a cylindrical surface of radius r per unit azimuthal angle and thus must be equal to the accretion rate divided by 2π. The same process applied to Eq. (B4) yields the energy conservation equation where factors proportional to H/r were ignored and we assume Π ≈ ρ to integrate the second term on the left hand side. is the average energy density measured on the LRF (see the discussion around Eq. (B7)). The first term on the right hand side is the viscous heating rate F heat and the second term is the cooling rate F cool . The last constitutive equation is obtained by replacing the density in Eq. (17) using Eq. (B12)

C. Equations of State
We consider that the main contribution to the restmass energy density of the disk is made up of neutrons, protons and ions. This way ρ = ρ B = n B m B with baryon number density n B and baryon mass m B equal to the atomic unit mass. The disk's baryonic mass obeys Maxwell-Boltzmann statistics and its precise composition is determined by the Nuclear Statistical Equilibrium (NSE). We denote the mass fraction of an ion i by X i = ρ i /ρ B (if i = p or n then we are referring to proton or neutrons) and it can be calculated by the Saha equation [116,117] With constraints In these equations T , are the temperature, atomic number, neutron number, proton number, electron fraction (electron abundance per baryon), ion abundance per baryon, nuclear partition function, chemical potential (including the nuclear rest-mass energy) and ion binding energy. The µ C i are the Coulomb corrections for the NSE state in a dense plasma (see Appendix C). The binding energy data for a large collection of nuclei can be found in [118] and the temperature-dependent partition functions are found in [119,120]. Even though we take into account Coulomb corrections in NSE we assume that the baryonic mass can be described by an ideal gas 1 and 1 Since bulk viscosity effects appear as a consequence of correlations between ion velocities due to Coulomb interactions and of large relaxation times to reach local equilibrium, the NSE and ideal gas assumptions imply that imposing Stokes' hypothesis becomes de rigueur [115,121,122] The disk also contains photons, electrons, positrons, neutrinos and anti-neutrinos. As it is usual in neutrino oscillations analysis, we distinguish only between electron (anti)-neutrinos ν e , (ν e ) and x (anti)-neutrinos ν x (ν x ), where x = µ + τ is the superposition of muon neutrinos and tau neutrinos. Photons obey the usual relations while, for electrons and positrons we have with ξ = T /m e and written in terms of the generalized Fermi functions In these equations η e ± = (µ e ± − m e ) /T is the electron (positron) degeneracy parameter without rest-mass contributions (not to be confused with η in Sec. (II B)). Since electrons and positrons are in equilibrium with photons due to the pair creation and annihilation processes (e − +e + → 2γ) we know that their chemical potentials are related by µ e + = −µ e − , which implies η e + = −η e − − 2/ξ. From the charge neutrality condition and we obtain For neutrinos, the story is more complicated. In the absence of oscillations and if the disk is hot and dense enough for neutrinos to be trapped within it and in thermal equilibrium, n ν , U ν , P ν can be calculated with Fermi-Dirac statistics using the same temperature T n trapped where it is understood that F(η) = F(y = 0, η) with η ν(ν) = µ ν(ν) /T and the ultra-relativistic approximation m ν 1 for any neutrino flavour is used. If thermal equilibrium is has not been achieved, Eq. (26) cannot be used. Nevertheless, at any point in the disk and for a given value of T and ρ, (anti)-neutrinos are being created through several processes. The processes we take into account are pair annihilation e − + e + → ν +ν, electron or positron capture by nucleons p + e − → n + ν e or n + e + → p +ν e , electron capture by ions A + e − → A + ν e , plasmon decayγ → ν +ν and nucleon-nucleon bremsstrahlung n 1 + n 2 → n 3 + n 4 + ν +ν. The emission rates can be found in Appendix D. The chemical equilibrium for these processes determines the values of η ν(ν) . In particular, satisfy all equations. Here, Q = (m n − m p )/m e ≈ 2.531. Once the (anti)-neutrino number and energy emission rates (R i , Q i ) are calculated for each process i, the (anti)neutrino thermodynamic quantities are given by Remember we are using Planck units so in these expressions there should be an H/c instead of just an H. The transition for each (anti)-neutrino flavour between both regimes occurs when Eq. (26b) and Eq. (28b) are equal and it can be simulated by defining the parameter With this equation, the (anti)-neutrino average energy can be defined as and the approximated number and energy density are Note that both Eq. (28c) and (31c) are approximations since they are derived from equilibrium distributions, but they help make the transition smooth. Besides, the neutrino pressure before thermal equilibrium is negligible.
This method was presented in [94] where it was used only for electron (anti)-neutrinos. The total (anti)-neutrino number and energy flux through one the disk's faces can be approximated witḣ Here, τ νi is the total optical depth for the (anti)neutrino ν i (ν i ). Collecting all the expressions we write the total internal energy and total pressure The (anti)-neutrino energy flux through the disk faces contributes to the cooling term in the energy conservation equation but it is not the only one. Another important energy sink is photodisintegration of ions. To calculate it we proceed as follows. The energy spent to knocking off a nucleon of an ion i is equal to the binding energy per nucleon B i /A i . Now, consider a fluid element of volume V whose moving walls are attached to the fluid so that no baryons flow in or out. The total energy of photodisintegration contained within this volume is the sum over i of (energy per nucleon of ion i)×(# of freed nucleons of ion i inside V ). This can be written as If we approximate B i /A i by the average binding energy per nucleonB (which is a good approximation save for a couple of light ions) the expression becomes n B VB i X f,i = n B VBX f = n B VB(X p + X n ). The rate of change of this energy as measured by an observer on the LRF with proper time λ is The derivative of n B V vanishes by baryon conservation. Transforming back to CF and taking the average we find the energy density per unit time used in disintegration of ions The average energy density measured on the LRF appearing in Eq. (18) is  Table I. Mixing and squared mass differences as they appear in [123]. Error values in parenthesis are shown in 3σ interval. The squared mass difference is defined as ∆m 2 = m 2 3 − m 2 2 + m 2 1 /2 and its sign depends on the hierarchy m1 < m2 < m3 or m3 < m1 < m2.
Finally, a similar argument allows us to obtain the equation of lepton number conservation.
For any lepton , the total lepton number density is ∈{e,µ,τ } (n − n¯ + n ν − nν ). So, with Eq. (25), calculating the rate of change as before, using Gauss' theorem and taking the average we get where the right hand side represents the flux of lepton number through the disk's surface.

III. NEUTRINO OSCILLATIONS
To study the flavour evolution of neutrinos within a particular system, a Hamiltonian governing neutrino oscillation must be set up. The relative strength of the potentials appearing in such Hamiltonian depends on four elements: geometry, mass content, neutrino content and neutrino mass hierarchy. Geometry refers to the nature of net neutrino fluxes and possible gravitational effects. Mass and neutrino content refers to the distribution of leptons of each flavour (e, µ, τ ) present in the medium. Finally, mass hierarchy refers to the relative values of the masses m 1 , m 2 , m 3 for each neutrino mass eigenstates (see Table I). We dedicate this section to a detailed derivation of the equations of flavour evolution for a neutrino dominated accretion disk. To maintain consistency with traditional literature of neutrino oscillations we will reuse some symbols appearing in previews sections. To avoid confusion we point out that the symbols in this section are independent of previews sections unless we explicitly draw a comparison.

A. Equations of Oscillation
The equations that govern the evolution of an ensemble of mixed neutrinos are the Boltzmann collision equations The collision terms should include the vacuum oscillation plus all possible scattering interactions that neutrinos undergo through their propagation. For free streaming neutrinos, only the vacuum term and the forwardscattering interactions are taken into account so that the equations become Here, H p,t (H p,t ) is the oscillation Hamiltonian for (anti)-neutrinos and ρ p,t (ρ p,t ) is the matrix of occupation numbers: (ρ p,t ) ij = a † j a i p,t for neutrinos and ((ρ p,t ) ij = ā † iā j p,t for anti-neutrinos), for each momentum p and flavours i, j. The diagonal elements are the distribution functions f νi(νi) (p) such that their integration over the momentum space gives the neutrino number density n νi of a determined flavour i at time t. The offdiagonal elements provide information about the overlapping between the two neutrino flavours. Taking into account the current-current nature of the weak interaction in the standard model, the Hamiltonian for each equation is [124][125][126] where G F is the Fermi coupling constant, Ω p,t is the matrix of vacuum oscillation frequencies, l p,t andl p,t are matrices of occupation numbers for charged leptons built in a similar way to the neutrino matrices, and v p,t = p/p is the velocity of a particle with momentum p (either neutrino or charged lepton).
As stated before, we will only consider two neutrino flavours: e and x = µ + τ . Three-flavour oscillations can be approximated by two-flavour oscillations as a result of the strong hierarchy of the squared mass differences |∆m 2 13 | ≈ |∆m 2 23 | |∆m 2 12 |. 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 ν e ν x andν e ν x , and each term in the Hamiltonian governing oscillations becomes a 2 × 2 Hermitian matrix. Now, consider an observer on the LRF (which is almost identical to the CRF due to Eq. (7) at a point r. In its spatial local frame, the unit vectorsx,ŷ,ẑ are parallel to the unit vectorsr,θ,φ of the CF, respectively. Solving Eq. (39) in this coordinate system would yield matrices ρ,ρ as functions of time t. However, in our specific physical system, both the matter density and the neutrino density vary with the radial distance from the BH. This means that the equations of oscillations must be written in a way that makes explicit the spatial dependence, i.e. in terms of the coordinates x, y, z. For a collimated ray of neutrinos, the expression dt = dr would be good enough, but for radiating extended sources or neutrino gases the situation is more complicated.
In Eq. (39) we must replace the matrices of occupation numbers by the space-dependent Wigner functions ρ p,x,t (andρ p,x,t ) and the total time derivative by the Liouville operator [127,128] ρ p,x,t = Explicit Time In this context, x represents a vector in the LRF. In the most general case, finding ρ p,x,t andρ p,x,t means solving a 7D neutrino transport problem in the variables x, y, z, p x , p y , p z , t. Since our objective is to construct a simple model of neutrino oscillations inside the disk, to obtain the specific form of Eq. (39) we must simplify the equations by imposing on it conditions that are consistent with the assumptions made in Sec. II.
• Due to axial symmetry, the neutrino density is constant along the z direction. Moreover, since neutrinos follow geodesics, we can setṗ z ≈ṗ φ = 0.
• Within the thin disk approximation (as represented by Eq. (11)) the neutrino and matter densities are constant along the y direction and the momentum change due to curvature along this direction can be neglected, that is,ṗ y ≈ 0.
• In the LRF, the normalized radial momentum of a neutrino can be written as p Hence, the typical scale of the change of momentum with radius is ∆r px,eff = d ln px , which obeys ∆r px,eff > r s for r > 2r in . This means that we can assumeṗ x ≈ 0 up to regions very close to the inner edge of the disk.
• We define an effective distance ∆r ρ,eff = . For all the systems we evaluated we found that is comparable to the height of the disk (∆r ρ,eff ∼ 2 − 5)r s . This means that at any point of the disk we can calculate neutrino oscillations in a small regions assuming that both the electron density and neutrino densities are constant.
• We neglect energy and momentum transport between different regions of the disk by neutrinos that are recaptured by the disk due to curvature. This assumption is reasonable except for regions very close to the BH but is consistent with the thin disk model [see e.g . 110]. We also assume that the neutrino content of neighbouring regions of the disk (different values of r) do not affect each other. As a consequence of the results discussed above, we assume that at any point inside the disk and at any instant of time an observer can describe both the charged leptons and neutrinos as isotropic gases around small enough regions of the disk. This assumption is considerably restrictive but we will generalize it in Sec. V.
The purpose of these approximations is twofold. (1) We can reduce the problem considerably since they allow us to add the neutrino oscillations dynamics of the disk by simply studying the behaviour of neutrinos at each point of the disk using the constant values of density and temperature at that point. As we will see in Sec. V, this corresponds to a transient state of the disk since, very fast, neighbouring regions of the disk start interacting. (2) Also, the approximations allow us to simplify the equations of oscillation considering that all but the first term in Eq. (41) vanish, leaving only a time derivative. In addition, both terms of the form v q,t · v p,t in Eq. (40) average to zero so that ρ p,x,t = ρ p,t andρ p,x,t =ρ p,t .
We are now in a position to derive the simplified equations of oscillation for this particular model. Let us first present the relevant equations for neutrinos. Due to the similarity between H p,t andH p,t , the corresponding equations for anti-neutrinos can be obtained analogously. For simplicity, we will drop the suffix t since the time dependence is now obvious. In the two-flavour approximation, ρ p is a 2 × 2 Hermitian matrix and can be expanded in terms of the Pauli matrices σ i and a polarization vector P p = (P x , P y , P z ) in the neutrino flavour space, such that where f p = Tr[ρ p ] = f νe (p) + f νx (p) is the sum of the distribution functions for ν e and ν x . Note that the z component of the polarization vector obeys Hence, this component tracks the fractional flavour composition of the system. Appropriately normalizing ρ p allows to define a survival and mixing probability The Hamiltonian can be written as a sum of three interaction terms: The first term is the Hamiltonian in vacuum [27]: where ω p = ∆m 2 /2p, B = (sin 2θ, 0, − cos 2θ) and θ is the smallest neutrino mixing angle in vacuum. The other two terms in Eqs. (40) are special since they make the evolution equations non-linear. Since we are considering that the electrons inside the form an isotropic gas, the vector v q in the first integral is distributed uniformly on the unit sphere and the factor v q · v p averages to zero. After integrating the matter Hamiltonian is given by where λ = √ 2G F (n e − − n e + ) is the charged current matter potential and L = (0, 0, 1). Similarly, the same product disappears in the last term and after integrating we get Clearly, P = P p dp/(2π) 3 . Introducing every Hamiltonian term in Eqs. (39), and using the commutation relations of the Pauli matrices, we find the equations of oscillation for neutrinos and anti-neutrinos for each momentum mode ṗ where we have assumed that the total neutrino distribution remains constantḟ p = 0. In this form, it is clear how the polarization vectors can be normalized. Performing the transformation P p /f p → P p andP p /f p →P p and, multiplying and dividing the last term by the total neutrino density Eqs. (49) can be written aṡ This is the traditional form of the equations in terms of the vacuum, matter and self-interaction potentials ω p , λ and µ with Different normalization schemes are possible [see e.g. 36,49,126,129]. By assuming that we can solve the equations of oscillation with constant potentials λ and µ simplifies the problem even further. Following [130], with the vector transformation (a rotation around the z axis of flavour space) eliminating the λ potential but making B time dependent. Defining the vector S p = P p +P p and, adding and subtracting Eq. (53a) and Eq. (53b) we geṫ The last approximation is true if we assume that the self-interaction potential is larger than the vacuum potential ω p /µ 1. We will show later that this is the case for thin disks (see Fig. 5). The first equation implies that all the vectors S p and their integral S evolve in the same way, suggesting the relation S p = f p +f p S. By replacing in Eq. (54b) and integratinġ where ω = ω p f p +f p dp/(2π) 3 is the average vacuum oscillation potential. The fact that in our model the equations of oscillations can be written in this way has an important consequence. Usually, as it is done in supernovae neutrino oscillations, to solve Eq. (50) we would need the neutrino distributions throughout the disk. If neutrinos are trapped, their distribution is given by Eq. (26). If neutrinos are free, their temperature is not the same as the disk's temperature. Nonetheless, we can approximate the neutrino distribution in this regime by a Fermi-Dirac distribution with the same chemical potential as defined by Eq. (27) but with an effective temperature T eff ν . This temperature can be obtained by solving the equation where ζ (3) is Apéry's constant (ζ is the Riemann zeta function) and Li s (z) is Jonquière's function. For convenience and considering the range of values that the degeneracy parameter reaches (see Sec. VI), we approximate the effective temperature of electron neutrinos and anti-neutrinos with the expressions with constants a = 0.0024, b = −0.085, c = 0.97. However, Eq. (55) allow us to consider just one momentum mode, and the rest of the spectrum behaves in the same way.

IV. INITIAL CONDITIONS AND INTEGRATION
In the absence of oscillations, we can use Eqs. (18), (16) and (37) to solve for the set of functions η e − (r), ξ (r), Y e (r) using as input parameters the accretion ratė M , the dimensionless spin parameter a, the viscosity parameter α and the BH mass M . From [86,94] we learn that neutrino dominated disks require accretion between 0.01 M s −1 and 1 M s −1 (this accretion rate range vary depending on the value of α). For accretion rates smaller than the above lower value the neutrino cooling is not efficient and, for rates larger than the upper value, the neutrinos are trapped within the flow. We also limit ourselves to the above accretion rate range since it is consistent with the one expected to occur in a BdHN (see e.g. [56,60,61]).
We also know that high spin parameter, high accretion rate and low viscosity parameter produces disks with higher density and higher temperature. This can be explained using the fact that several variables of the disk, like pressure, density and height are proportional to a positive power of the quotientṀ /α. To avoid this semidegeneracy in the system, reduce the parameter space and considering that we want to study the dynamics of neutrino oscillations inside the disk, we fix the BH mass at M = 3M , the viscosity parameter at α = 0.01 and the spin parameter at a = 0.95 while changing the accretion rate.
Eqs. (18) and (37) are first order ordinary differential equations and since we perform the integration from an external (far away) radius r out up to the innermost stable circular orbit r in we must provide two boundary conditions at r out . Following the induced gravitational collapse (IGC) paradigm of gamma-ray bursts (GRBs) associated with type Ib/c supernovae we assume that at the external edge of the disk, the infalling matter is composed mainly by the ions present in the material ejected from an explosion of a carbon-oxygen core, that is, mainly oxygen and electrons. This fixes the electron fraction Y e (r out ) = 0.5. The second boundary condition can be obtained by the relation (T η + m B ) √ g tt = constant [131][132][133], with η the degeneracy parameter of the fluid. If we require the potentials to vanishes at infinity and invoking Euler's theorem we arrive at the relation in the weak field limit For a classical gas composed of ions and electrons this relation becomes Eq. (59) can be used with Eq. (16) and Eq. (33) to solve for η e − (r out ), ξ (r out ). The value of r out is chosen to be at most the circularization radius of the accreting material as described in [59,134]. We can estimate this radius by solving for r in the expression of the angular momentum per unit mass for a equatorial circular orbits. So using Eq. (6) we need to solve where x = r/M which yields r out ∼ 1800r s and the expression is in geometric units. Finally, for the initial conditions to be accepted, they are evaluated by the gravitational instability condition [135] Rθtθt θ=π/2 Integration of the equations proceeds as follows, with the initial conditions we solve Eq. (37) to obtain the electron fraction in the next integration point. With the new value of the electron fraction we solve the differentialalgebraic system of Eqs. (18) and (16) at this new point. This process continues until the innermost stable circular orbit r in is reached.
To add the dynamics of neutrino oscillations we proceed same as before but at each point of integration, once the values of Y e , η and ξ are found, we solve Eq. (50) for the average momentum mode to obtain the survival probabilities as a function of time. We then calculate the new neutrino and anti-neutrino distributions with the conservation of total number density and the relations Since the disk is assumed to be in a steady-state, we then perform a time average of Eq. (62) as discussed in Sec. II. With the new distributions, we can calculate the new neutrino and anti-neutrino average energies and use them to re-integrate the disk equations.
Neutrino emission within neutrino-cooled disks is dominated by electron and positron capture which only produces electron (anti)-neutrinos. The second most important process is electron-positron annihilation but it is several orders of magnitude smaller. In Fig. 2 we show the total number emissivity for these two processes for an accretion rate ofṀ = 0.1M s −1 . Other cases behave similarly. Moreover, although the degeneracy parameter suppresses the positron density, a high degeneracy limit does not occur in the disk and the degeneracy is kept low at values between ∼ (0.2-3), as shown in Fig. 3. The reason for this is the effect of high degeneracy on neutrino cooling. Higher degeneracy leads to a lower density of positrons which suppresses the neutrino production and emission, which in turn leads to a lower cooling rate, higher temperature, lower degeneracy and higher positron density. This equilibrium leads, via the lepton number conservation Eq. (37), to a balance between electronic and non-electronic neutrino densities within the inner regions of the disk. Given this fact, to solve the equations of oscillations, we can approximate the initial conditions of the polarization vectors with P =P ≈ (0, 0, 1).
V. RESULTS AND ANALYSIS In Fig. 3 and Fig. 4 we present the main features of accretion disks for the parameters M = 3M , α = 0.01, a = 0.95, and three different accretion ratesṀ = 1M s −1 ,Ṁ = 0.1M s −1 andṀ = 0.01M s −1 . It exhibits the usual characteristics of thin accretions disks. High accretion rate disks have higher density, temperature and electron degeneracy. Also, for high accretion rates, the cooling due to photodisintegration and neutrino emission kicks in at larger radii. For all cases, as the disk heats up, the number of free nucleons starts to increase enabling the photodisintegration cooling at r ∼ (100-300)r s . Only the disintegration of alpha particles is important and the nucleon content of the infalling matter is of little consequence for the dynamics of the disk. When the disk reaches temperatures ∼ 1.3 MeV, the electron capture switches on, the neutrino emission becomes significant and the physics of the disk is dictated by the energy equilibrium between F heat and F ν . The radius at which neutrino cooling becomes significant (called ignition radius r ing ) is defined by the condition F ν ∼ F heat /2. For the low accretion rateṀ = 0.01M s −1 , the photodisintegration cooling finishes before the neutrino cooling becomes significant, this leads to fast heating of the disk. Then the increase in temperature triggers a strong neutrino emission that carries away the excess heat generating a sharp spike in F ν surpassing F heat by a factor of ∼ 3.5. This behaviour is also present in the systems studied in [94], but there it appears for fixed accretion rates and high viscosity (α = 0.1). This demonstrates the semidegeneracy mentioned in Sec. V. The evolution of the fluid can be tracked accurately through the degeneracy parameter. At the outer radius, η e − starts to decrease as the temperature of the fluid rises. Once neutrino cooling becomes significant, it starts to increase until the disk reaches the local balance between heating and cooling. At this point, η e − stops rising and is maintained (approximately) at a constant value. Very close to r in , the zero torque condition of the disk becomes important and the viscous heating is reduced drastically. This is reflected in a sharp decrease in the fluid's temperature and increase in the degeneracy parameter. For the high accretion rate and additional effect has to be taken into account. Due to high ν e optical depth, neutrino cooling is less efficient, leading to an increase in temperature and a second dip in the degeneracy parameter. This dip is not observed in low accretion rates because τ νe does not reach significant values. With the information in Fig. 3 we can obtain the oscillation potentials which we plot in Fig. 5. Since the physics of the disk for r < r ign is independent of the initial conditions at the external radius and for r > r ign the neutrino emission is negligible, the impact of neutrino oscillations is important only inside r ign . We can see that the discussion at the end of Sec. III A is justified since for r in < r < r ign the potentials obey the relation Generally, the full dynamics of neutrino oscillations is a rather complex interplay between the three potentials, yet it is possible to understand the neutrino response in the disk using some numerical and algebraic results obtained in [33,36,126] and references therein. Specifically, we know that if µ ω , as long as the MSW condition λ ω is not met (precisely our case), collective effects should dominate the neutrino evolution even if λ µ. On the other hand, if µ ω , the neutrino evolution is driven by the relative values between the matter and vacuum potentials (not our case). With Eq. (55) we can build a very useful analogy. These equations are analogous to the equations of motion of a simple mechanical pendulum with a vector position given by S, precessing around with angular momentum D, subjected to a gravitational force ω µB with mass µ −1 . Using Eq. (63) obtains the expression |S| = S ≈ 2+O( ω /µ). Calculating ∂ t (S · S) it can be checked that this value is conserved up to fluctuations of order ω /µ. The analogous angular momentum is D = P −P = 0. 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 ϕ between S and the z-axis such that S = S (sin ϕ, 0, cos ϕ) .
Note that the only non-zero component of D is ycomponent. From Eq. (55) we finḋ The above equations can be equivalently written as where we have introduced the inverse characteristic time k by which is related to the anharmonic oscillations of the pendulum. The role of the matter potential λ is to logarithmically extend the oscillation length by the relation [126] τ = −k −1 ln k θ (k 2 + λ 2 ) The total oscillation time can then be approximated by the period of an harmonic pendulum plus the logarithmic extension The initial conditions of Eq. (63) imply so that ϕ is a small angle. The potential energy for a  simple pendulum is If k 2 > 0, which is true for the normal hierarchy ∆m 2 > 0, we expect small oscillations around the initial position since the system begins in a stable position of the potential. The magnitude of flavour conversions is of the order ∼ ω /Sµ 1. We stress that normal hierarchy does not mean an absence of oscillations but rather imperceptible oscillations in P z . No strong flavour oscillations are expected. On the contrary, for the inverted hierarchy ∆m 2 < 0, k 2 < 0 and the initial ϕ indicates that the system begins in an unstable position and we expect very large anharmonic oscillations. P z (as well asP z ) oscillates between two different maxima passing through a minimum −P z (−P z ) several times. This implies total flavour conversion: all electronic neutrinos (anti-neutrinos) are converted into non-electronic neutrinos (anti-neutrinos) and vice-versa. This has been called bipolar oscillations in the literature [44]. If the initial condition are not symmetric as in Eq. (63), the asymmetry is measured by a constant ς =P z /P z ifP z < P z or ς = P z /P z ifP z > P z so that 0 < ς < 1. Bipolar oscillations are present in an asymmetric system as long as the relation is obeyed [126]. If this condition is not met, instead of bipolar oscillation we get synchronised oscillations.
Since we are considering constant potentials, synchronised oscillations is equivalent to the normal hierarchy case. From Fig. 5 we can conclude that in the normal hierarchy case, neutrino oscillations have no effects on neutrino-cooled disks under the assumptions we have made. On the other hand, in the inverted hierarchy case, we expect extremely fast flavour conversions with periods of order t osc ∼ (10 −9 − 10 −5 ) s for high accretion rates and t osc ∼ (10 −8 − 10 −5 ) s for low accretion rates, between the respective r in and r ign . For the purpose of illustration we solve the equations of oscillations for theṀ = 0.1M s −1 case at r = 10r s . The electronic (anti)-neutrino survival probability at this point is shown in Fig. 6 for inverted hierarchy and normal hierarchy, respectively. On both plots, there is no difference between the neutrino and anti-neutrino survival probabilities. This should be expected since for this values of r the matter and self-interaction potentials are much larger than the vacuum potential, and there is virtually no difference between Eq. (50a) and Eq. (50b). Also, as mentioned before, note that the (anti)-neutrino flavour proportions remain virtually unchanged for normal hierarchy while the neutrino flavour proportions change drastically. The characteristic oscillation time of the survival probability in inverted hierarchy found on the plot is which agree with the ones given by Eq. (70) up to a factor of order one. Such a small value suggests extremely quick ν eνe → ν xνx oscillations. A similar ef-   57), the (anti)-neutrino spectrum for both flavours can be constructed. But, more importantly, this means that the local observer at that point in the disk measures, on average, an electron (anti)-neutrino loss of around 8% which is represented by an excess of non-electronic (anti)-neutrinos. In Sec. III A we proposed to calculate neutrino oscillations assuming that small neighbouring regions of the disk are independent and that neutrinos can be viewed as isotropic gases in those regions. However, this cannot be considered a steady-state of the disk. To see this consider Fig. 4. The maximum value of the neutrino optical depth is of the order of 10 3 for the highest accretion rate, meaning that the time that takes neutrinos to travel a distance of one Schwarzschild inside the disk radius obeys t rs Max (τ ν ) r s ≈ 10 −2 s. (75) which is lower than the accretion time of the disk as discussed in Sec. II but higher than the oscillation time. Different sections of the disk are not independent since they, very quickly, share (anti)-neutrinos created with a non-vanishing momentum along the radial direction. Furthermore, the oscillation pattern between neighbouring regions of the disk is not identical. In Fig. 7 we show the survival probability as a function of time for different (but close) values of r forṀ = 0.1M s −1 . The superposition between neutrinos with different oscillation histories has several consequences: (1) It breaks the isotropy of the gas because close to the BH, neutrinos are more energetic and their density is higher producing a radially directed net flux, meaning that the factor v q,t · v p,t does not average to zero. This implies that realistic equations of oscillations include a multi-angle term and a radially decaying neutrino flux similar to the situation in SN neutrinos. (2) It constantly changes the neutrino content ant any value of r independently of the neutrino collective evolution given by the values of the oscillation potentials at that point. This picture plus the asymmetry that electron and non-electron neutrinos experience through the matter environment (electron (anti)-neutrinos can interact through n + ν e → p + e − and p +ν e → n + e + ), suggests that the disk achieves complete flavour equipartition (decoherence). We can identify two competing causes, namely, quantum decoherence and kinematic decoherence.
Quantum decoherence is the product of collisions among the neutrinos or with a thermal background medium can be understood as follows [136]. From Appx. D 2 we know that different (anti)-neutrino flavours posses different cross-sections and scattering rates Γ νi,νi . In particular, we have Γ νx ≈ Γν x < Γν e < Γ νe . An initial electron (anti)-neutrino created at a point r will begin to oscillate into ν x (ν x ). The probability of finding it in one of the two flavors evolves as previously discussed. However, in each interaction n+ν e → p+e − , the electron neutrino component of the superposition is absorbed, while the ν x component remains unaffected. Thus, after the interaction the two flavors can no longer interfere. This allows the remaining ν x oscillate and develop a new coherent ν e component which is made incoherent in the next interaction. The process will come into equilibrium only when there are equal numbers of electronic and nonelectronic neutrinos. That is, the continuous emission and absorption of electronic (anti)-neutrinos generates a non-electronic (anti)-neutrinos with an average probability of P νe→νe in each interaction and once the densities of flavours are equal, the oscillation dynamic stops. An initial system composed of ν e ,ν e turns into an equal mixture of ν e ,ν e and ν x ,ν x , reflected as an exponential damping of oscillations. For the particular case in which non-electronic neutrinos can be considered as sterile (do not interact with the medium), the relaxation time of this process can be approximated as [137,138] t Q = 1 2l νν ω 2 sin 2 2θ + 2l νν λ 2 ω 2 sin 2 2θ (76) where l νν represents the (anti)-neutrino mean free path. Kinematic decoherence is the result of a non-vanishing flux term such that at any point, (anti)-neutrinos travelling in different directions, do not experience the same self-interaction potential due to the multi-angle term in the integral of Eq. (40). Different trajectories do not oscillate in the same way, leading to a de-phasing and a decay of the average P ν→ν and thus to the equipartition of the overall flavour content. The phenomenon is similar to an ensemble of spins in an inhomogeneous magnetic field. In [35] it is shown that for asymmetric νν gas, even an infinitesimal anisotropy triggers an exponential evolution towards equipartition, and in [36] it was shown that if the symmetry between neutrinos and anti-neutrinos is not broken beyond the limit of 25%, kinematic decoherence is still the main effect of neutrino oscillations. As a direct consequence of the νν symmetry present within the ignition radius of accretion disks (see Fig. 3), equipartition among different neutrino flavours is expected. This multi-angle term keeps the order of the characteristic time t osc of Eq. (70), unchanged and kinematic decoherence happens within a few oscillation cycles. The oscillation time gets smaller closer to the BH due to the 1/µ 1/2 dependence. Therefore, we expect that neutrinos emitted within the ignition radius will be equally distributed among both flavours in about few microseconds. Once the neutrinos reach this maximally mixed state, no further changes are expected. We emphasize that kinematic decoherence does not mean quantum decoherence. Figs. 6 and Fig. 7 clearly show the typical oscillation pattern which happens only if quantum coherence is still acting on the neutrino system. Kinematics decoherence, differently to quantum decoherence, is just the result of averaging over the neutrino intensities resulting from quick flavour conversion. Therefore, neutrinos are yet able to quantum oscillate if appropriate conditions are satisfied.
Simple inspection of Eq. (70) and Eq. (76) with Fig. 4 yields t osc t Q . Clearly the equipartition time is dominated by kinematic decoherence. These two effects are independent of the neutrino mass hierarchy and neutrino flavour equipartition is achieved for both hierarchies. Within the disk dynamic, this is equivalent to imposing the condition P νe→νe = Pν e→νe = 0.5. Fig. 8 shows a comparison between disks with and without neutrino flavour equipartition for the three accretion rates considered. The role of equipartition is to increase the disk's density, reduce the temperature and electron fraction, and further stabilize the electron degeneracy for regions inside the ignition radius. The effect is mild for low accretion rates and very pronounced for high accretion rates. This result is in agreement with our understanding of the dynamics of the disk and can be explained in the following way. In low accretion systems the neutrino optical depth for all flavors is τ νν 1 and the differences between the cooling fluxes, as given by Eq. (32) are small. Hence, when the initial (mainly electron flavour) is redistributed among both flavours, the total neutrino cooling remains virtually unchanged and the disk evolves as if equipartition had never occurred save the new emission flavour content. On the other hand, when accretion rates are high, the optical depth obeys τ νx ≈ τν x τν e < τ ν ∼ 10 3 . The ν e cooling is heavily suppressed while the others are less so. When flavours are redistributed, the new ν x particles a free to escape, enhancing the total cooling and reducing the temperature. As the temperature decreases, so do the electron and positron densities leading to a lower electron fraction. The net impact of flavour equipartition is to make the disk evolution less sensitive to ν e opacity and, thus, increase the total cooling efficiency. As a consequence, once the fluid reaches a balance between F + and F ν , this state is kept without being affected by high optical depths and η e − stays at a constant value until the fluid reaches the zero torque condition close to r in . Note that for every case, inside the ignition radius, we have τν e ≈ τ νx = τν x so that equipartition enhances, mainly, neutrino cooling F ν (and not anti-neutrino cooling Fν).
The quotient between neutrino cooling with and without equipartition can be estimated with This relation exhibits the right limits. From Fig. 3 we see that E νe ≈ E νx . Hence, If 1 τ νe > τ νx , then F eq ν = F ν and the equipartition is unnoticeable. But if 1 < τ νx < τ νe then F eq ν /F ν > 1. In our simulations, this fraction reaches values of 1.9 forṀ = 1M s −1 to 2.5 forṀ = 0.01M s −1 .
The disk variables at each point do not change beyond a factor of order 5 in the most obvious case. However, these changes can be important for cumulative quantities, e.g. the total neutrino luminosity and the total energy deposition rate into electron-positron pairs due to neutrino anti-neutrino annihilation. To see this we perform a Newtonian calculation of these luminosities following [86,87,99,[139][140][141][142], and references therein. The neutrino luminosity is calculated by integrating the neutrino cooling flux throughout both faces of the disk  Table II. Comparison of total neutrino luminosities Lν and annihilation luminosities Lνν between disk with and without flavour equipartition.
The factor 0 < C cap < 1 is a function of the radius (called capture function in [108]) that accounts for the proportion of neutrinos that are re-captured by the BH and, thus, do not contribute to the total luminosity. For a BH with M = 3M and a = 0.95, the numerical value of the capture function as a function of the dimensionless distance x = r/r s is well fitted by with a relative error smaller than 0.02%. To calcu-late the energy deposition rate, the disk is modeled as a grid of cells in the equatorial plane. Each cell k has a specific value of differential neutrino luminosity ∆ k νi = F k νi r k ∆r k ∆φ k and average neutrino energy E νi k . If a neutrino of flavour i is emitted from the cell k and an anti-neutrino is emitted from the cell k , and, before interacting at a point r above the disk, each travels a distance r k and r k , then, their contribution to the energy deposition rate at r is (see Appx. D 3 for details) The total neutrino annihilation luminosity is simply the sum over all pairs of cells integrated in space where A is the entire space above (or below) the disk. In Table II we show the neutrino luminosities and the neutrino annihilation luminosities for disks with and without neutrino collective effects. In each case, flavour equipartition induces a loss in L νe by a factor of ∼3, and a loss in Lν e luminosity by a factor of ∼2. At the same time, L νx and Lν e are increased by a factor ∼10. This translates into a reduction of the energy deposition rate due to electron neutrino annihilation by a factor of ∼7 while the energy deposition rate due to non-electronic neutrinos goes from being negligible to be of the same order of the electronic energy deposition rate. The net effect is to reduce the total energy deposition rate of neutrino annihilation by a factor of ∼ (3−5) for the accretion rates considered. In particular we obtain a for a factor 3.03 and 3.66 forṀ = 1 M s −1 andṀ = 0.01 M s −1 , respectively and a factor 4.73 forṀ = 0.1 M s −1 . The highest value correspond to and intermediate value of the accretion rate because, for this case, there is a ν e cooling suppression (τ νe > 1) and the quotient τ νe /τ νx is maximal. By Eq. (77), the difference between the respec-tive cooling terms is also maximal. In Fig. 9 we show the energy deposition rate per unit volume around the BH for each flavour with accretion ratesṀ = 1 M s −1 anḋ M = 0.1 M s −1 . There we can see the drastic enhancement of the non-electronic neutrino energy deposition rate and the reduction of the electronic deposition rate. Due to the double peak in the neutrino density forṀ = 0.01 M s −1 case (see Fig. 3), the deposition rate per unit volume also shows two peaks. One at r s < r < 2r s and the other at 10 r s < r < 11 r s . Even so, the behaviour is similar to the other cases.

VI. CONCLUDING REMARKS
The generation of a seed, energetic e − e + plasma seems to be a general prerequisite of GRB theoretical models for the explanation of the prompt gamma-ray emission. The e − e + pair annihilation produce photons leading to an opaque pair-photon plasma that self-accelerates, expanding to ultrarelativistic Lorentz factors of the order of 10 2 -10 3 (see, e.g., [143][144][145]). The reaching of transparency of MeV-photons at large Lorentz factor and corresponding large radii is requested to solve the so-called compactness problem posed by the observed non-thermal spectrum in the prompt emission [146][147][148]. There is a vast literature on this subject and we refer the reader to  Figure 9. Comparison of the neutrino annihilation luminosity per unit volume ∆Qν iνi = k,k ∆Q ν iνi kk between disk without (left column) and with (right column) flavour equipartition. [149][150][151][152][153][154], and references therein, for further details. Neutrino-cooled accretion disks onto rotating BHs have been proposed as a possible way of producing the above-mentioned e − e + plasma. The reason is that such disks emit a large amount of neutrino and antineutrinos that can undergo pair annihilation near the BH [87][88][89][90][91][92][93][94][95][96][97][98][99]. The viability of this scenario clearly depends on the energy deposition rate of neutrino-antineutrinos into e − e + and so on the local (anti)-neutrino density and energy.
We have here shown that, inside these hyperaccreting disks, a rich neutrino oscillations phenomenology is present due to the high neutrino density. Consequently, the neutrino/antineutrino emission and the corresponding pair annihilation process around the BH leading to electron-positron pairs, are affected by neutrino flavour conversion. Using the thin disk and α-viscosity approximations, we have built a simple stationary model of general relativistic neutrino-cooled accretion disks around a Kerr BH, that takes into account not only a wide range of neutrino emission processes and nucleosynthesis but also the dynamics of flavour oscillations. The main assumption relies on considering the neutrino oscillation behaviour within small neighbouring regions of the disk as independent from each other. This, albeit being a first approximation to a more detailed picture, has allowed us to set the main framework to analyze the neutrino oscillations phenomenology in inside neutrino-cooled disks.
In the absence of oscillations, a variety of neutrinocooled accretion disks onto Kerr BHs, without neutrino flavour oscillations, have been modelled in the literature (see e.g. [87,94,99,106] and [86] for a recent review). The physical setting of our disk model follows closely the ones considered in [94], but with some extensions and differences in some aspects: (i) The equation of vertical hydrostatic equilibrium, Eq. (16), can be derived in several ways [106,109,113]. We followed a particular approach consistent with the assumptions in [109], in which we took the vertical average of a hydrostatic Euler equation in polar coordinates. The result is a an equation that leads to smaller values of the disk pressure when compared with other models. It is expected that the pressure at the centre of the disk is smaller than the average density multiplied by the local tidal acceleration at the equatorial plane. Still, the choice between the assortment of pressure relations is tantamount to a fine-tuning of the model. Within the thin disk approximation, all these approaches are equivalent since they all assume vertical equilibrium and neglect self-gravity.
(ii) Following the BdHN scenario for the explanation of GRBs associated with Type Ic SNe (see Sec. II), we considered a gas composed of 16 O at the outermost radius of the disk and followed the evolution of the ion content using the Saha equation to fix the local NSE. In [94], only 4 He is present and, in [99], ions up to 56 Fe are introduced. The affinity between these cases implies that this particular model of disk accretion is insensible to the initial mass fraction distribution. This is explained by the fact that the average binding energy for most ions is very similar, hence any cooling or heating due to a redistribution of nucleons, given by the NSE, is negligible when compared to the energy consumed by direct photodisintegration of alpha particles. Additionally, once most ions are dissociated, the main cooling mechanism is neutrino emission that is similar for all models, modulo the supplementary neutrino emission processes included in addition to electron and positron capture. However, during our numerical calculations, we noticed that the inclusion of non-electron neutrino emission processes can reduce the electron fraction by up to ∼ 8%. This effect is observed again during the simulation of flavour equipartition alluding to the need for detailed calculations of neutrino emissivities when establishing NSE state. We obtain similar results to [94] (see Fig. 3), but by varying the accretion rate and fixing the viscosity parameter. This suggests that a more natural differentiating set of variables in the hydrodynamic equations of an α-viscosity disk is the combination of the quotientṀ /α and eitheṙ M or α. This result is already evident in, for example, Fig. 11 and Fig. 12 of [94], but was not mentioned there.
Concerning neutrino oscillations, we showed that the conditions inside the ignition radius, the oscillation potentials follow the relation ω µ λ, as it is illustrated by Fig. 5. We also showed that the within this region the number densities of electron neutrinos and anti-neutrinos are very similar. As a consequence of this particular environment very fast pair conversions ν eνe ν xνx , induced by bipolar oscillations, are obtained for the inverted mass hierarchy case with oscillation frequencies between 10 9 s −1 and 10 5 s −1 . For the normal hierarchy case no flavour changes are observed (see Fig. 6 and Fig. 7). Bearing in mind the magnitude of these frequencies and the low neutrino travel times through the disk, we conclude that an accretion disk under our main assumption cannot represent a steady state. However, using numerical and algebraic results obtained in [33,35,36], and references therein, we were able to generalize our model to a more realistic picture of neutrino oscillations. The main consequence of the interaction between neighbouring regions of the disk is the onset of kinematic decoherence in a timescale of the order of the oscillation times. Kinematic decoherence induces fast flavour equipartition among electronic and non-electronic neutrinos throughout the disk. Therefore, the neutrino content emerging from the disk is very different from the one that is usually assumed (see e.g. [100,155]). The comparison between disks with and without flavour equipartition is summarized in Fig. 8 and Table II. We found that flavour equipartition, while leaving anti-neutrino cooling practically unchanged, it enhances neutrino cooling by allowing the energy contained (and partially trapped inside the disk due to high opacity) within the ν e gas to escape in the form of ν x , rendering the disk insensible to the electron neutrino opacity. We give in Eq. (77) a relation to estimate the change in F ν as a function of τ νe τ νx that describes correctly the behaviour of the disk under flavour equipartition. The variation of the flavour content in the emission flux implies a loss in L νe and an increase in L νx and Lν e . As a consequence, the total energy deposition rate of the process ν +ν → e − + e + is reduced. We showed that this reduction can be as high 80% and is maximal whenever the quotient τ νe /τ νx is also maximal and the condition τ νe > 1 is obtained.
At this point we can identify several issues which have still to be investigated: (1) Throughout the accretion disk literature, several fits to calculate the neutrino and neutrino annihilation luminosity can be found (see e.g. [86] and references therein). However, all these fits were calculated without taking into account neutrino oscillations. Since we have shown that oscillations directly impact luminosity, these results need to be extended.
(2) Additionally, the calculations of the neutrino and neutrino annihilation luminosities we have performed ignore general relativistic effects and the possible neutrino oscillations from the disk surface to the annihilation point. In [156], it has been shown that general relativistic effects can enhance the neutrino annihilation luminosity in a neutron star binary merger by a factor of 10. In [87], however, it is argued that in BHs this effect has to be mild since the energy gained by falling into the gravitational potential is lost by the electron-positron pairs when they climb back up. Nonetheless, this argument ignores the bending of neutrino trajectories and neutrino capture by te BH which can be significant for r 10r s .
In [157], the increment is calculated to be no more than a factor of 2 and can be less depending on the geometry of the emitting surface. But, as before, they assume a purely ν eνe emission and ignore oscillations after the emission. Simultaneously, the literature on neutrino oscillation above accretion disks (see e.g. [100]) do not take int account oscillations inside the disk and assume only ν eνe emission. A similar situation occurs in works studying the effect of neutrino emission on r-process nucleosynthesis in hot outflows (wind) ejected from the disk (see e.g. [158]). It is still unclear how the complete picture (oscillations inside the disk → oscillations above the disk + relativistic effects) affect the final energy deposition. We are currently working on the numerical calculation of the annihilation energy deposition rate using a ray tracing code and including neutrino oscillations both inside the disk and after their emission from the disk surface. These results will be the subject of a future publication.
This articles serves as a primer that has allowed us to identify key theoretical and numerical features involved in the study of neutrino oscillations in neutrino-cooled accretion disks. The unique conditions inside the disk and its geometry lend themselves to varied neutrino oscillations that can have an impact on a wide range of astrophysical phenomena: from e − e + plasma production above BHs in GRB models, to r-process nucleosynthesis in disk winds and possible MeV neutrino detectability. As such, this topic deserves appropriate attention since it paves the way for new astrophysical scenarios for testing neutrino physics.

Appendix A: Transformations and Christoffel symbols
For the sake of completeness, here we give the explicitly the transformation used in Eq. (6) and the Christoffel symbols used during calculations. The coordinate transformation matrices between the CF and the LNRF on the tangent vector space is [105] so that the basis vectors transform as ∂ν = e µν ∂ µ , that is, with e T . For clarity, coordinates on the LNRF have a caret (xμ), coordinates on the CRF have a tilde (xμ) and coordinates on the LRF have two (xμ). An observer on the LNRF sees the fluid elements move with an azimuthal velocity βφ. This observer then can perform a Lorentz boost L βφ to a new frame. On this new frame an observer sees the fluid elements falling radially with velocity βr, so it can perform another Lorentz boost L βr to the LRF. Finally, the transformation between the the LRF and the CF coordinates x µ = e μ ρ (L βφ )ρ α (L βr )α ν xν = A μ ν xν, where the components of A and A −1 are Since Lorentz transformations do not commute, the transformation A raises the question: what happens if we invert the order? In this case, we would not consider a corotating frame but a cofalling frame on which observers see fluid elements, not falling, but rotating. The new transformation velocities β r , β φ are subject to the conditions β φ = γ r βφ, β r = βr/γφ and γ r γ φ = γrγφ. Although both approaches are valid, considering that the radial velocity is an unknown, the first approach is clearly cleaner. To obtain the coordinate transformation between the CF and the CRF A μ ν and Aν µ we can simply set βr = 0 in Eqs. (A2). With this we can calculate The non-vanishing Christoffel symbols are Using the connection coefficients and the metric, both evaluated at the equatorial plane we can collect several equations for averaged quantities. The expansion of the fluid world lines is There are several ways to obtain an approximate version of the shear tensor [e.g. 106,159,160] but by far the simplest one is proposed by [109]. On the CRF the fluid four-velocity can be approximated by uμ = (1, 0, 0, 0) by Eq. (7). Both the fluid four-acceleration a ν = u µ ∇ µ u ν and expansion parameter, Eq. (A6), vanish so that the shear tensor reduces to 2σμν = ∇μuν + ∇νuμ. In particular, the r-φ component is where cα µν are the commutation coefficients for the CRF. Finally, of particular interest is theθ component of the Riemann curvature tensor which gives a measurement of the relative acceleration in theθ direction of nearly equatorial geodesics.
Finally, we reproduce the zero torque at the innermost stable circular orbit condition that appears in [110]. Using the killing vector fields ∂ φ , ∂ t we can calculate u φ + 4rHησ r φ = 2Hu φ after integrating vertically and using Eq. (17) Analogously for ∂ t , ∂ r Ṁ 2π u t − 4rHΩησ r φ = 2Hu t after integrating vertically and using Eq. (12).
The vertical integration of the divergence of the heat flux is as follows: Since, on average, q = q θ ∂ θ , we have ∇ µ q µ = ∂ θ q θ and by Eq. (A2), q θ = rqθ. Vertically integrating yields where qθ is the averaged energy flux radiating out of a face of the disk, as measured by an observer on the LRF, which we approximate as the half-thickness of the disk H times the average energy density per unit proper time lost by the disk. With the variable change z = 8πrHησ r φ /Ṁ and y = 4πH /Ṁ the equations reduce to Using the relation ∂ r u t = −Ω∂ r u φ [see Eq. (10.7.29) in 161] and ∂ r (u t + Ωu φ ) = u φ ∂ r Ω we can combine the previous equations to obtain with A = y/∂ r Ω and B = u t + Ωu φ . To integrate these equations we use the zero torque condition z(r = r * ) = 0 where r * is the radius of the innermost stable circular orbit, which gives the relation or, equivalently, Using Eq. (6), the approximation γr ≈ 1 and the variable change r = xM 2 the integral can be easily evaluated by partial fractions where x 1 , x 2 , x 3 are the roots of the polynomial x 3 − 3x + 2a.   [162].

Appendix C: Nuclear Statistical Equilibrium
The results in this section appear in [162]. We include them here since they are necessary to solve Eq. (20). Neutrino dominated accretion disks reach densities above ∼ 10 7 g cm −3 and temperatures above ∼ 5 × 10 9 K. For these temperatures, forward and reverse nuclear reactions are balanced and the abundances in the plasma are determined by the condition µ i = Z i µ p + N i µ n , that is, the Nuclear Statistical Equilibrium. However, for densities above 10 6 g cm −3 , the electron screening of charged particle reactions can affect the nuclear reaction rates. For this reason, to obtain an accurate NSE state it is necessary to include Coulomb corrections to the ion chemical potential. The Coulomb correction to the i-th chemical potential is given by The ion coupling parameter is written in terms of the electron coupling parameter as Γ i = Γ e Z 5/3 i with Γ e = e 2 T 4πY e n B 3 where e is the electron charge. The parameters K i , C i are given in table (III).

Appendix D: Neutrino Interactions and cross-sections
In this appendix we include the neutrino emission rates and neutrino cross-sections used in the accretion disk model. These expressions have been covered in [163][164][165][166][167][168][169]. We also include the expression energy emission rate for νν annihilation into electron-positron pairs. Whenever possible we write the rates in terms of generalized Fermi functions since some numerical calculations were done following [170]. Before proceeding we list some useful expressions and constants in Planck units that will be used The numerical values can be found in [123].

Neutrino Emissivities
• Pair annihilation e − + e + → ν +ν This process generates neutrinos of all flavours but around 70% are electron neutrinos [see 62]. This is due to the fact that the only charged leptons in the disk are electrons and positrons, so creation of electron neutrinos occurs via either charged or neutral electroweak currents while creation of non-electronic neutrinos can only occur through neutral currents. Using the electron or positron four-momentum p = (E, p), the Dicus' cross-section for a particular flavour i is [163] σ D = G 2 F 12πE e − E e + C +,i m 4 e + 3m 2 e p e − ·p e + + 2 (p e − ·p e + ) 2 + 3C −,i m 4 e + m 2 e p e − ·p e + .