A multi-dimensional implementation of the Advanced Spectral neutrino Leakage scheme
Davide Gizzi Evan O'Connor Stephan Rosswog Albino Perego Ruben Cabezón Lorenzo Nativi


firstpage–lastpage
Fashion Retail: Forecasting Demand for New Items
    Aruna Rajan
    

================================================





We present a new, multi-dimensional 
implementation of
the Advanced Spectral Leakage (ASL) scheme with
the purpose of modelling neutrino-matter interactions in
neutron star mergers. A major challenge is the neutrino
absorption in the semi-transparent regime, which 
is responsible for driving winds from the merger remnant.
The composition of such winds is crucial in the understanding of the electromagnetic emission in the recently observed macronova following GW170817. 
Compared to the original version, 
we introduce an optical-depth-dependent flux factor to
model the average angle of neutrino propagation, and
a modulation that accounts for flux anisotropies in
non-spherical geometries. We scrutinise our approach by
first comparing the new scheme against the original one
for a spherically symmetric core-collapse supernova
snapshot, both in 1D and in 3D, and additionally
against a two-moment (M1) scheme 
as implemented in 1D into the code GR1D. 
The luminosities and
mean energies agree to a few percents in most tests.
Finally, for the case of a binary merger remnant snapshot we compare the new ASL 
scheme with the M1 scheme that is implemented in the 
Eulerian adaptive mesh refinement
code FLASH.
We find that the neutrino absorption distribution
in the semi-transparent regime is overall well reproduced.
Both approaches agree to within ≲ 15% for the
average energies and to better than ∼ 35 % in the total luminosities.






neutrinos, radiative transfer, hydrodynamics, star: neutron, stars: supernovae: general








§ INTRODUCTION


The first multi-messenger detection of a neutron star merger
<cit.> has brought major leaps forwards for many areas
of (astro)physics. For example, the 1.7s delay between the 
gravitational wave (GW) peak and the gamma-rays from an event
detected by the Fermi satellite <cit.> allowed to
constrain the deviations of the GW propagation speed from the 
speed of light to 1 part in 10^15 <cit.>. 
The detection further allowed for an independent measure 
of the Hubble parameter <cit.> as suggested by
<cit.>. The GW signal was followed by emission all 
across the electromagnetic (EM) spectrum
<cit.>.
The intensity of the EM emission 
detected in the aftermath of the event
decayed with a power-law exponent close to
-1.3  <cit.> as expected for
a distribution of freshly synthesised r-process elements
<cit.>.  Estimates of the involved ejecta
masses point to ∼ 0.02 for the early blue emission 
component and ∼ 0.04 for the later emerging red component <cit.>.
The early blue component requires lanthanide-free ejecta which, in turn,
are the r-process nucleosynthesis result of matter with Y_e ≳ 0.25 <cit.> (ejected
at velocities of ∼ 0.3c).  The later emerging, red component 
stems from matter with electron fractions below this threshold value.
Since the original neutron stars are in β-equilibrium
they contain only about 10^-4 of matter with 
Y_e > 0.25. Therefore, the
observed ∼ 2% of a solar mass in the
blue component point to a major re-processing of a large
fraction of the ejecta by weak interactions, raising 
Y_e via e^+ + n → p + ν̅_e
and ν_e + n → p + e^-. With GW170817 and its EM emission we have thus
witnessed weak interaction "in flagranti".
This underlines the paramount importance of
carefully modelling weak interactions and neutrino 
physics in a neutron star merger
for reliable predictions of their EM
signature.

Addressing the neutrino transport problem by solving the full 
multi-dimensional Boltzmann equation <cit.> is 
computationally very demanding and for most astrophysical problems
it is prohibitively expensive. Therefore, most multi-dimensional 
hydrodynamic studies, both for supernovae and compact binary mergers
resort to transport approximations
<cit.>.
Our particular focus here is on 
the Advanced Spectral Leakage (hereafter ASL)
<cit.>, that has 
recently been scrutinised against more 
expensive neutrino treatments <cit.> in
a core-collapse supernova context. 
Since supernovae are roughly spherically
symmetric, they allow for approximations that
are not admissible in a neutron star merger
context. In this paper we extend the original 
ASL scheme to multi-dimensional applications,
while keeping the general structure 
of the equations as presented in
the original papers
<cit.>. We 
examine the modified scheme in typical core-collapse supernova and
neutron star merger remnant
snapshots. We focus on the modelling of the
absorption in the semi-transparent regime, 
which is responsible for the neutrino-driven winds 
<cit.>, 
one of the possible
ejection channels related to
the observed blue EM component.
We rely on a spectral treatment of 
the neutrino-matter interactions, a 
key ingredient for capturing the
composition of the polar ejecta 
<cit.>. The ASL 
presented here allows for a 
computationally-inexpensive, 
spectral treatment of the neutrino
absorption in the semi-transparent
regime, and it is therefore suitable
for long-term binary merger simulations,
where more detailed neutrino treatments
require larger computational resources.

The paper is structured as follows: in Sec. <ref> we 
describe the ASL methodology
both in its original 1D version 
(Sec. <ref>)
and in our new 
multi-D implementation (Sec. <ref>). 
Simulation results are presented 
in Sec. <ref>. In
Sec. <ref> we start with
a one-dimensional core-collapse 
profile.
We then move to three-dimensional configurations
in Sec. <ref> where we first use the same
core-collapse supernova profile to inspect our 
multi-D implementation of the 
ASL scheme. 
Finally, we apply the ASL to 
a neutron star merger remnant. 
In all cases, we scrutinise the ASL scheme by
comparison with a two-moment (M1) scheme, 
and we neglect relativistic effects
everywhere.
In Sec. <ref> we summarise our results.



§ THE ADVANCED SPECTRAL LEAKAGE




 §.§ 1D implementation


We first summarise the most relevant 
features of the ASL scheme for spherically 
symmetric systems, as they are described in <cit.>. In their work,
the ASL scheme is explored both in 1D 
and multi-D spherically symmetric
core-collapse setups,
showing flexibility and overall agreement
with other neutrino transport models.

At the heart of the ASL approach is 
a spectral (i.e. energy-dependent) 
description of 
neutrino transport in 
which a neutrino energy spectrum
is initially set up to 
account for the energy-squared  
dependence of neutrino-matter
interactions.
As in most approximate treatments, 
we model neutrinos as three independent species:
electron neutrinos ν_e, electron
anti-neutrinos ν_e, and
a collective species for 
heavy-lepton neutrinos and 
anti-neutrinos ν_x.
For the interactions between neutrinos and 
matter we consider the production and absorption
of electron neutrinos and anti-neutrinos via charged current processes involving 
nucleons and nuclei, neutrino emission
by bremsstrahlung and pair processes,
and finally the scattering off nucleons and
nuclei. These reactions enter the 
computation of
the local optical depth τ_ν(E,x)
for an energy E at position x, 
which is a measure of the 
average number of interactions
a neutrino experiences before 
escaping to infinity
and defined as integral of the inverse
local mean free path λ_ν(E,x')
over a path γ


    τ_ν(E,x)= ∫_γ:x→+∞1/λ_ν(E,x'(s)) ds.


Two different optical depths are defined: 
the first is the total optical depth 
τ_ν,tot 
where both absorption and elastic scattering 
interactions are equally considered in the inverse mean free 
path calculation.
The second is the energy optical depth
τ_ν,en, which is related to
the mean free path over which neutrinos
can exchange energy with the fluid. 
An analytical estimate of the latter is given 
by computing the geometric mean between the total
and the absorption inverse mean free paths:


    λ_ν,en(E,x)= √(c λ_ν,tot(E,x)/∑_sχ_ν,ab,s(E,x)),


where χ_ν,ab,s(E,x) is the absorpitivity 
of the absorption process 's' and c is the speed of light.
Each optical depth defines 
a neutrino surface at 
τ_ν= 2/3, where neutrinos begin to
decouple from matter.

The net specific[To be explicit:
we always use "specific" for quantities 
on a per mass basis.],
spectral[We always use
"spectral" for quantities in units of E^2dE.]
neutrino emission rate
(units of s^-1g^-1erg^-3)
is initially calculated as a smooth 
interpolation between the 
production r_ν,prod(E,x) 
and diffusion
r_ν,diff(E,x) rate


    r_ν(E,x)= r_ν,prod(E,x)  r_ν,diff(E,x)/r_ν,prod(E,x)+r_ν,diff(E,x).


where r_ν,prod(E,x) 
depends on the production timescale t_ν,prod,
which in turn is set by the local 
emissivity, while 
r_ν,diff(E,x)
depends on the timescale over
which neutrinos diffuse out
of the system, t_ν,diff. 
This timescale
is set by the local opacity via 
∼τ_ν,tot^2(E,x).
Eq. (<ref>) 
favours r_ν,diff(E,x) 
in optically thick conditions
(τ_ν,tot(E,x) ≫ 1) and 
r_ν,prod(E,x) in
optically thin conditions 
(τ_ν,tot(E,x) ≲ 1).
We add two further corrections. 
First, when a large amount of neutrinos
is emitted at the neutrino surface or 
is locally produced, Pauli blocking
occurs as a consequence
of the fermionic nature of neutrinos.
Second, emission in optically
thin regimes provided by r_ν,prod
is assumed isotropic, and a fraction 
of neutrinos are emitted toward the optically
thick regime. 
To account for these effects,
neutrino emission is
reduced by introducing a Pauli 
blocking parameter α_ν,blk:
r_ν→ (1-α_ν,blk)r_ν.
Second, during the 
diffusion process in the optically thick 
regime, neutrinos thermalize to lower energies 
and therefore the spectrum at the neutrino
surface is softened. The softening of the spectrum 
is included via the term 1/Ψ_ν(x)exp(-τ_ν,en(E,x)/τ_cut),
with Ψ_ν(x) defined as


    Ψ_ν(x)= ∫_0^+∞r_ν(E,x)e^-τ_ν,en(E,x)/τ_cutE^2 dE/∫_0^+∞r_ν(E,x)E^2 dE,


where τ_cut parametrizes the typical number of interactions required to thermalize neutrinos.
The equation for the neutrino emission rate finally becomes
(we will occasionally refer to this
emission as cooling)

    r_ν(E,x)= (1 - α_ν,blk)r_ν(E,x)1/Ψ_ν(x)exp(-τ_ν,en(E,x)/τ_cut).

The values of 
α_ν,blk and τ_cut
are geometry-dependent, and must ideally be calibrated
every time the system geometry changes significantly over
the time of the simulation. A good trade-off is to provide
fixed values for these parameters that are able to approximately
reproduce the neutrino properties dynamically in comparison to 
other transport approaches. So far, 
α_ν,blk and 
τ_cut have been calibrated in 
the context of spherically symmetric 
core-collapse supernovae simulations 
against full Boltzmann neutrino transport
<cit.>.
Electron neutrinos and anti-neutrinos 
have α_ν,blk∼ 0.55.
Heavy-lepton neutrinos are generally
subdominant and their emission in optically
thin regime is negligible, therefore 
α_ν,blk∼ 0.
For the thermalization coefficient we adopt τ_cut=20
for all neutrino species.
Although a new calibration in the context of binary merger simulations would be preferred, for the time being we
assume the same values adopted for core-collapse
simulations, leaving the detailed 
binary merger calibration task to a future work.

The absorption of neutrinos in the optically thin
regime is hereafter referred to as heating, 
and the specific, spectral absorption rate
(units of s^-1g^-1erg^-3)
is defined as


    h_ν(E,x)= 1/ρ(x) n_ν,τ≲ 1 χ_ν,ab F_e^∓ H,


where ρ(x) is the mass density of the fluid 
at position x, χ_ν,ab the 
absorpitivity, n_ν,τ≲ 1 the 
neutrino number density in optically thin regime, H=
 exp(-τ_ν,tot) 
ensures the heating to be applied only in 
the optically thin regime.
All quantities on the RHS of Eq. (<ref>) are functions
of energy E and position
x.
The Pauli blocking factor
for final state electrons and positrons 
is given by


    F_e^∓= 1 - 1/exp((E± Q ∓μ_e)/k_BT)+1,


where k_B is the 
Boltzmann constant, Q ≈ 1.293 MeV
is the difference between 
neutron and proton rest mass energy, 
μ_e is the
electron chemical potential and T is the
fluid temperature.
The form of n_ν,
τ≲ 1 for a spherically symmetric 
heating is


    n_ν,τ≲ 1(E,R)= l_ν(E,R)/4π R^2 c μ_ν(E,R),


where l_ν(E,R) is the total, 
spectral number 
rate (in s^-1erg^-3)
at radius R obtained as solution of a differential 
equation that accounts for both 
emission and absorption of neutrinos while they 
propagate from the centre of the system to a distance R

    dl_ν(E,R)/dR= 4π R^2ρ(R)r_ν(E,R) - χ_ab(E,R)/cH(E,R)l_ν(E,R).

In Eq. (<ref>) μ_ν(E,R) 
is called flux factor.
It corresponds to the average 
of the cosine of the propagation angle for the
free streaming neutrinos. 
An analytic approximation is given by
<cit.>


    μ_ν(E,R)= 1/2(1+√(1-(R_ν(E)/max(R,R_ν(E)))^2)),


where R_ν(E) is the neutrino surface 
radius for energy E. Far from the neutrino surface 
(R ≫ R_ν) the neutrino flux points 
toward the observer direction and the propagation  
angle is 0, i.e. μ_ν(E,R)= 1.
Close to the neutrino surface (R ∼ R_ν)
and assuming isotropic neutrino emission
above the plane 
tangential to it
μ_ν(E,R) ∼ 1/2.

Given the spectral, specific rates 
r_ν(E,x) and 
h_ν(E,x) at each point from
Eqs. (<ref>) and (<ref>),
the energy-integrated emission and 
absorption specific rates are


    R_ν^k(x)   = ∫_0^+∞ r_ν(E,x) E^2+k  dE,
        
    
        H_ν^k(x)   = ∫_0^+∞ h_ν(E,x) E^2+k  dE,


respectively, where k = 0 
specifies the number rate 
(g^-1s^-1) and k = 1 
the energy rate (erg g^-1s^-1). Eqs. (<ref>) and
(<ref>) define 
the specific number and energy net rates Q_ν^k = 0(x) and  Q_ν^k = 1(x)


    Q_ν^k = 0(x)= R_ν^k = 0(x) - H_ν^k = 0(x),
    Q_ν^k = 1(x)= R_ν^k = 1(x) - H_ν^k = 1(x),


from which the total neutrino number net rate L_ν^k = 0 
and the neutrino luminosity L_ν^k = 1 
can be derived by integrating over the volume V
of the fluid


    L_ν^k = 0   = ∫_VQ_ν^k = 0(x)ρ(x) dV,
    
        L_ν^k = 1   = ∫_VQ_ν^k = 1(x)ρ(x) dV.


From the last two equations
the neutrino average energy is
calculated as


    ⟨E_ν⟩= L_ν^k = 1/L_ν^k = 0.


The root-mean squared (rms) energy 
can be defined too as


    E_rms= √(L_ν^k = 2/L_ν^k = 0).


In Sec (<ref>) we will mainly
refer to Eq. (<ref>) to describe
the neutrino energy, but we 
additionally provide values for the rms
energies for completeness.
From Eqs. (<ref>) and (<ref>) 
we can also recover the local net 
change in the total lepton number fraction Y_l(x)


    Y_l(x)= m_b (H_ν_e^k = 0(x) - H_ν̅_e^k = 0(x) + R_ν̅_e^k = 0(x) - R_ν_e^k = 0(x))


where m_b is the baryon mass, 
and the change of
the total specific 
matter internal energy u(x)


    u(x)=    
        - (R_ν_e^k = 1(x) +  R_ν̅_e^k = 1(x) +  4R_ν̅_μ,τ^k = 1(x))  + 
          
        + (H_ν_e^k = 1(x) + H_ν̅_e^k = 1(x)).


Both Y_l(x) and u(x)
contain variations in the lepton number fraction and
specific internal energy by local emission and 
absorption of neutrinos and by neutrino diffusion
that would dominate in the optically thick regime. 
Denoting the variation of the number 
and energy of trapped neutrinos per baryon
as Y_ν(x) and 
Z_ν(x) driven 
by diffusion, we can recover the change 
in the electron 
fraction Y_e(x) 


    Y_e(x)= Y_l(x) - Y_ν_e(x) + Y_ν̅_e(x)


and the rate of change of 
the specific internal energy
due to local neutrino 
emission and absorption 
e(x)


    e(x)= 
        u(x) - 1/m_b (Z_ν_e(x) 
        + Z_ν̅_e(x) + 4Z_ν_μ,τ(x)).


Y_ν(x) and 
Z_ν(x) 
are evaluated at first order in time as


    Y_ν(x)   = Y_ν,t+Δ t(x)-Y_ν(x)/Δ t,
    Z_ν(x)   = Z_ν,t+Δ t(x)-Z_ν(x)/Δ t,


where Δ t is the current time step. 
An adequate time step should ensure that the 
variation of all the  variables for which 
the ASL provide a source term in the 
hydrodynamic equations is less than a given, 
small percentage (namely ≲ 1%).
The number 
and energy of trapped neutrinos per baryon
at time t and location x are related to
the neutrino-trapped distribution function
f_ν,t^tr(E,x) by 


    Y_ν,t(x)   =
        4π/(hc)^3m_b/ρ(x)∫ f_ν,t^tr(E,x)
        E^2 dE,
    
        Z_ν,t(x)   =
        4π/(hc)^3m_b/ρ(x)∫ f_ν,t^tr(E,x)
        E^3 dE,


where h is the Planck constant.
Starting from Y_ν,t(x) and Z_ν,t(x), 
f_ν,t^tr(E,x) is first recovered
on the basis of equilibrium arguments. In particular, we assume a distribution of the form


    f_ν^tr(E,x)= 
         f_ν^eq(E,x)(1-e^ 
         -τ_ν,en(E,x)),


at time t, where f_ν^eq(E,x) is a
Fermi-Dirac distribution


    f_ν^eq(E,x)=
         1/e^(E/ ( k_ B T_ν(x) )-η_ν(x))+1,


with T_ν(x) being the neutrino
temperature, which is assumed to be equal to the matter
temperature, and 
η_ν(x) the 
degeneracy parameter, evaluated by assuming
weak equilibrium.
Second, f_ν,t^tr(E,x) is 
evolved between t and t+Δ t
considering production and diffusion of neutrinos
as two competing processes.
Namely, we integrate the equation


    d f_ν^tr/dt=
      ḟ_ν,prod^tr + 
      ḟ_ν,diff^tr,


where


    ḟ_ν,prod^tr=
    f_ν^eq - f_ν^tr/max(t_ν,prod,Δ t)exp(-t_ν,prod/t_ν,diff)


and 


    ḟ_ν,diff^tr=
    -f_ν^tr/max(t_ν,diff,Δ t)exp(-t_ν,diff/t_ν,prod).


At last, Y_ν,t + Δ t and Z_ν,t + Δ t
are recovered by using f_ν,t + Δ t^tr in Eqs. (<ref>) and (<ref>).
For more details, see <cit.>.

Note that in Eqs. (<ref>) 
and (<ref>)
we have neglected absorption by 
heavy-lepton neutrinos in 
the semi-transparent regime,
because at the decoupling surfaces 
they do not have
enough energy to produce muons and taus
by charged-current interactions.
Moreover, heating by heavy-lepton 
neutrino annihilation and 
inverse nucleon-nucleon bremsstrahlung
provides negligible contributions 
in the semi-transparent regime, 
since the opacities
are small
<cit.>.
Neglecting heating
by heavy-lepton neutrinos is therefore a
reasonable assumption.



 §.§ Multi-D implementation


All physical quantities shown 
in Sec. <ref>
are local and independent 
of the system geometry, except for Eqs. (<ref>),(<ref>) and 
(<ref>).
In particular, 
Eq. (<ref>) is straightforward 
to use only in cases where the
neutrino surface is easy to reconstruct. 
While this argument is certainly 
valid for a spherically 
symmetric configuration, for a more 
complex geometry 
like a neutron star merger remnant is not. 
Indeed, given
the presence of a torus
around the central compact object
the neutrino decoupling surface has 
larger radii on the equatorial 
plane than along the polar axis <cit.>. 
In the following, we describe
our implementation of 
Eqs. (<ref>),(<ref>)
and (<ref>) to a multi-D 
configuration.



  §.§.§ Optical depth

The computation of the 
optical depth is performed by taking the minimum among 
values of the optical depth calculated by integrating 
the neutrino mean free path over a set of predefined 
radii. In particular, given a point (x,y,z) 
we consider the following outgoing paths:

    
  * fixed (y,z), path along x
    
  * fixed (x,z), path along y
    
  * fixed (x,y), path along z
    
  * fixed x, diagonal path along y,z
    
  * fixed y, diagonal path along x,z
    
  * fixed z, diagonal path along x,y
    
  * diagonal path along x,y,z.

The likelihood of being close to
the true minimum optical depth 
at a point increases by increasing number of paths.
For this reason, the choice of diagonal paths ensures 
a more accurate calculation of the optical depth by
avoiding local
overestimates that would arise otherwise.
However, it is important to
stress that this algorithm
leads inevitably to an overestimation
of the local optical depth, 
as a consequence of 
the limited number of paths
that can be practically chosen for calculations.


  §.§.§ Flux factor

To construct a more general form for 
the flux factor that still 
resembles the general properties of 
Eq. (<ref>) we borrow the linear dependence of the inverse
flux factor from the optical depth from equation (31) 
of <cit.>. Although their equation (31) does 
not consider a spectral distribution of energies, we 
take that form to get an approximate 
expression of the flux factor 
at any energy by just extending the grey interpolation 
formula to a spectral form by adding the energy dependence.
Therefore, we use

    1/μ_ν(E,x)=  1.5 τ_ν,tot(E,x) + 1    if  τ_ν,tot(E,x) ≤ 2/3
    
    2    otherwise.

Using this expression we mimic Eq. (<ref>) 
with a flux factor μ_ν(E,x) 
tending to 1 for small optical depths 
and having its minimum value
at τ_ν,tot(E,x)= 2/3 equal to 
1/2. 
Moreover, we enforce a similar
constraint as in Eq. (<ref>) 
by setting the value of the flux 
factor to be 1/2 for any optical depth larger than 2/3.
Eq. (<ref>) 
is more suitable than Eq. (<ref>) for 
a general geometry since it only depends on
the local optical depth, and no previous 
knowledge of the radius of the neutrino surface 
is required. However, a different choice of the 
flux factor might lead to noticeable variations of 
the amount of heating. 
We will quantify the effects of this choice of flux 
factor in the next section.



  §.§.§ Flux anisotropy

Modelling the neutrino density distribution in space 
accurately requires the 
knowledge of the path along which neutrinos 
propagate from the optically thick to 
the optically thin region. While this is trivial 
in spherical symmetry, it is not so 
for a general geometry
such as a merger remnant. 
Sophisticated 
algorithms have been 
designed <cit.>, 
but come at large 
computational expense.
Here we take a simpler approach.
We modify Eq. (<ref>) 
according to the 
equation (3) of <cit.>.
In their work the total,
axially symmetric
neutrino flux an observer 
receives far from the 
neutrino emitting region 
is approximated by


    J_ν(θ,R)= 3(1+β_ν cos^2θ)/3+β_νL_ν^k = 1(θ,R)/4π R^2⟨E_ν⟩,


where θ is the angle of the observer with
respect to the source pole and β_ν= 
J_ν(θ=0,R)/J_ν(θ=π/2,R) - 1
measures the degree of flux anisotropy
which depends on the geometry of the source.
In particular, a spherically symmetric source
would have β_ν= 0. On the other hand,
in the context of mergers
the presence of a torus implies
more escaping neutrinos along 
the polar region rather 
than along the optically thick 
equatorial region, and 
therefore β_ν > 0
<cit.>.
In order to adapt Eq. (<ref>) to our
problem of reconstructing the neutrino density
distribution in our domain, we assume a 
similar modulation ∼cos^b(θ),
but we gauge the exponent b by 
comparison with an M1 neutrino
transport approach.
Moreover, we make the assumption that
the value of β_ν that would in
principle be a function of the distance to the 
source of neutrino emission 
is constant and equal to 
the value measured at large distances.
Furthermore, at distances close to the neutrino
surface the effects of the neutrino angular 
distribution arise and the hypothesis of pure
radial fluxes breaks down. We thus 
calculate the neutrino
density n_ν,τ≲ 1(E,x)
at each radius R=|x| 
as in Eq. (<ref>) by retaining 
a spherically symmetric neutrino spectral
number rate l_ν(E,|x|), but 
on one side we add the flux factor 
μ_ν(E,x) of 
Eq. (<ref>) and on the other
side we include 
the modulation factor ∼cos^b(θ)
to keep track of
the effect of the geometry of the system
on the neutrino fluxes.
Therefore we replace Eq. (<ref>) with


    n_ν,τ≲ 1(E,x)= (1+b)(1+β_νcos^b(θ))/1+b+β_νl_ν(E,|x|)/4π |x|^2cμ_ν(E,x),


where θ≡θ(x)
and we rewrite the angular
dependent pre-factor such
that the integral over the solid
angle is 4π.
Note that unlike Eq. (<ref>)
we use μ_ν(E,x) from 
Eq. (<ref>) and therefore
points in space at same distance R=|x|
to the centre of the system can 
potentially have different values of 
the flux factor.
Only for spherically symmetric systems
they are the same and 
we recover Eq. (<ref>) 
by setting β_ν=0. 
Note that in Eq. (<ref>)
we keep the assumption of axial symmetry when
calculating the flux modulation. 
A high degree of axisymmetry is 
generally expected after few tens of ms 
irrespective of the mass ratio
<cit.>.
This timescale reduces even more 
if angular momentum
transport by turbulent magnetic viscosity
is effective. Therefore, except the first few ms
after the merger, our assumption of axis symmetry
is reasonable.
As last step we estimate 
β_ν considering
an approach similar 
to the one of <cit.>. 
We divide neutrinos in two groups: 
diffusive and free streaming neutrinos. 
The former propagate outward following 
the direction of the
gradient n̂_ρ= 
-∇ρ/|∇ρ|, while the latter
are emitted isotropically. The luminosity 
at a polar angle θ coming from 
the diffusive neutrinos is selected by choosing 
those neutrinos emitted within a ring of 
width Δθ around the θ direction. 
Therefore, the luminosity per solid angle 
Λ_ν(θ)=Δ L_ν^k = 1/ΔΩ is


    Λ_ν(θ)= ∑_i R_i,ν, diff
        ^k = 1(θ)m_i/2 πsin(θ) Δθ +
        ∑_j R_j,ν,prod^k = 1m_j/4π,


 
where the index i in the sum 
is limited to those fluid points 
for which θ - Δθ/2 <
θ_i < θ + Δθ/2 
and cos(θ_i)= n̂_ρ·ê_z, ê_z being
the unit vector along the z-axis,
and the index j 
extends over the whole volume.
m_i and m_j are the masses of the fluid 
elements i and j.
The diffusive and free 
streaming contributions to
the luminosity
are calculated as


    R_i,ν, diff^k = 1(θ)   = ∫_0^+∞ 
    f_i,ν, diff(E)r_i,ν(E)E^3 dE 
    
    R_j,ν, prod^k = 1   = ∫_0^+∞ 
    f_j,ν, prod(E)r_j,ν(E)E^3 dE


with f_i,ν,diff(E) and f_j,ν,prod(E) 
being the fractions of diffusive and 
free streaming contribution to the emission
rates r_i,ν(E) and r_j,ν(E) at i 
and j points respectively,
which can be approximated by 
[We neglect the Pauli blocking and
the thermalization correction for this
calculation.]


    f_i,ν,diff(E)    ≈r̃_i,ν(E)/r_i,ν,diff(E), 
    
    f_j,ν,prod(E)    ≈r̃_j,ν(E)/r_j,ν,prod(E).


The value of β_ν is then estimated from
the ratio of the fluxes ∼Λ_ν/R^2 
an observer close to pole and equator would see
at a fixed distance R to the source
[We do not exactly choose the angles 
θ=0 and θ= 90 for two reasons:
first, the solid angle corresponding to θ=0
is small and the value of Λ_ν(θ=0 )
would be associated to a small region that would
not be representative of the flux at 
the pole. Second, the cosine dependence 
is only meant as an approximate trend 
of the flux between pole and equator.
We therefore set our fiducial
angles close to pole and
equator to θ=10
and θ=80 
respectively.]



    β_ν= Λ_ν(θ≈ 0 )/Λ_ν(θ≈ 90 ) - 1.


It is important to note that in 
the computation of β_ν we distinguish 
between neutrino species as the relevant 
neutrino-matter interactions differ between them
and they also have different decoupling surfaces
(see Figure <ref> 
for example), but for each species we do not  
consider different β_ν for 
different neutrino energies. Nevertheless,
our algorithm accounts for the different
contributions to β_ν from different
energies and therefore we assume Eq. (<ref>)
as indicative for a suitable estimate of 
the spectral neutrino density of Eq. (<ref>).



§ RESULTS


In this section we summarise
our results. We begin by applying the ASL to 
a 1D core-collapse supernova profile.
Subsequently we map this profile on a 3D
grid and apply our multi-D implementation of
the ASL, see Sec. <ref>.
Finally, we use the ASL to extract 
the neutrino physics from a neutron 
star merger remnant. Our tests 
are performed by taking snapshots of density,
temperature and electron fraction from
dynamical simulations <cit.> 
as a background on which we evolve the neutrino
quantities until they achieve a steady state.
For the core collapse supernovae 
tests, we use the Lattimer-Swesty Equation 
of State (EoS) <cit.> with nuclear
incompressibility K= 220 MeV, 
a standard choice widely used 
throughout the literature.
In light of the recent constraints
on the EoS from 
gravitational and electromagnetic
observations of 
GW170817 <cit.> 
for the binary merger remnant we 
consider the SFHo EoS, which provides masses in 
agreement with the highest neutron star masses  
measured until very recently <cit.>,
but in tension with the most recent 
measurement of a 2.17^+0.11_-0.10 neutron 
star <cit.>.
The neutrino transport is run with a 
spectrum of 20 geometrically increasing energy 
groups from 3 MeV to 300 MeV. We also explore 
for fixed energy interval the convergence 
of our results by changing the number 
of energy bins, and for fixed number of bins
we also increase and decrease 
the energy interval to test the 
dependence of our results 
on the energy spectrum.
Our ASL results are 
compared with a two moment scheme 
(M1) <cit.>.
In particular, for the core-collapse 
supernovae tests we compare with 
the M1 scheme of the spherically 
symmetric Eulerian hydrodynamics
code GR1D <cit.>.
Since we take the 1D core-collapse 
supernova snapshot 
from the dynamical simulations of
<cit.>, which uses the same ASL
described in Sec. <ref>, 
we perform a dynamical simulation
of the same progenitor with GR1D and 
take the outcome from the simulation at
the same time post-bounce to make 
a consistent comparison.
For the binary neutron star merger case
we compare with the 
M1 scheme implemented in the
Eulerian hydrodynamics code FLASH <cit.>.
Unlike the core-collapse case, 
we take a snapshot of a post-merger remnant from the 
simulations of <cit.>,
which uses a grey leakage scheme,
and apply both the ASL implementation
described in Sec. <ref>
and the M1 scheme in FLASH to it.
Applying a different transport approach
to the one adopted for obtaining
the snapshot introduces inconsistencies
that will be quantified.



 §.§ ASL in 1D core-collapse
supernovae




We take a snapshot at 275 ms post-bounce 
of the 15  progenitor of <cit.>,
whose dynamical evolution has 
been simulated with
the spherically symmetric hydrodynamics
code Agile <cit.>.
The radial profile has a variable resolution 
with radius ranging from ≲ 1 km to 
tens of km moving from 
the inner to the outer regions,
with a maximum radial extension of 6832 km.
In order to describe the profile properties, 
it is worth summarising the previous dynamics.
During the collapse phase the deleptonization 
of the iron core reduces the fraction of 
electrons Y_e in the core by producing 
electron neutrinos. 
Neutrinos freely stream out initially and 
therefore the total lepton number decreases.
Once neutrinos get trapped, 
the decrease in Y_e is compensated by the 
fraction of trapped neutrinos. At
core bounce a shock 
forms and propagates outward.
During the shock propagation, 
iron nuclei falling into the
shock are photo-dissociated into neutrons and 
protons. At this stage, neutrinos of all flavours
are largely produced by charged current interactions 
and pair processes. Dissociation of iron 
nuclei into nucleons, the ram pressure of the still 
infalling outer layers, and the energy losses due to
the neutrino burst when it surpasses the neutrinosphere,
cause the shock to stall. Neutrino
absorption behind the shock in the so called gain
region helps the shock to revive and in some cases
leads to the final explosion.
In the top row 
panels of Figure <ref> 
we show density, temperature and electron fraction as 
a function of the radius, while in
the second and third rows we show the 
number and energy of trapped neutrinos
per baryon, denoted
as Y_ν and Z_ν respectively.

We evolve such fractions until steady state is
reached while keeping density, 
temperature and electron fraction of the fluid
fixed. Once in equilibrium, the  rates of the
trapped neutrino  components vanish, which 
translates into u=e 
and Y_l= Y_e, 
see Eqs. (<ref>) and (<ref>). 
The distribution of  u and Y_l
(from a simulation using the flux factor of Eq. (<ref>))
is shown in Figure <ref>.
All the results we are going to show 
are referred to this equilibrium state.


  §.§.§ Neutrino rates



In Figure 
<ref> we show the specific
emission rate along the radial profile 
for different energies,
obtained from Eq. (<ref>) 
by multiplying for each energy
by E^2Δ E of the corresponding bin,
as well as the cumulative contribution from all energies. 
The main contribution to 
the neutrino emission comes from neutrino
energies below ∼ few tens of MeV for all 
neutrino species. In particular, 
all species show a decreasing 
contribution to the emission with 
increasing energy at low radii because of the 
∼τ_ν,tot^-2∼ E_ν^-2
dependence of the diffusion rate 
(equation (31) of <cit.>) 
that causes high energy neutrinos to have 
a lower diffusion rate and therefore 
to diffuse out less efficiently 
than low energy neutrinos.
The production rate dominates over 
the diffusion rate at large radii.
Given the emissivity dependence ∼ j_ν
(equation (30) of <cit.>), neutrino 
emission is governed by the low temperatures 
and is therefore suppressed at all energies.
Similarly, in Figure <ref> we show 
the specific absorption rate
for electron neutrinos and anti-neutrinos
(heavy-lepton neutrinos are not included
in the heating). Neutrino absorption is 
calculated locally for each
neutrino energy if the condition τ_en,ν(E,x) ≤ 1 is satisfied,
and depends on the local amount 
of available neutrinos, 
which in turn depends on both local
production and on neutrinos coming 
from innermost regions. 
Local neutrino production 
decreases as the radius increases,
as a consequence of the decrease
in temperature and density. 
In addition, the amount of 
neutrinos at a given radius 
coming from innermost 
regions is reduced both by absorption
occurring at smaller radii and
by the ∼ 1/R^2 dependence
in Eq. (<ref>).
Therefore, all energies show 
a decrease in the local neutrino
absorption with increasing distance
from the centre. We also notice 
in the same way of the emission rate 
that the larger contribution 
to the heating rate comes from neutrinos of ∼
few tens of MeV, as a result of
the dependence of the neutrino 
absorption on the emission via
the number rate l_ν(E,x).



  §.§.§ Flux factor and heating

In Figure <ref> 
we compare the inverse
of the flux factor for the two prescriptions of Eqs. (<ref>) and 
(<ref>), 
labelled as old and new 
respectively, for electron neutrinos and anti-neutrinos
at different energies, 
as it enters Eq. (<ref>) and it 
therefore affects the local heating.
The largest differences are at energies 
≳ 40 MeV with discrepancies 
up to ∼ 25-30%.
Such differences are expected since the flux
factor from Eq. (<ref>)
depends on the optical depth rather than 
the radius as in Eq. (<ref>), and
therefore the trend in Figure 
<ref>
resembles the trend of the optical depth.
Nevertheless, as we have seen in 
Figure <ref>
the largest contribution to the
heating comes from neutrinos 
of energies of few tens of MeV. 
We therefore do not expect such 
differences to contribute 
sensitively to the global neutrino 
luminosities for our snapshot calculations.
 In fact, in the bottom panels of 
Figure <ref> we show 
the quantity -Q_ν^k=1 in 
units of 10^20 erg g^-1 s^-1
calculated
from Eq. (<ref>),
which is a measure of the 
local heating rate. Overall, 
we see a very good agreement between  
the old (blue line) and new (red line) 
flux factor prescriptions, with 
differences at the order of a few percent
in the region where heating gets important 
(-Q_ν^k=1 > 0).
More detailed quantification of the heating
with our new flux factor prescription
during full dynamical simulations will be
the subject of future work.

As an additional test, we show in Figure <ref>
a comparison at the same time post-bounce of 
the same quantity -Q_ν^k=1
provided by a dynamical 
evolution of the same progenitor 
starting at the onset of collapse
performed with GR1D (black-dashed line). 
Compared to
our test with the ASL, the net rate 
from GR1D provides less cooling 
(less negative -Q_ν^k=1) in the 
ranges ∼ 30-50 km and ∼ 65-80 km,
with a larger one only in the range
∼ 50-65 km for electron neutrinos.
This is the result of a dynamical
evolution with a different neutrino transport, 
where at the same time after bounce the structure
of the star shows significant
differences with respect
to the ASL run, clearly visible in the density and 
temperature profiles on the top panels of Figure <ref>. 
In particular, the GR1D profile has a colder and 
less compact layer within ∼ 30-50 km where
most of the emission occurs.
Note also the different location of the shock, 
which is located at ∼ 105 km for the GR1D case. 
On the other hand, we do not see a 
difference in the peak heating
rate near 100 km, with -Q_ν^k=1∼ 2.5·10^20 erg 
g^-1s^-1 for both runs.



  §.§.§ Neutrino luminosities and average-rms energies

Given the specific net rates we calculate the
total neutrino luminosities and average 
energies given by Eqs. (<ref>)
and (<ref>).
The results are shown in the first row of 
Table <ref> for both the old and 
the new flux factor prescriptions.
The latter reduces the luminosities 
compared to the former by ≲ 5%. 
The average energies are 
less affected
with differences ≲ 1%.
For completeness, we also 
calculate the rms neutrino energies 
from Eq. (<ref>) and find 
E_rms,ν_e= 14.21 MeV,
E_rms,ν̅_e= 17.07 MeV, 
E_rms,ν_x= 25.33 MeV, 
in agreement with <cit.>. We
further show the results obtained 
by performing the same calculations with GR1D in the
last row of Table <ref>. 
The electron neutrino and anti-neutrino
luminosities are 
lower by ≲ 10% compared 
to the ASL runs due to the combination of an 
overall weaker neutrino cooling 
and comparable heating. 
Heavy lepton neutrinos have instead 
∼ 7% lower luminosity. 
By looking at the average energies,
differences are at the level of
≲ 5% for electron neutrinos 
and anti-neutrinos, and of ∼ 7%
for heavy-lepton neutrinos.
A similar trend is seen in 
the rms energies.
Overall, beside the different grid
setups, these discrepancies
are a consequence of the usage of
different neutrino transport schemes.



  §.§.§ Number of energy bins

We finally determine the 
number of energy bins that 
are needed for trustworthy results for the luminosities and
average energies. We first vary the number of 
energy bins in the fixed energy 
interval [3, 300] MeV. 
Table <ref> summarises our 
findings. A number of energy bins ≲ 10 causes 
sensitive deviation from a regime of convergence
that is visible at larger numbers
(which includes our preliminary
choice of 20 energy groups in [3, 300] MeV). In particular,
we notice a decrease in the luminosities and to 
a minor extent in the average energies, as a result 
of poorly resolved energy-integrated emission and 
absorption rates. We therefore choose 20 energy bins 
and vary the energy interval. In this way, variations
in the simulation outcome  by a reduction 
of the spectrum size would provide information 
on those energy ranges that 
are too small regardless
of the number of energy bins. Moreover,
we can assess whether the regime of convergence 
with 20 bins is satisfied for wider intervals of
energies or if it is strictly bound to certain 
energy ranges. Results are shown in Table 
<ref>. 
We see that cutting the energy 
spectrum [3, 300] MeV at high energies leads to
a significant decrease 
of L_ν_x^1 that starts
appearing from spectra with upper energies below
∼ 75 MeV and that worsen by up to about 
50% decrease for the smallest
range [3, 30] MeV.
The average energy ⟨E_ν_x⟩ is
also reduced to 14.37 MeV. This reduction is due 
to the lack of contribution coming from energies 
≳ 40 MeV, that can still be relevant 
to the emission from heavy-leptons 
(see left panel of Figure <ref>).
In contrast, electron neutrino 
and anti-neutrino values remain 
almost stationary. In the same way, cutting the 
energy spectrum at low energies from 3 MeV has a
strong impact on the luminosities, with a reduction
of up to 70-80% for electron neutrinos and 
anti-neutrinos in the case of the [25, 300] MeV
range. The average energies show the 
opposite trend, 
increasing as a result
of the high energies 
giving contribution to the
luminosity L_ν^1
via Eqs. (<ref>),
(<ref>) and
(<ref>) and thus affecting the
mean of Eq. (<ref>). 
On the other hand, we notice a convergence in
the values of L_ν^1 and  
⟨E_ν⟩ 
for spectra 
spanning a range from 3 
MeV to hundreds of MeV.
In particular, no sensitive 
variations are seen
by extending the interval 
of energies above 300 MeV.
The smallest interval above which we start 
seeing convergence is [3, 75] MeV. 
However, it is important to specify 
that we are here basing our 
convergence tests by just looking at
global neutrino luminosities
and neglecting the convergence
of the fractions of neutrino 
trapped components Y_ν and 
Z_ν which are 
crucial in the modelling of the
diffusion while performing
dynamical simulations.
Such convergence requires
neutrinos of energies 
at least equal to
the neutrino chemical potential
in the core (set by the 
beta-equilibrium condition),
i.e. ≳ 100 MeV.
Moreover, we want
to stress that our
convergence tests are performed 
over a snapshot, but the way the 
selected energy interval affects 
the emission can change for different
stages of a dynamical evolution.
Summarising,
a number of 20 energy bins for the neutrino
spectrum is a reasonable choice 
despite the assumed interval of energies, provided that 
there are neither cuts in the spectrum at low 
energies nor at energies that would exclude 
neutrinos of ≳ 100 MeV. 








 §.§ ASL in 3D applications




  §.§.§ Comparison in spherical symmetry between 1D and 3D

To scrutinise our multi-D implementation,
we start by taking the snapshot analyzed in 
Sec. <ref>, map the initial 
data on a 3D grid with uniform resolution
of 1 km and evolve the neutrino 
transport part until equilibrium. 
Instead of mapping the whole radial
profile which extends up to 6832 km at 
densities below 10^6 g/cm^3 and 
temperatures of ∼ 0.1 MeV, we take 
the profile information only up to 150 km
and neglect the remaining part to save 
computational time. Indeed, beyond such 
distances densities and temperatures are 
low enough that the neutrino contribution
to the outcome of the simulations is 
negligible, ≲ 1% in the total luminosities 
and average energies. We choose the energy interval
[3, 300] MeV, see Sec. <ref>. 
We show the mapped initial conditions 
on the y=0 plane in the bottom 
panel of Figure <ref>.


3D optical depth

To test our 3D implementation 
of the optical depths, we 
first calculate the optical
depth on the 3D grid
as explained in Sec. <ref>.
We then create a 1D profile 
equivalent to the one 
used in Sec. <ref>
but with uniform resolution 
of 1 km, where we calculate
the optical depth by 
doing a simple integration
over the radial path. We finally map
such optical depth on the 
3D grid and calculate 
the relative error ε_ν=
τ_1D-τ_3D/τ_1D
on the y=0 plane.
In the left panel of 
Figure <ref>
we show our result.
As reference case we
take the total optical 
depth for electron
neutrinos of energy 
E_ν_e= 10.08 MeV, 
the other energies 
and neutrino species show
a similar behaviour. 
The differences between the
1D and the 3D calculations 
are at the level of 
∼ 8% at the most, with the
3D implementation providing
larger values overall. 
To get an idea of the distribution of the 
neutrino surfaces at different energies, we
show the location of the total and energy 
neutrinosphere radii on the density map at 
y=0 in Figure <ref>. 
The selected set of energies are 
the same used to show the 
plots in Figures <ref> and
<ref>. 
The neutrino surfaces are almost
perfectly spherical, suggesting
that our optical depth algorithm 
overall captures the actual path that 
minimises the optical depth
and that neutrinos preferentially
cross, i.e. the radial path.
Obviously the larger the energy
the larger the radius of 
the neutrinosphere, since  
τ_ν∼ E_ν^2. Comparing 
between species, heavy-lepton 
neutrinos have smaller energy neutrinospheres
because the only interactions where
they exchange energy with the
fluid are pair processes and bremsstrahlung.
On the contrary, elastic scattering on nucleons
and nuclei makes the total neutrino surfaces 
comparable with the other species. Electron neutrinos and anti-neutrinos show similar energy and total neutrinospheres 
as a result of the comparable amount 
of emission and absorption interactions involving both
species at this time of the post-bounce phase. Moreover, for each of these neutrino species we notice comparable 
radii of total and energy neutrino surface at 
each energy, indicating that 
neutrino emission and absorption reactions, 
efficient in thermalizing neutrinos,
provide also an important opacity contribution
to the total optical depth.
Overall, the total neutrino surfaces extend from 
∼ 27 km to ∼ 121 km for electron
neutrinos, from ∼ 25 km to ∼ 121 
km for electron anti-neutrinos, and from 
∼ 25 km to ∼ 118 km for heavy- 
lepton neutrinos. Accordingly, the energy 
surfaces extend from ∼ 26 km to ∼ 121
km for electron neutrinos, from ∼ 24 km 
to ∼ 121 km for electron anti-neutrinos,
and from ∼ 23 km to ∼ 69 km for 
heavy-lepton neutrinos.


Heating

Calculation of the absorption rates on the
grid is done by applying Eq. (<ref>).
Computation of the neutrino density is
performed by means of Eq. (<ref>) 
with β_ν= 0 (i.e. Eq. (<ref>)) 
and of Eq. (<ref>).
Given the spherical symmetry 
(see Figure <ref>) 
the value of l_ν(E,R) over the
grid is approximately 
recovered by integrating 
Eq. (<ref>) over
a reference radial path from the origin,
and by mapping the obtained values on the path
over the rest of the grid.

Figure <ref>
shows the heating rate 
for electron neutrinos and anti-neutrinos on 
the plane y=0 with the new flux factor prescription 
in units of 10^20 erg g^-1 s^-1, resembling the values 
from Figure <ref>. Moving outward from 
the centre, the net rate decreases and reaches
its minimum at ∼ 50-60 km from the centre,
then increasing and getting to positive values
where neutrino heating dominates between 80 km 
and 120 km. Unlike Figure <ref> we 
notice a slightly lower minimum that goes below
-17.5·10^20 erg g^-1 s^-1
and -15·10^20 erg g^-1 s^-1 for electron neutrinos and 
anti-neutrinos respectively, and a 
larger maximum above ∼ 2.5·10^20 erg g^-1 s^-1 for the heating.
This is expected because 
the original
profile of Sec. <ref>
has decreasing
resolution with increasing distance
from the centre, leading to larger
optical depths in the transition
between the optically thick and 
the optically thin regime. 
This is clearly visible in the right panel
of Figure <ref>, where
we plot the total optical depth
for electron neutrinos of energy 10.08 MeV
(taken as reference)
calculated by integrating along the 
profile of Sec. <ref>
(blue line), together with values of
the optical depth calculated
along several paths
of the 3D grid. In the range
∼ 40-60 km the 1D 
optical depth is ≲ 10% larger
than the 3D ones.
The lower values from the 3D
calculations lead to  
lower diffusion timescales
and therefore to 
stronger emission and
absorption rates.


Neutrino luminosities and average-rms energies

Given the heating rate, we
calculate the total neutrino luminosities 
and average energies. The results are shown 
in the third row of Table <ref>.
Differences with respect to the 1D case
(Sec. <ref>) 
are of the order of ∼ 7-8% 
for the electron neutrino and anti-neutrino 
luminosities, and of ∼ 1% for the
heavy-lepton neutrinos. 
Average energies all differ by 
≲ 1%. Considering the results
with the old flux factor prescription,
differences in the luminosities compared to the 
1D case are of the order of ∼ 10% for 
electron neutrinos and anti-neutrinos and 
∼ 1% for heavy lepton neutrinos. 
In addition, 
we repeat the calculation
for the case of a 1D radial 
profile with uniform 
resolution of 1 km. 
Results are shown in the second 
row of Table <ref>.
Errors in the luminosities between 
the 3D and the 1D implementation 
reduce to ∼ 3-4% and to
∼ 1-2% with the old 
and new flux factor prescription 
respectively, confirming that the 
resolution contributes as 
source of variability to the 
outcome of the simulations. 
As we have similarly seen already 
in Sec. <ref>, 
the choice of the new flux 
factor in our 3D simulations provides 
luminosities ≲ 6-7% lower respect 
to the old prescription, while average energies are 
less affected with discrepancies 
of ≲ 2%. 
Overall, we state that the
modification adopted to the standard ASL of 
<cit.> perform well in the context
of core-collapse supernova, although more 
precise assessments require full dynamical evolutions. 



  §.§.§ Neutron star merger remnants in 3D


As a final 3D test we apply our new scheme
to the remnant of a neutron star merger. We start 
from a snapshot of a merger remnant (t= 38 ms after merger;
<cit.>) that has been obtained using the SPH method
with the grey leakage scheme of <cit.>. For
general reviews of the SPH method we refer to recent reviews <cit.>.
We use the snapshot as a background 
on which to evolve the 
neutrino properties until a steady state is reached.
Unlike the core-collapse supernovae case, 
we do not have any
information on the neutrino trapped components at this time and therefore we perform our tests by assuming
Y_ν=Z_ν=0 as initial condition for neutrinos. 
Although strictly speaking 
inconsistent as in the configuration
trapped neutrinos 
should already be present, 
our choice of such initial condition 
is justified by several arguments.
First, the role of the neutrino
trapped components is mainly 
important when doing full 
dynamical evolutions and in 
particular it has been shown recently 
<cit.> that they only
marginally affect
the thermodynamical properties 
of the remnant, 
and only close to the merger 
time where neutrinos 
are produced in the first place. 
Second, setting Y_ν=Z_ν=0 
implies a remnant configuration 
out of equilibrium, and given the
large temperature dependence of 
the neutrino-matter cross sections
<cit.>
the system rapidly reaches
a new state by refilling the 
neutrino fractions over a
timescale Y_ν/Y_ν≲ 10^-6 s,
which is much less than the typical 
dynamical timescale of the 
remnant ∼ (Gρ̅)^-1/2≈ 2·10^-4 s.
Third, absorption under optically thin 
conditions (that we are modeling here) 
is led by non-trapped neutrinos. For all 
the above reasons, we do not expect the 
modeling of the trapped neutrinos to 
have a significant impact on our calculations.
The neutron stars have been discretized 
with N ∼ 1 million particles and the initial
conditions are mapped on a 3D grid whose 
borders are defined when densities go
below ∼ 10^8 g cm^-3. 
Figure <ref> shows the
density, temperature and electron fraction 
on the planes y=0 and z=0. 
The central object has a mass 
≈ 2.3, 
densities ρ∼ 3·10^14g cm^-3, temperatures T∼ 20 MeV 
and electron fractions Y_e∼ 0.1. 
Around it is a torus with Y_e < 0.1, an 
inner region with densities ρ∼ 10^12-10^13g cm^-3 and temperatures 
T∼ 3-5 MeV, and an outer region 
with densities ρ∼ 10^10-10^11g cm^-3 and temperatures T∼ 1 MeV. 
Electron fractions Y_e∼ 0.2-0.3 are 
located along the low-density polar region.




3D optical depth

Computation of the neutrino properties 
requires the calculation of the optical
depth on the grid as explained in 
Sec. <ref>.

Figure <ref> 
shows the location of the total and 
energy neutrino surfaces for the 
same sets of energies used for the core-collapse 
supernovae tests 
(see Figures <ref> and <ref>). Overall, we see an 
agreement with the distribution of the surfaces 
described in <cit.>.

As already noticed for the core-collapse supernovae case, the higher the energy of neutrinos the more extended 
the neutrino surface. Accordingly,
the heavy-lepton neutrinos
have less extended energy 
than total neutrino surfaces
because they can only exchange
energy by pair processes and 
bremsstrahlung.
Elastic scattering on nuclei and nucleons make the 
heavy lepton total neutrino surfaces comparable with 
the ones of the other species. Electron neutrinos again 
show comparable energy and total neutrino surfaces due to
the neutron rich environment that favours 
energy exchange in processes like neutrino absorption on neutrons. However, we notice that compared to the core-collapse 
supernova case electron 
anti-neutrinos show less extended
total and energy neutrino surfaces compared
to the corresponding electron neutrino ones as
a consequence of the low abundance of free protons.



Heating

The heating is modelled by Eqs. (<ref>) and (<ref>) and by estimating β_ν as described in Sec. <ref>. 
To save computational time, we limit our 
transport calculation to those regions
where density is above 10^9 g cm^-3.
Indeed, we find that the contribution 
at lower densities affect the transport 
quantities by less than 1%.
For the computation of l_ν(E,x) we
create a 1D profile of 1 km of resolution where 
to each bin of radius R_bin we assign a neutrino 
emission by summing up the contribution
of all grid points with radial distance from
the centre within that bin. We then solve
Eq. (<ref>) and assign the 
same l_ν(E,R_bin) to all
grid points inside that bin.
In this way we create a spherically
symmetric neutrino emission by coupling 
fluid points from the torus with 
fluid points along the poles, 
and we then leave to 
β_ν the task of 
approximately recovering 
the degree of anisotropy of the system 
and consequently the degree of decoupling 
between points at same distance 
from the centre but at different polar angles.
The determination of the exponent b in
Eq. (<ref>) is 
performed by comparing with an M1 
calculation of the heating from FLASH. 
In particular, we find b=8 
to overall best recover the 
electron neutrino 
contribution to the heating
(which constitutes more than 
50% of the total)
and we therefore assume the same
value for the corresponding 
anti-neutrinos as well.

In Figure <ref> we show
the angular dependence of the
neutrino flux, i.e. 
Λ_ν(θ) vs θ, 
calculated from Eq. (<ref>)
with Δθ= 10
[We have tested different 
Δθ and no appreciable
variations appears in the computation 
of β_ν, therefore we set 
Δθ= 10.]. 

We notice that the modulation of the
flux with the polar angle is more pronounced
for electron neutrinos than for electron 
anti-neutrinos.
This is due to the fact that 
the neutron rich torus 
is more opaque to electron
neutrinos than to electron anti-neutrinos.
Therefore, the electron
neutrino flux points mostly along
the z-direction. 
In terms of neutrino emission,
the largest contribution to the 
cooling for the electron neutrinos
(∼ 10^21 erg g^-1 s^-1) 
is confined 
within radii ≲ 20 km from the centre,
and the remaining subdominant part 
(∼ 10^19-10^20 erg g^-1 s^-1) occurs
inside the torus.
We find that 
electron anti-neutrinos in contrast
are mostly emitted in the torus
(cooling ∼ 10^20 erg g^-1 s^-1)
and no relevant emission is found 
at radii ≲ 20 km from the centre.
We obtain β_ν_e≈ 16 and 
β_ν̅_e≈ 2.
In the first row of Figure <ref>
we show 2D maps of the 
resulting heating rate on the plane y=0.
In the second row we show
the same maps for 
calculations performed with the M1 
scheme implemented in FLASH.

The major contribution to the heating
is located within ∼45 from the
pole for both species. 
The largest differences of the ASL 
compared to M1 are in the region with
θ≲ 45, 
and in particular 
above the central object
at θ≲ 15,
where the electron
neutrino and anti-neutrino heating
are respectively lower by a factor of 
∼ 1.5 and larger by 
a factor of ∼ 2.
Moreover, unlike the M1 implementation our
ASL algorithm provides a residual electron
neutrino heating 
(∼ 10^19-10^20erg g^-1 s^-1) at 
θ∼ 60-90 in regions 
where τ_tot≤ 1.
We notice that both ASL and M1
provide an electron neutrino 
heating larger by one order of 
magnitude compared to electron 
anti-neutrinos. This is due to two 
effects: first, our snapshot
calculations show an electron
anti-neutrino cooling which is
at the most ∼ 10^20erg g^-1 s^-1,
i.e. one order of magnitude lower
than the electron neutrino one,
which reaches ∼ 10^21erg g^-1 s^-1
close to the central remnant
where the anti-neutrino emission
is negligible in comparison.
Second, the neutron rich environment
favours the electron neutrino
absorption over the corresponding
anti-neutrinos.



Neutrino luminosities and average-rms
energies

The overall low anti-neutrino cooling
leads to almost equal
electron neutrino and 
anti-neutrino total luminosities and 
average energies, contrarily to
what is typically expected from
merger simulations
<cit.>. In particular, 
a summary of these values is reported 
in Table <ref>. The reason
for this is most likely due to
the application
of a different
neutrino transport 
than the one
adopted for generating 
the snapshot.
In fact, we further calculate
the neutrino cooling with our ASL 
for a similar snapshot 
of the merging of two 1.4 neutron
stars taken from 
the simulations of 
<cit.>
at ∼ 66 ms
after the first contact of the two stars. 
The ASL of <cit.>
is similar to ours, and 
in particular is spectral.
In this case we get
L_ν_e^cool= 1.45·10^52 erg s^-1,
L_ν̅_e^cool= 1.85·10^52 erg s^-1,
L_ν_x^cool= 1.34·10^52 erg s^-1 ,
⟨E_ν_e⟩= 10.81 MeV, ⟨E_ν̅_e⟩= 14.97 MeV, ⟨E_ν_x⟩= 16.00 MeV,
E_rms,ν_e= 12.02 MeV,
E_rms,ν̅_e= 16.77 MeV, 
E_rms,ν_x= 19.74 MeV,
that is, the largest cooling
is from the electron anti-
neutrinos being  
the most luminous, and we 
also clearly recover the 
expected hierarchies, i.e.
L_ν̅_e^cool > L_ν_e^cool > L_ν_x^cool,
⟨E_ν_x⟩ > 
⟨E_ν̅_e⟩ >
⟨E_ν_e⟩, 
E_rms,ν_x > E_rms,ν̅_e
> E_rms,ν_e
[Although we
recover the expected hierarchies,
we do not find 
the same values of 
<cit.> as 
they perform a dynamical 
simulation assuming neither
blocking nor 
thermalization corrections 
in Eq. (<ref>).].
In particular,
we observe that the configuration at
∼ 66 ms shows larger densities 
and temperatures at radii 
≲ 20 km from the centre
compared to the ∼ 38 ms one,
producing a dominant
electron anti-neutrino 
cooling ∼ 10^21 erg g^-1 s^-1
(see Figure <ref>
for a comparison done with ASL 
along θ=0)
that is instead missing in
the configuration at ∼ 38 ms.
Moreover, the two
snapshots have been 
previously evolved
with different neutrino transport
schemes: the one at 
∼ 66 ms with the scheme of
<cit.> which is similar
to our version of the ASL,
the ∼ 38 ms one with the
grey scheme of <cit.>.
Applying a different
neutrino transport than
the one used for generating
the snapshot introduces
inconsistencies that are
likely to be the reason
for our luminosity values.
From Table <ref> we notice that
the ASL provides electron 
neutrino and anti-neutrino luminosities
lower by ∼ 35% and ∼ 25%
respectively in comparison to M1. 
Average-rms energies 
agree within 5-15%.
The lower luminosity 
values from the ASL can be 
due to several reasons. First, the 
ASL has excess heating
at θ∼ 60-90
for the electron neutrino and 
at θ∼ 15 for the
electron anti-neutrinos compared to M1.
Second, 
we have kept the values of 
the ASL parameters
entering Eq. (<ref>) to the
ones calibrated for core-collapse
supernovae simulations.
New calibrations are crucial as 
they might impact 
the neutrino cooling and 
consequently the amount of neutrino
heating. In particular, the comparison with
M1 suggests a lower α_ν,  blk.
This can be explained in the following way.
While in core-collapse supernovae 
the quasi-spherical symmetry implies that
neutrinos move preferentially along the
radial direction, in binary mergers
neutrinos have more directional 
freedom in escaping the system
and consequently the blocking is
expected to be less effective.
Third and above all, 
a more consistent comparison
between two neutrino transport
approaches would require dynamical 
evolutions rather than 
snapshot calculations.
Putting together 
all these uncertainties 
we find our
ASL-M1 luminosity discrepancies, 
lower than a factor of 2,
a promising initial step
toward future
developments.
We also want to stress that
the choice we have made
for the flux modulation ∼cos^8(θ)
is only meant to be a preliminary step
toward explorations in full dynamical
evolution of binary mergers.
Finally, we stress the fact that
we have based our heating calculations
on the comparison with a moment
approach with analytical closure.
However, as pointed out by <cit.>
the M1 closure can overestimate
the neutrino density by up to ∼ 50%
along the polar regions,
thus limiting the possibility of
properly modeling the neutrino driven-winds.
Future simulations of binary mergers
will definitely need improvements
in moment schemes as well for 
more reliable assessments of the
neutrino physics in such systems.




§ SUMMARY


In this paper we have presented an extension of the
Advanced Spectral Leakage scheme originally introduced
by <cit.>. Our main goal was to  adapt the scheme
so that it can be conveniently  used for neutron star merger
simulations. The main advantage compared to simpler leakage
schemes is that the ASL includes neutrino heating processes
and can therefore, with a reasonable computational effort, also
model the emergence of neutrino-driven winds in a 3D merger
simulation.
The main novelty compared to the original approach is the usage
of an optical depth-dependent flux factor, and a modification
to the equation of the neutrino density in the semi-transparent regime,
both designed for the multi-dimensional modeling of compact
binary mergers.

We scrutinized the new scheme on the case of a 15 
core-collapse supernova snapshot taken from <cit.>
at 275 ms after bounce. For this spherically symmetric case we
first tested in 1D our new flux factor and found that it agrees well
with the original choice (≲ 2% for the average energies
and ≲ 5% in the neutrino luminosities). As a further
1D test, we have compared the new scheme to the M1 implementation of the GR1D code <cit.>
 and here we found 
 agreement, beside differences 
 arising from the usage of 
 different transport schemes. 
We also mapped the 1D case
onto a spherically symmetric, but three-dimensional grid. Here the agreement is slightly
worse, but overall still very good: ≲ 2% for the
average energies and ≲ 7% for the neutrino luminosities.

We have finally explored the ASL for an SPH snapshot
of a 1.4-1.4  binary neutron star merger. As a reference
we compared against the results obtained with a M1 scheme that
is implemented in FLASH <cit.>.
Here, in ASL the anisotropy in the neutrino fluxes is taken into
account by an anisotropy parameter β_ν
which is estimated from the  ratio between the neutrino fluxes
at pole and equator in a similar way as it has been
done in <cit.>. The neutrino density is modelled
∝cos^b(θ), where the value of b= 8 
is obtained by a comparison of the neutrino heating rate with the
M1-FLASH results. Overall, we find good agreement in the 
neutrino heating distribution between both approaches. 
The average energies
agree within 5-15%. The ASL total luminosities 
are lower by 25-35% compared to M1. 
This discrepancy may suggest that some of the free parameters
need to be calibrated more specifically for the case of binary 
compact mergers. While specific questions may require more 
sophisticated neutrino
transport methods, we are confident that this enhanced ASL scheme
delivers reasonably accurate bulk neutrino properties. This will be applied and further tested in future dynamical simulations.

In the tests shown, relativistic effects 
were neglected.
The inclusion of general relativity enhances the gravitational well, leading to more compact 
and hotter remnants. The energy of the radiation
field measured by an observer at a given distance
from the emission region is redshifted as neutrinos 
climb out of the gravitational well.
In addition, for relativistic fluid motions 
the neutrino energy is Doppler shifted 
depending on the relative velocity
between the source and the observer. 
Simulations of core-collapse supernovae
including relativistic effects
show an increase in 
the neutrino luminosities and average energies
with respect to the Newtonian case, as well
as larger energy-deposition rates by 
neutrino absorption
<cit.>.
In the context of compact binary mergers,
a larger heating can affect both
the amount of neutrino-driven
wind ejecta and  
the reprocessing of the 
electron fraction in the ejecta.
At last, it is worth keeping 
in mind that 
most of the 
neutrino emission in a 
merger remnant comes 
from the torus.
The orbital 
velocities of the matter inside
of it can be mildly 
relativistic (e.g. <cit.>).
As a consequence, the angular distribution
of the neutrino fluxes emitted at the 
neutrino decoupling surfaces can be
sensitively affected 
by relativistic beaming. 
In particular, the neutrino emission 
at the neutrino decoupling surface
seen by the observer 
will be concentrated in a cone   
directed toward the direction of 
the fluid motion,
rather than being isotropic as seen in 
the rest frame of the fluid.
This can reduce the amount of 
received flux
along the poles with 
respect to the Newtonian case.
The inclusion of all these effects will
be the subject of future work.



§ ACKNOWLEDGEMENTS

This work has been supported by the Swedish Research 
Council (VR) under grant number 2016- 03657_3, by 
the Swedish National Space Board under grant number 
Dnr. 107/16 and by the research environment grant 
"Gravitational Radiation and Electromagnetic Astrophysical
Transients (GREAT)" funded by the Swedish Research 
council (VR) under Dnr 2016-06012. 
We would like to 
thank Sean Couch for major contributions to the
development of FLASH.
Moreover, we gratefully 
acknowledge support from COST Action CA16104 
"Gravitational waves, black holes and fundamental 
physics" (GWverse), from COST Action CA16214 
"The multi-messenger physics and astrophysics of 
neutron stars" (PHAROS), and from COST Action MP1304 "Exploring 
fundamental physics with compact stars (NewCompStar)".

The simulations were performed on resources
provided by the Swedish National Infrastructure for
Computing (SNIC) at PDC (Center 
for High Performance Computing) and NSC
(National Supercomputer Centre) 
and on the resources of the Northern German
Supercomputing Alliance (HLRN).






mnras



















