language=Python, basicstyle=, basewidth=0.53em,0.44em, numbers=none,
tabsize=2, breaklines=true, escapeinside=@@, showstringspaces=false,
numberstyle=, keywordstyle=, stringstyle=, identifierstyle=, commentstyle=,
emphstyle=, frame=single, rulecolor=, rulesepcolor=, literate = 1 as  as 3
.set.set4

TTK-23-001, TTP23-012

PEANUTS :
a software for the automatic computation of solar neutrino flux and its
propagation within Earth

Tomás E. Gonzalo†,1 and Michele Lucente⋆,2
1 Institute for Theoretical Particle Physics (TTP), Karlsruhe Institute of
Technology (KIT), 76128 Karlsruhe, Germany
2 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen
University, D-52056 Aachen, Germany
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present PEANUTS (Propagation and Evolution of Active NeUTrinoS), an open-
source Python package for the automatic computation of solar neutrino spectra
and active neutrino propagation through Earth. PEANUTS is designed to be
_fast_ , by employing analytic formulae for the neutrino propagation through
varying matter density, and _flexible_ , by allowing the user to input
arbitrary solar models, custom Earth density profiles and general detector
locations. It provides functionalities for a fully automated simulation of
solar neutrino fluxes at a detector, as well as access to individual routines
to perform more specialised computations. The software has been extensively
tested against the results of the SNO experiment, providing excellent
agreement with their results. In addition, the present text contains a
pedagogical derivation of the relations needed to compute the oscillated solar
neutrino spectra, neutrino propagation through Earth and nadir exposure of an
experiment.

###### Contents

  1. 1 Introduction
  2. 2 Solar neutrino flux
    1. 2.1 Neutrino survival probability at Sun surface
    2. 2.2 Neutrino propagation from the Sun
  3. 3 Neutrino propagation through Earth
    1. 3.1 Probability of transition through matter
    2. 3.2 Evolutor operator
    3. 3.3 Earth matter regeneration
  4. 4 Neutrino propagation Hamiltonian
    1. 4.1 Perturbative expansion of the neutrino propagation Hamiltonian
  5. 5 Time integration and exposure
  6. 6 Quick start guide
    1. 6.1 Validation
  7. 7 Conclusions and Outlook

##  1 Introduction

Solar neutrinos represent the most abundant source of neutrinos on Earth, with
a flux of the order of $~{}6\times 10^{10}\ \text{cm}^{-2}\ s^{-1}$ [1]. Even
though only a small fraction of this flux is actually detectable, due to many
production channels resulting in neutrinos with energies below the typical
experimental detection thresholds [2], solar neutrinos provide an invaluable
source of information for the study of neutrino properties, solar dynamics and
Earth internal structure. Historically, solar neutrinos provided the first
hints for the non-conservation of leptonic flavours in neutrino propagation
[3], when in 1968 the Homestake experiment reported far less solar electron
neutrino events [4] than the number expected from the recently developed solar
models [5, 6, 7, 8, 9, 10, 11, 12, 13]. This so-called _solar neutrino
problem_ was only solved in 2001, when the SNO collaboration released [14] the
measurement of the total flux of active 8B neutrinos, that resulted in a close
agreement with solar models predictions and implied that the origin of the
discrepancy had to be tracked to new physics effects in the neutrino sector.
It is now firmly established that the deficit is due to neutrino oscillations
[15, 16, 17], with oscillation parameters as inferred by solar neutrino
experiments well in agreement with data from other neutrino sources [18, 19,
20].

Nowadays, solar neutrinos continue providing invaluable information, being the
only known probes that can directly test the interior structure and dynamics
of the Sun, thereby providing strong constraints on the solar models
parameters [21, 22, 23], or probe the internal solar dynamics on long time-
scales [24]. More exotic production mechanisms such as neutrino emission
during intense solar flares [25, 26], neutrino production from cosmic-rays
scattering in the solar atmosphere [27] and neutrinos from dark matter
annihilation in the Sun [28] can also be probed. On the other hand, the
theoretically established solar models allow to study and constrain neutrino
properties, from standard oscillation parameters [29, 30] to more hypothetical
scenarios such as non-standard interactions [31] or finite magnetic moment
[32].

In all these scenarios, an important wealth of statistical information is
encoded in the energy spectra for the individual neutrino flavours. Due to
neutrino oscillations, however, the spectra on Earth generally differ from the
ones at neutrino production, and a proper account of the oscillation dynamics
is mandatory in order to compare a theoretical model with data. The
oscillation dynamics is in general non-trivial, including three different
regimes and corresponding matter effects [33, 34]: propagation within the
slowly-changing111 _Slow_ or _fast_ here refers to the importance of variation
of matter density over a neutrino oscillation length scale. matter density
within the Sun, propagation in vacuum between Sun surface and Earth, and
propagation within Earth featuring a fastly-evolving matter density profile.

In this work we present PEANUTS (Propagation and Evolution of Active
NeUTrinoS), an open-source, fast and flexible package to compute the neutrino
oscillation dynamics in all the above-mentioned regimes. The emphasis in
developing the software has been put on both _performance_ and _flexibility_ :
PEANUTS computes the coherent neutrino propagation inside Earth analytically,
completely removing the need for time-consuming numerical integrations.
Moreover, the user can input any arbitrary solar model, as well as any custom
Earth matter density profile, and simulate experiments at any latitude and
underground depth. The software can perform the full chain of computations to
simulate the expected neutrino spectra for a given solar model, Earth matter
density and detector location, or its modules can be called individually, to
compute for instance the evolved neutrino state after Earth crossing given an
arbitrary initial (coherent or incoherent) state, or the solar angle
distribution for an experiment taking data between two arbitrary days of the
year. PEANUTS 222https://github.com/michelelucente/PEANUTS is provided as an
open-source Python package, under the GPL-3.0 license. It can be run
standalone in various different ways, as well as interfaced from other
frameworks. Details of the various software requirements and dependencies, as
well as instructions on how to use PEANUTS can be found in Section 6, whereas
the specific functions and classes that perfom the computations will be
scattered throught whenever the physical quantity computed is introduced. We
have extensively validated our implementation against the results of the SNO
experiment [35, 36], details of which will also be located where useful,
including a final comparison of likelihood contours for our simulation of the
experiment with the probabilities computed by PEANUTS.

This document is thus structured as follows. Section 2 describes the
theoretical background for the computation of the solar neutrino flux and the
propagation of the neutrinos from the Sun. In Section 3 we describe the effect
of Earth regeneration of the neutrino flux, i.e. the propagation of neutrinos
through the Earth, followed by a detailed explanation of our fast
approximation for the neutrino propagation Hamiltonian in Section 4. Section 5
details the time integration over the exposure of a given experiment. For
interested users, Section 6 provides a quick start guide to PEANUTS and a
explanations on how to reproduce our validation procedure. Lastly, we provide
our conclusions and outlook in Section 7.

##  2 Solar neutrino flux

###  2.1 Neutrino survival probability at Sun surface

Solar neutrinos are produced over a wide region within the Sun, making the
incertitude over the production point much larger than the typical detectors
size (in fact, much larger than the Earth itself. See e.g. Fig. 1) [37]. This
feature implies that the solar neutrino flux at Earth is given by an
_incoherent_ superposition of neutrino mass eigenstates [38, 39], whose
composition remains constant as the flux propagates in vacuum333The
oscillation dynamics would be different in case of a coherent flux [40].. If
the neutrino oscillation length-scale is much smaller than the scale over
which matter density varies significantly, neutrino oscillations in matter
proceed adiabatically [41]; the values of the neutrino oscillation parameters
inferred by global fits of neutrino data [18, 19, 20] imply that the adiabatic
regime is realised for solar neutrinos as they propagate from the production
point towards the Sun surface. In the adiabatic approximation, the neutrino
flavour composition at the Sun surface only depends on the matter effects at
neutrino production point.

Given $U=U(\theta_{12},\theta_{13},\theta_{23},\delta)$ the PMNS mixing matrix
in vacuum (cf. eq.s (15, 16)), it is possible to define an analogous matrix
$T$ that diagonalises the neutrino Hamiltonian in matter by simply replacing
the vacuum values of $\theta_{12},\theta_{13}$ with their matter-rotated ones
[42, 43]444See also [44, 45].:

$\displaystyle\cos 2\widetilde{\theta}_{13}$ $\displaystyle=$
$\displaystyle\frac{(\cos 2\theta_{13}-A/\Delta m^{2}_{ee})}{\sqrt{(\cos
2\theta_{13}-A/\Delta m^{2}_{ee})^{2}+\sin^{2}2\theta_{13}~{}}},$ (1)
$\displaystyle\cos 2\widetilde{\theta}_{12}$ $\displaystyle=$
$\displaystyle\frac{(\cos 2\theta_{12}-A^{\,\prime}/\Delta
m^{2}_{21})}{\sqrt{(\cos 2\theta_{12}-A^{\,\prime}/\Delta
m^{2}_{21})^{2}~{}+~{}\sin^{2}2\theta_{12}\cos^{2}(\widetilde{\theta}_{13}-\theta_{13})~{}~{}}},$
(2)

where $\theta_{12},\theta_{13}$ are the PMNS mixing angles in vacuum,
$\widetilde{\theta}_{12},\widetilde{\theta}_{13}$ the corresponding ones in
matter,

$A=2EV=2\sqrt{2}EG_{F}n_{e}$ (3)

represents the matter potential for a neutrino with energy $E$ travelling in a
medium with electron density $n_{e}$, $G_{F}$ is the Fermi constant,

$\displaystyle A^{\,\prime}\equiv A\,\cos^{2}\widetilde{\theta}_{13}+\Delta
m^{2}_{ee}\sin^{2}(\widetilde{\theta}_{13}-\theta_{13})$ (4)

is the $\theta_{13}$-modified matter potential and

$\displaystyle\Delta m^{2}_{ee}\equiv\cos^{2}\theta_{12}\Delta
m^{2}_{31}+\sin^{2}\theta_{12}\Delta m^{2}_{32}.$ (5)

The matrix $T$ is then simply defined as

$\displaystyle
T\left(E,n_{e}\right)=U(\widetilde{\theta}_{12},\widetilde{\theta}_{13},\theta_{23},\delta).$
(6)

Notice that the matrix $T$ depends both on local matter density and energy of
the produced neutrino.

Numerically (cf. eq.s 4.17, 4.18 in [46])

$\displaystyle V$ $\displaystyle=$
$\displaystyle\sqrt{2}G_{F}n_{e}=\frac{3.868\times
10^{-7}}{m}\times\frac{n_{e}}{\text{mol}/\text{cm}^{3}},$ (7) $\displaystyle
k$ $\displaystyle=$ $\displaystyle\frac{\Delta
m^{2}}{2E}=\frac{2.533}{m}\times\frac{\Delta
m^{2}}{\text{eV}^{2}}\times\frac{\text{MeV}}{E},$ (8)

where the wavenumber $k$ is useful in the relation

$\frac{A}{\Delta m^{2}}=\frac{2\sqrt{2}EG_{F}n_{e}}{\Delta
m^{2}}=\frac{V}{k}.$ (9)

The probability of producing a neutrino mass eigenstate $i$ in matter is thus
$\left|T_{\alpha i}\left(E,n_{e}\right)\right|^{2}$, where $\alpha$ is the
flavour of the charged lepton entering the neutrino production vertex.

PEANUTS assumes the validity of the adiabatic regime when computing the solar
neutrino flux, as this guarantees a fast computation and is an excellent
approximation for the neutrino oscillation parameters realised in the standard
3-flavour mixing scheme [47, 48, 49, 50, 35]. Hence, it provides a function
for the computation of the matrix $T$ in eq.(6), with signature

⬇

Tei(th12,th13,DeltamSq21,DeltamSq3l,E,ne)

where the arguments correspond, respectively, to $\theta_{12}$, $\theta_{13}$,
$\Delta m_{21}^{2}$ (eV2), $\Delta m_{3\ell}^{2}$ (eV2), $E$ (MeV) and $n_{e}$
(mol/cm3). Note that in PEANUTS the input variable $\Delta m_{3\ell}$ has
different meanings according to the ordering of neutrino mass eigenstates. In
normal ordering (NO) $l=1$, and thus $\Delta m_{31}^{2}=\Delta m_{3\ell}^{2}$
and $\Delta m_{32}^{2}=\Delta m_{3\ell}^{2}-\Delta m_{21}^{2}$, whereas for
inverted ordering (IO) $l=2$, so $\Delta m_{32}^{2}=\Delta m_{3\ell}^{2}$ and
$\Delta m_{31}^{2}=\Delta m_{3\ell}^{2}+\Delta m_{21}^{2}$. PEANUTS also
provides accessible functions for the mixing angles in matter
$\widetilde{\theta}_{12}$ and $\widetilde{\theta}_{13}$, as well as useful
quantities such as $\Delta m_{ee}^{2}$ and the ratio $V/k$, in the form of the
functions

⬇

th12_M(th12,th13,DeltamSq21,DeltamSq3l,E,ne)

th13_M(th12,th13,DeltamSq21,DeltamSq31,E,ne)

DeltamSqee(th12,DeltamSq21,DeltamSq3l)

Vk(Deltam2,E,ne)

In the adiabatic approximation neutrinos evolve as pure mass eigenstates
within the Sun. For a fixed value of neutrino energy $E$, the flux composition
at Sun surface is given by the average over the neutrino production points
inside the Sun. Assuming spherical symmetry, if $f(r)$ is the fraction of
neutrinos produced at point $r\equiv R/R_{\odot}$, where $R_{\odot}$ is the
solar radius and $R$ the distance from the center of the Sun, the probability
of a solar neutrino with energy $E$ to emerge as mass eigenstate $i$ is

$P_{\nu_{e}\rightarrow\nu_{i}}^{\odot}(E)=\int_{0}^{1}\text{d}r\left|T_{ei}\left(E,n_{e}(r)\right)\right|^{2}f(r),$
(10)

with the normalization

$\int_{0}^{1}\text{d}rf(r)=1.$ (11)

Figure 1: Solar density (left) and fraction of solar neutrinos produced
(right), as a function of relative radius $r\equiv R/R_{\odot}$, from Solar
model BS05(AGS,OP) [2].

The computation of $f(r)$ assumes a specific solar model. In our validation of
PEANUTS we use the BS05(AGS,OP) model [2]555Numerical data are available at
http://www.sns.ias.edu/~jnb/, which is one of the models assumed by the SNO
collaboration in their neutrino oscillation fit [35]. Fig. 1-left shows the
value of solar matter density $n_{e}(r)$ as a function of solar radius $r$ in
the BS05(AGS,OP) model, while Fig. 1-right shows the neutrino production
fractions for different production chains. Note that the electron density in
[2] is given in units of
$\frac{1}{\text{cm}^{3}N_{A}}=\frac{\text{mol}}{\text{cm}^{3}\ 6.022\cdot
10^{23}}.$

The choice of a particular solar model is encoded in PEANUTS as the Python
class SolarModel. Through the constructor of the class, of signature,

⬇

SolarModel(solar_model_file=None, flux_file=None, spectrum_files=None,
fluxrows=None, fluxcols=None, fluxscale=None, distrow=None, radiuscol=None,
densitycol=None, fractioncols=None)

one can select the location of a file describing the chosen solar model with
the optional argument solar_model_file. By default, when no other file
location is provided, PEANUTS will assume the use of the BS16(AGSS09) model
[51]. At the time of publication, PEANUTS can also work out of the box with
the B16 (GS98) [51] model, the BS05(AGS,OP) and BS05(OP) models [2], and the
BP00 [37] solar model. The user is nevertheless encouraged to implement their
own solar model (and thus neutrino fraction distribution $f(r)$), by providing
a custom solar model file via the solar_model_file optional argument and/or
flux file with flux_file (they will be assumed to be the same if the latter
argument is missing), but in such case one must also specify the rows and
columns where the relevant information in the files can be located. The
options fluxrows and fluxcols sets the rows and columns in flux_file where one
can find the total neutrino flux per fraction, which can be either
dictionaries or real numbers, but at least one of the two must be a dictionary
with the names of the fractions of interest and the corresponding row or
column. If the fluxes must be rescaled, the option fluxscale can be provided,
which can be a real number, for blanket rescaling of all fluxes, or a
dictionary with different rescaling for each fraction. The option distrow
points to the first row of the table in solar_model_file containing the
distributions, i.e. the radius, density and fraction samples, whose columns
must also be specified with the options radiuscol, densitycol and
fractioncols, respectively, with the latter a dictionary of names and columns
of the fractions of interest. Naturally all of these optional arguments must
be provided if the solar model file is not known to PEANUTS.

In addition to allowing the selection of different solar models, PEANUTS also
allows the use of different energy spectra for the various neutrino fractions.
By default the ${}^{8}B$ and $hep$ spectra are assumed and taken from [52] and
[53] respectively. Spectra for other neutrino fractions or different spectra
can be provided via the optional argument spectrum_files in the constructor
for the SolarModel class. It should be noted that spectra are assumed in
PEANUTS to be normalised to 1, so the user should make sure to renormalise the
spectrum accordingly, as it is the case, for instance for the ${}^{8}B$
spectrum from [54], which is normalized to 1000.

In addition to the neutrino fraction distributions and fluxes, the SolarModel
class in PEANUTS provides the density of the Sun at various radius samples,
necessary for computing the probability of neutrinos at Sun exit from eq.
(10). This probability, or equivalenty the weight of neutrino mass eigenstates
in the surface of the Sun, is computed by PEANUTS with the function

⬇

solar_flux_mass(th12, th13, DeltamSq21, DeltamSq3l, E, radius_samples,
density, fraction)

###  2.2 Neutrino propagation from the Sun

Being it an incoherent flux, the fraction of mass eigenstates within the solar
neutrino flux remains constant as long as neutrinos propagate in the vacuum,
on their path from Sun surface to Earth. Here we adopt the convention that the
unitary PMNS mixing matrix $U$ defines the change of basis between mass and
flavour neutrino _fields_

$\nu_{\alpha}(x)=U_{\alpha i}\ \nu_{i}(x),$ (12)

with $\alpha=e,\mu,\tau$ the flavour indexes, $i=1,2,3$ the mass ones. With
this convention, a neutrino _state_ produced at the origin results from the
linear superposition of mass states666See [46] for a comprehensive discussion
on the relations between fields, plane waves and momentum eigenstates.

$\left|\nu_{\alpha},x=0\right>=U_{\alpha i}^{*}\left|\nu_{i},x=0\right>,$ (13)

implying that the probability of observing a neutrino of flavour $\alpha$ from
a mass eigenstate $i$ is given by $|U_{\alpha i}^{*}|^{2}$. Thus, the
probability for a solar neutrino to manifest as flavour $\alpha$ is given by

$P^{S}_{\nu_{e}\rightarrow\nu_{\alpha}}(E)=|U_{\alpha
i}^{*}|^{2}P^{\odot}_{\nu_{e}\rightarrow\nu_{i}}(E),$ (14)

and we assume throughout this paper that repeated indexes are summed.

In PEANUTS the PMNS matrix is implemented as the PMNS class, which is
constructed from the mixing angles $\theta_{ij}$ and CP phase
$\delta_{\textrm{CP}}$ as

⬇

PMNS(th12,th13,th23,d)

This PMNS provides accessor functions for all mixings parameters, as well as
other useful quantities such as the orthogonal/unitary matrices $R_{12}$,
$R_{13}$, $R_{23}$ and $\Delta$, which allow the PMNS mixing matrix $U$ to be
expressed as

$U=R_{23}\Delta R_{13}\Delta^{*}R_{12},$ (15)

and are defined as

$\begin{array}[]{ll}R_{23}=\left(\begin{array}[]{ccc}1&0&0\\\
0&c_{23}&s_{23}\\\
0&-s_{23}&c_{23}\end{array}\right),&R_{13}=\left(\begin{array}[]{ccc}c_{13}&0&s_{13}\\\
0&1&0\\\ -s_{13}&0&c_{13}\end{array}\right),\\\
R_{12}=\left(\begin{array}[]{ccc}c_{12}&s_{12}&0\\\ -s_{12}&c_{12}&0\\\
0&0&1\end{array}\right),&\Delta=\left(\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
0&0&e^{i\delta}\end{array}\right).\end{array}$ (16)

This decomposition of the PMNS matrix will be useful further down when
computing the propagation Hamiltonian through Earth. Finally, in order to
compute the solar probability for flavour eigenstates with PEANUTS, one can
use the function

⬇

PSolar(pmns,DeltamSq21,DeltamSq3l,E,radius_samples,density,fraction)

which simply implements eq. (14) using a PMNS object and calling the
solar_flux_mass function, and returns a list of the probability for each
flavour eigenstate.

##  3 Neutrino propagation through Earth

###  3.1 Probability of transition through matter

If the neutrino flux from the Sun crosses the Earth (or any finite density
matter in general) the probabilities are modified, since the propagation
eigenstates in matter differ from the vacuum ones. In general the flavour
eigenstate at time777Since we always assume ultrarelativistic neutrinos, we
can interchangeably use traveled distance $x$ or elapsed time $t$ to identify
the evolved neutrino state. $t$ is obtained by evolving the one at $t=0$ with
the appropriate evolutor operator $\mathcal{U}$

$\left|\nu_{\alpha},t\right>=\mathcal{U}_{\alpha\beta}(t)\left|\nu_{\beta},0\right>,$
(17)

where we assume a pure flavour (vacuum) eigenstate at $t=0$,

$\left|\nu_{\alpha},0\right>=U_{\alpha i}^{*}\left|\nu_{i},0\right>.$ (18)

Notice that the propagation of a single neutrino within Earth is given by the
_coherent_ propagation of incoming mass eigenstates. The determination of the
evolutor operator $\mathcal{U}$ is in general a non-trivial problem, and will
be discussed in detail in the following sections; for the moment let us assume
we know a closed form expression for it.

A general neutrino state $\left|\nu\right>$ (not necessarily an eigenstate)
can always be expressed as a linear superposition of flavour or mass
eigenstates. Working in the flavour basis

$\left|\nu,t=0\right>=c_{\alpha}(t=0)\left|\nu_{\alpha},t=0\right>,$ (19)

where $\left|\nu_{\alpha},0\right>$ stands for a vacuum flavour eigenstate.
This state will evolve as

$\left|\nu,t\right>=c_{\alpha}(0)\
\mathcal{U}_{\alpha\beta}(t)\left|\nu_{\beta},0\right>,$ (20)

which allows to define the linear coefficients at time $t$

$c_{\alpha}(t)=\left<\nu_{\alpha},0|\nu,t\right>=\mathcal{U}_{\alpha\beta}^{T}(t)\
c_{\beta}(0),$ (21)

coinciding with the probability amplitude of the state $\left|\nu,t\right>$ to
interact as flavour $\alpha$ at time $t$.

A mass eigenstate expressed as linear combination of flavour eigenstates is

$\left|\nu_{i},0\right>=U_{i\alpha}^{T}\left|\nu_{\alpha},0\right>=U_{\alpha
i}\left|\nu_{\alpha},0\right>,$ (22)

which implies a transition amplitude from (evolved) mass to flavour eigenstate

$\left<\nu_{\alpha},0|\nu_{i},t\right>=\mathcal{U}_{\alpha\beta}^{T}(t)U_{\beta
i}.$ (23)

Putting everything together, the final probability for a solar neutrino to
manifest as $\alpha$ flavour while crossing the Earth is given by

$P_{\alpha}^{SE}(t,E)=\left|\mathcal{U}_{\alpha\beta}^{T}(t)U_{\beta
i}\right|^{2}P_{\nu_{e}\rightarrow\nu_{i}}^{\odot}(E),$ (24)

where $t=0$ is defined at the time of neutrino crossing the Earth surface. The
interpretation of Eq. (24) is the following: $U_{\beta i}$ are the
coefficients of the mass eigenstate $i$ expressed as linear combination of
flavour eigenstates, $\left|\nu_{i},0\right>=U_{\beta
i}\left|\nu_{\beta},0\right>$, and $\mathcal{U}^{T}_{\alpha\beta}(t)U_{\beta
i}=\left<\nu_{\alpha},0|\nu_{i},t\right>$ is the transition amplitude from the
evolved mass eigenstate $i$ to interaction eigenstate $\alpha$. Finally, each
probability $|\left<\nu_{\alpha},0|\nu_{i},t\right>|^{2}$ is multiplied by the
weight of the mass eigenstate $i$ in the incoherent solar flux,
$P_{\nu_{e}\rightarrow\nu_{i}}^{\odot}(E)$.

The probability of oscillation for each flavour eigenstate in eq. (24) is
implemented in PEANUTS by the function Pearth, with signature

⬇

Pearth(nustate, density, pmns, DeltamSq21, DeltamSq3l, E, eta, depth,
mode=”analytical”,

basis=”mass”, full_oscillation=False, antinu=False)

which takes as arguments the neutrino state, nustate, an instance of the Earth
density class (see Section 3.3 below), density, an instance of the PMNS class,
pmns, the mass splitting parameters, DeltamSq21 and DeltamSq3l, the neutrino
energy, E, the nadir angle of the incoming neutrino, eta, and the depth of the
experiment at which the probability is to be computed, depth.

The optional argument basis defines the basis of the neutrino state. If
basis=”flavour” the initial state is assumed to be a _coherent_ one, expressed
in flavour basis, with nustate defining the complex coefficients $c_{\alpha}$
in eq. (19), and final probabilities computed by squaring the coefficients
evolved as in eq. (21). If basis=”mass”, the initial state is assumed to be an
_incoherent_ superposition of mass eigenstates, with nustate defining the real
weights $P_{\nu_{e}\rightarrow\nu_{i}}^{\odot}(E)$ and final probabilities
computed as in eq. (24). To compute the probability for an antineutrino, one
can set the optional argument antinu=True (False by default).

One can also optionally select the evolution mode to be either numerical or
analytical (default) by providing the optional argument mode with either
option. The function Pearth thus splits into two functions for each of the
methods, Pearth_numerical and Pearth_analytical, whose details and differences
will be described below in Section 4. Lastly, one can request the full
evolution of the probability with the optional argument full_oscillation
(defaults to False), returned as a list of probability values for each flavour
eigenstate and certain discrete coordinate locations along the path of the
neutrino. Note that the full oscillation can only be provided via the
numerical evolution mode, so if selected along with the analytical mode,
PEANUTS will provide a warning and simply compute the final probabilities.

In addition, if one wishes to know the final evolved complex coefficients
$c_{\alpha}(t)$ from a coherent (flavour basis) initial neutrino state, from
eq. 21, they can be obtained with the function

⬇

evolved_state(nustate, density, pmns, DeltamSq21, DeltamSq3l, E, eta, depth,
mode=”analytical”, full_oscillation=False, antinu=False)

which has the same arguments as Pearth, except the basis as this is only
available for neutrino states in the flavour basis. As with the probability
function, this function splits into evolved_state_numerical and
evolved_state_analytical according to the selected mode.

###  3.2 Evolutor operator

The evolutor operator $\mathcal{U}$ is defined by the equations

$\left|\nu_{\alpha},t\right>=\mathcal{U}_{\alpha\beta}(t)\left|\nu_{\beta},0\right>,\hskip
28.45274pt\text{ with }\hskip 28.45274pt\mathcal{U}(0)=\mathbb{1}.$ (25)

It is possible to define an effective Hamiltonian $H$ for neutrinos [46], such
that the flavour states obey to the Schrödinger equation

$\frac{\text{d}}{\text{d}t}\left|\nu_{\alpha},t\right>=-iH_{\alpha\beta}(t)\left|\nu_{\beta},t\right>,$
(26)

thus giving the differential equation for $\mathcal{U}$

$\displaystyle\frac{d}{\text{d}t}\mathcal{U}_{\alpha\beta}(t)\left|\nu_{\beta},0\right>$
$\displaystyle=$ $\displaystyle-i\ H_{\alpha\gamma}(t)\
\mathcal{U}_{\gamma\beta}(t)\left|\nu_{\beta},0\right>$ (27)
$\displaystyle\Rightarrow\hskip
5.69046pt\frac{\text{d}}{\text{d}t}\mathcal{U}(t)$ $\displaystyle=$
$\displaystyle-i\ H(t)\ \mathcal{U}(t),$ (28)

where the second equation follows from the states
$\left|\nu_{\alpha},0\right>$ forming a complete orthonormal basis.

Eq. (28) allows to express the evolutor as a function of the Hamiltonian of
the system. Its formal solution is

$\mathcal{U}(t,t_{0})=\mathcal{T}\left[e^{-i\int_{t_{0}}^{t}\text{d}t^{\prime}H(t^{\prime})}\right],$
(29)

where $\mathcal{T}$ is the time-order operator. Eq. (29) does not generally
admit an analytic closed form, except for very special cases, for instance if
the Hamiltonian at different times does commute.

A well known approach to the problem is the Dyson series [55]

$\mathcal{U}(t,t_{0})=\mathbb{1}+\sum_{n=1}^{\infty}\frac{\left(-i\right)^{n}}{n!}\int_{t_{0}}^{t}\text{d}t_{1}\int_{t_{0}}^{t}\text{d}t_{2}\cdots\int_{t_{0}}^{t}\text{d}t_{n}\mathcal{T}\left[H(t_{1})H(t_{2})\cdots
H(t_{n})\right],$ (30)

which allows for an approximate solution obtained by truncating Eq. (30) at
finite values of $n$, if the series is expected to be perturbative. We will
return to this approximation in Section 4, where we will derive an
approximated analytical expression for the evolutor.

###  3.3 Earth matter regeneration

Due to the incoherent nature of the solar neutrino flux, its relative weights
in mass eigenstates components, $P_{\nu_{e}\rightarrow\nu_{i}}^{\odot}(E)$,
are constant while travelling from Sun to Earth in vacuum. However, before
arriving at a terrestrial detector, a solar neutrino can cross the Earth
itself, along which path the matter potential makes the propagation
eigenstates different from the vacuum ones. This results in coherent neutrino
oscillations inside the Earth that, on average, result in a regeneration of
$\nu_{e}$ with respect to the vacuum case [38].

The electron density inside Earth can be parametrised by 5 shells, within
which the density itself varies smoothly as [39]

$N_{j}(r)=\alpha_{j}+\beta_{j}r^{2}+\gamma_{j}r^{4},\hskip
28.45274pt\text{with}\hskip
28.45274pt\left[N\right]=\text{mol}/\text{cm}^{3},$ (31)

and where $r$ is the radial distance normalised to the Earth radius. The
numerical values of the parameters are reported for convenience in Table 1.

$j$ | Shell | $[r_{j-1},\,r_{j}]$ | $\alpha_{j}$ | $\beta_{j}$ | $\gamma_{j}$
---|---|---|---|---|---
1 | Inner core | $[0,\,0.192]$ | 6.099 |  $-$4.119 | 0.000
2 | Outer core | $[0.192,\,0.546]$ | 5.803 |  $-$3.653 |  $-$1.086
3 | Lower mantle | $[0.546,\,0.895]$ | 3.156 |  $-$1.459 | 0.280
4 | Transition Zone | $[0.895,\,0.937]$ |  $-$5.376 | 19.210 |  $-$12.520
5 | Upper mantle | $[0.937,\,1]$ | 11.540 |  $-$20.280 | 10.410
Table 1: Values of the parameters for the electron density expressed as
$N_{j}(r)=\alpha_{j}+\beta_{j}r^{2}+\gamma_{j}r^{4}$ with
$[N]=\text{mol}/\text{cm}^{3}$, for each of the Earth internal shells, as
derived in [39]. The radial distance $r$ is normalised to the radius of Earth.
Figure 2: Earth section showing the different shells and the path of a solar
neutrino having nadir angle $\eta$. The trajectory coordinate $x$ is also
schematised.

The parametrisation in Eq. (31) is valid for radial trajectories, i.e. paths
crossing the center of the Earth. For a path forming a nadir angle $\eta$ with
the radial trajectory (cf. Fig. 2), the parametrisation is functionally
invariant, with modified coefficients

$\displaystyle N_{j}(x)$ $\displaystyle=$
$\displaystyle\alpha^{\prime}_{j}+\beta^{\prime}_{j}x^{2}+\gamma^{\prime}_{j}x^{4},$
(32) $\displaystyle\alpha^{\prime}_{j}$ $\displaystyle=$
$\displaystyle\alpha_{j}+\beta_{j}\sin^{2}\eta+\gamma_{j}\sin^{4}\eta,$ (33)
$\displaystyle\beta^{\prime}_{j}$ $\displaystyle=$
$\displaystyle\beta_{j}+2\gamma_{j}\sin^{2}\eta,$ (34)
$\displaystyle\gamma^{\prime}_{j}$ $\displaystyle=$ $\displaystyle\gamma_{j},$
(35)

where the _trajectory coordinate_ $x$ is defined as the distance from the
trajectory mid-point.

The Earth density profiles for some example values of the nadir angle $\eta$
are reported in Fig. 3.

Figure 3: Earth density profile for different values of the nadir angle
$\eta$, following the parametrisation in [39].

This parametrisation assumes a nadir angle defined for a detection at the
surface of Earth. If the detector is located underground, we can define a
“detector shell” at the detector radial distance. For instance, SNO was placed
$H=2$ km underground, so we can define
$r_{det}=r_{\text{SNO}}=1-H/R_{\text{E}}\equiv 1-h$, where $R_{\text{E}}$ is
the (not rescaled) Earth radius, $R_{E}=6.371\cdot 10^{3}$ km. There are two
modifications for the scenario of an underground detector with respect to the
treatment in [39]: the first is that the measured nadir angle $\eta$ differs
from the angle $\eta^{\prime}=\arcsin(r_{det}\sin\eta)$ that one would measure
at the Earth surface, for same neutrino trajectory, cf. Figure (4). This
implies that it is the $\eta^{\prime}$ angle that should be used in Eq. (32)
to compute the value of the electron density profile along the neutrino
trajectory, and not the value $\eta$ measured by the experiment. The second
modification is that matter effects are present even for values of
$\eta\geq\pi/2$, if the detector is underground. The contribution to the
trajectory from the outer layer (between Earth surface and detector shell) is
given by

$\Delta x(\eta)=\left\\{\begin{array}[]{l c
c}-r_{det}\cos\eta+\sqrt{1-r_{det}^{2}\sin^{2}\eta}&\text{for}&0\leq\eta\leq\frac{\pi}{2},\\\
r_{det}\cos\eta+\sqrt{1-r_{det}^{2}\sin^{2}\eta}&\text{for}&\frac{\pi}{2}\leq\eta\leq\pi.\end{array}\right.$
(36)

$\Delta x$ reduces to $h$ for $\eta=0,\pi$, and attains the maximum value
$\sqrt{h(2-h)}$ for $\eta=\pi/2$ (this is approximately 160 km for SNO).

Figure 4: Schematic representation of Earth shell and “detector shell” for an
underground detector.

In general, not every shell is crossed by solar neutrinos: a shell $i$ is
crossed by a neutrino trajectory with nadir angle $\eta$ if
$r_{i}>r_{det}\sin\eta$. The value of the trajectory coordinate at each shell
crossing is given by

$x_{i}=\sqrt{r_{i}^{2}-r_{det}^{2}\sin^{2}\eta},\hskip 28.45274pt\text{ for
$i$ such that}\hskip 28.45274ptr_{i}>r_{det}\sin\eta.$ (37)

PEANUTS implements the Earth density as a Python class called EarthDensity,
with signature

⬇

EarthDensity(density_file=None)

where the optional argument density_file allows the user to supply a custom
Earth density profile, provided it is given as table with columns
$\\{r_{j},\alpha_{j},\beta_{j},\gamma_{j}\\}$; notice that PEANUTS can work
with an arbitrary number of shells, however the density profile within each
shell is assumed to follow the functional relation in eq. (31). This class
provides methods to compute the modified coefficients $\alpha^{\prime}_{j}$,
$\beta^{\prime}_{j}$ and $\gamma^{\prime}_{j}$, the radii and nadir angles
corresponding to each Earth shell, as well as the value of the density at
given coordinate $x$ and nadir angle $\eta$ using eq. (32).

##  4 Neutrino propagation Hamiltonian

The propagation Hamiltonian for an ultrarelativistic neutrino propagating in a
medium with electron density $n_{e}(x)$ is, in the flavour basis888We follow
in this work the notation and conventions of [46].

$H_{\nu}=\underbrace{U^{*}\text{diag}(k)U^{T}}_{H_{\nu}^{0}}+\underbrace{V(x)\
\text{diag}(1,0,0)}_{V_{\nu}},$ (38)

with

$\displaystyle k_{i}$ $\displaystyle=$ $\displaystyle\frac{m_{i}^{2}}{2E},$
(39) $\displaystyle V(x)$ $\displaystyle=$
$\displaystyle\sqrt{2}G_{F}n_{e}(x).$ (40)

For antineutrinos, the same Eq. (38) holds, with the replacements

$\displaystyle U$ $\displaystyle\rightarrow U^{*},$ (41) $\displaystyle V$
$\displaystyle\rightarrow-V.$ (42)

To streamline computations and statistical analysis, it is more convenient to
redefine the Hamiltonian by subtracting a constant term

$H_{\nu}\rightarrow H_{\nu}-U^{*}\left(k_{j}\mathbb{1}\right)U^{T},$ (43)

with $j=1$ for normal ordering (NO) and $j=2$ for inverted ordering (IO), such
that

$\displaystyle\left(\begin{array}[]{ccc}k_{1}&&\\\ &k_{2}&\\\
&&k_{3}\end{array}\right)-k_{1}\mathbb{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2E}\left(\begin{array}[]{ccc}0&&\\\ &\Delta
m_{21}^{2}&\\\ &&\Delta m_{31}^{2}\end{array}\right)\hskip 5.69046pt\text{ for
NO},$ (50) $\displaystyle\left(\begin{array}[]{ccc}k_{1}&&\\\ &k_{2}&\\\
&&k_{3}\end{array}\right)-k_{2}\mathbb{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2E}\left(\begin{array}[]{ccc}-\Delta m_{21}^{2}&&\\\
&0&\\\ &&\Delta m_{32}^{2}\end{array}\right)\hskip 5.69046pt\text{ for IO},$
(57)

so that we can use as free parameters $\Delta m_{21}^{2}>0$ and $\Delta
m_{3\ell}^{2}$, with $\ell=1$ for NO ($\Delta m_{31}^{2}>0$) and $\ell=2$ for
IO ($\Delta m_{32}^{2}<0$). In the following we keep using the notation
$\text{diag}(k)$ for the general expressions, valid for any choice of mass
ordering. A specific scenario can then be easily recovered by specifying the
structure of this diagonal matrix, e.g. in Eq.s (50, 57).

With the parametrisation in Eq. (15) and using
$\left[\Delta,R_{12}\right]=\left[V_{\nu},R_{23}\right]=\left[V_{\nu},\Delta\right]=0$,
the propagation Hamiltonian can be rewritten as

$H_{\nu}=R_{23}\Delta^{*}\tilde{H}\Delta R_{23}^{T},$ (58)

with

$\tilde{H}=R_{13}R_{12}\text{diag}(k)R_{12}^{T}R_{13}^{T}+V(x)\text{diag}(1,0,0).$
(59)

Notice that $\tilde{H}$ does not depend on $\theta_{23}$, $\delta$.

Given that $R_{23}$ and $\Delta$ do not depend on position, they can be
factorised in the time-ordered definition of the evolutor operator

$\displaystyle\mathcal{U}$ $\displaystyle=$
$\displaystyle\mathcal{T}e^{-i\int\text{d}xH_{\nu}(x)}=\mathbb{1}+R_{23}\Delta^{*}\mathcal{T}\left[-i\int\text{d}x\tilde{H}_{\nu}(x)\right]\Delta
R_{23}^{T}+$ (62)
$\displaystyle\frac{1}{2}R_{23}\Delta^{*}\mathcal{T}\left[\left(\left(-i\right)^{2}\int\text{d}x_{1}\text{d}x_{2}\tilde{H}_{\nu}(x_{1})\underbrace{\Delta
R_{23}^{T}R_{23}\Delta^{*}}_{=\mathbb{1}}\tilde{H}_{\nu}(x_{2})\right)\right]\Delta
R_{23}^{T}$ $\displaystyle+\dots(\text{ other terms of the Dyson series for
}n\rightarrow\infty)$ $\displaystyle=$ $\displaystyle
R_{23}\Delta^{*}\mathcal{T}\left[e^{-i\int\text{d}x\tilde{H}(x)}\right]\Delta
R_{23}^{T}=R_{23}\Delta^{*}\tilde{\mathcal{U}}\Delta R_{23}^{T}.$ (63)

A numerical solution for the evolutor can be obtained by solving eq. (63) or
by resolving the differential equation in eq. (28). PEANUTS offers a numerical
evaluation of the evolutor, via the function

⬇

Pearth_numerical(nustate, density, pmns, DeltamSq21, DeltamSq3l, E, eta,
depth, basis=”mass”, full_oscillation=False, antinu=False)

and the evolved coefficients from a coherent neutrino state, with

⬇

evolved_state_numerical(nustate, density, pmns, DeltamSq21, DeltamSq3l, E,
eta, depth, full_oscillation=False, antinu=False)

These computations, however, can be extremely time-consuming, and thus it is
convenient to find an approximated analytical expression by performing a
perturbative expansion of the Hamiltonian.

###  4.1 Perturbative expansion of the neutrino propagation Hamiltonian

We are interested in an expression for the operator in Eq. (63)

$\tilde{\mathcal{U}}=\mathcal{T}\left[e^{-i\int\text{d}l\tilde{H}(l)}\right],$
(64)

where $l$ is the coordinate along the neutrino path. We normalise distances to
the Earth radius $R_{E}$, by defining $x=\l/R_{E}$. The Hamiltonian
$\tilde{H}$ can be divided in a kinetic and a matter dependent terms

$\tilde{H}(x)=\tilde{H}_{k}+\sqrt{2}G_{F}n_{e}(x)\text{diag}(1,0,0),$ (65)

where $\tilde{H}_{k}$ does not depend on $x$.

To work out a perturbative expression for $\mathcal{U}$999For now we drop the
tilde from the “reduced” evolutor $\tilde{\mathcal{U}}$ for covenience, and to
avoid confusion with the average $\bar{\mathcal{U}}$. We will recover the
notation at the end of the section when computing the full evolutor. it is
convenient to express the electron density as a perturbation along its mean
value along the path [39]

$n_{e}(x)=\bar{n}_{e}+\delta n(x),\hskip
28.45274pt\bar{n}_{e}=\frac{1}{x_{2}-x_{1}}\int_{x_{1}}^{x_{2}}\text{d}x\
n_{e}(x),$ (66)

from which it follows

$\int_{x_{1}}^{x_{2}}\text{d}x\ \delta n(x)=0.$ (67)

We can analogously divide the Hamiltonian into a zeroth order term and a
perturbation

$\tilde{H}(x)=\underbrace{\tilde{H}_{k}+\sqrt{2}G_{F}\bar{n}_{e}\text{diag}(1,0,0)}_{\tilde{H}_{0}}+\underbrace{\sqrt{2}G_{F}\delta
n(x)\text{diag}(1,0,0)}_{\delta\tilde{H}(x)},$ (68)

where again $\tilde{H}_{0}$ does not depend on $x$.

The evolutor can thus be expressed as [39]

$\mathcal{U}(x_{2},x_{1})=\bar{\mathcal{U}}(x_{2},x_{1})-i\int_{x_{1}}^{x_{2}}\text{d}x\
\bar{\mathcal{U}}(x_{2},x)\ \delta\tilde{H}(x)\
\bar{\mathcal{U}}(x,x_{1})+\mathcal{O}(\delta\tilde{H}^{2}),$ (69)

where $\bar{\mathcal{U}}$ is the evolutor for constant matter density
$\bar{n}_{e}$.

The evolutor for a constant Hamiltonian $\bar{H}$ can generally be expressed
in a closed form [56]

$e^{-i\bar{H}x}=\phi\sum_{a=1}^{3}e^{-ix\lambda_{a}}\frac{1}{3\lambda_{a}^{2}+c_{1}}\left[\left(\lambda_{a}^{2}+c_{1}\right)\mathbb{1}+\lambda_{a}T+T^{2}\right]\equiv\phi\sum_{a=1}^{3}e^{-ix\lambda_{a}}M_{a},$
(70)

where $T=\bar{H}-\text{Tr}(\bar{H})\mathbb{1}/3$ is a traceless matrix and
$\lambda_{a}$ are the roots of the characteristic equation

$\lambda^{3}+c_{1}\lambda+c_{0}=0,$ (71)

with

$\displaystyle c_{1}$ $\displaystyle=$ $\displaystyle
T_{11}T_{22}-T_{12}T_{21}+T_{11}T_{33}-T_{13}T_{31}+T_{22}T_{33}-T_{23}T_{32},$
(72) $\displaystyle c_{0}$ $\displaystyle=$ $\displaystyle-\det T.$ (73)

Finally

$\phi=e^{-ix\frac{\text{Tr}(\bar{H})}{3}}.$ (74)

By noticing that the full dependence on $x$ in $e^{-i\bar{H}x}$ is now
contained within the scalar functions $e^{-i\lambda_{a}x}$, the first order
correction in Eq. (69) can be computed as

$\displaystyle\mathcal{U}^{(1)}(x_{2},x_{1})$ $\displaystyle=$
$\displaystyle-i\int_{x_{1}}^{x_{2}}\text{d}x\ \bar{\mathcal{U}}(x_{2},x)\
\delta\tilde{H}(x)\ \bar{\mathcal{U}}(x,x_{1})$ (75) $\displaystyle=$
$\displaystyle-i\sum_{a,b=1}^{3}\int_{x_{1}}^{x_{2}}\text{d}xe^{-i\tilde{\lambda}_{a}(x_{2}-x)}M_{a}\text{diag}\left(\sqrt{2}G_{F}\delta
n(x),0,0\right)M_{b}e^{-i\tilde{\lambda}_{b}(x-x_{1})}$ (76) $\displaystyle=$
$\displaystyle-i\sum_{a,b=1}^{3}M_{a}\text{diag}\left(\sqrt{2}G_{F}I_{ab}(x_{2},x_{1}),0,0\right)M_{b},$
(77)

where $\tilde{\lambda}_{a}=\lambda_{a}+\text{Tr}(\bar{H})/3$ and we defined

$I_{ab}(x_{2},x_{1})=\int_{x_{1}}^{x_{2}}\text{d}x\
e^{-i\tilde{\lambda}_{a}(x_{2}-x)}\ \delta n(x)\
e^{-i\tilde{\lambda}_{b}(x-x_{1})}.$ (78)

For a path fully contained within one shell we can parametrise

$\delta
n(x)=\tilde{\alpha}^{\prime}+\beta^{\prime}x^{2}+\gamma^{\prime}x^{4},$ (79)

where $\tilde{\alpha}^{\prime}=\alpha^{\prime}-\bar{n}_{e}$, implying that
$I_{ab}(x_{2},x_{1})$ can be expressed analytically in closed form.

Summarising, we can perturbatively expand the evolutor operator as

$\displaystyle\mathcal{U}(x_{2},x_{1})$ $\displaystyle=$
$\displaystyle\mathcal{U}^{(0)}(x_{2},x_{1})+\mathcal{U}^{(1)}(x_{2},x_{1})+\mathcal{O}(\delta\tilde{H}^{2}),$
(80) $\displaystyle\mathcal{U}^{(0)}(x_{2},x_{1})$ $\displaystyle=$
$\displaystyle
e^{-i\bar{H}(x_{2}-x_{1})}=\phi\sum_{a=1}^{3}e^{-i(x_{2}-x_{1})\lambda_{a}}M_{a},$
(81) $\displaystyle\mathcal{U}^{(1)}(x_{2},x_{1})$ $\displaystyle=$
$\displaystyle-i\sum_{a,b=1}^{3}M_{a}\
\text{diag}\left(\sqrt{2}G_{F}I_{ab}(x_{2},x_{1}),0,0\right)\ M_{b}.$ (82)

In PEANUTS we implement these perturbative expressions to compute
$\mathcal{U}$ in Eq. (64) at first order in perturbation theory, using a
“reduced” mixing matrix $\tilde{U}=R_{13}R_{12}$, with he function

⬇

Upert(DeltamSq21, DeltamSq3l, E, x2, x1, a, b, c, antiNu)

which depends on the neutrino mass differences squared DeltamSq21 and
DeltamSq3l, the PMNS matrix pmns, neutrino energy E, the start and end points
of the shell along the neutrino path x2 and x1, as well as the density
parameters of the traversed shell, with a = $\alpha$, b = $\beta$ and c =
$\gamma$. The flag anitNu labels whether the computation is to be done for a
neutrino (antiNu=False) or antineutrino (antiNu=True).

The procedure outlined above allows to express the evolutor at 1st order in
$\delta H$, for a path fully contained within one shell. In general, the full
evolutor on a generic path $(x_{1},x_{2})$ can be expressed as a time-ordered
product of evolutors along the same path

$\mathcal{U}(x_{2},x_{1})=\mathcal{U}(x_{2},x_{i})\mathcal{U}(x_{i},x_{1}),$
(83)

where $x_{i}$ is a generic point $x_{1}<x_{i}<x_{2}$ contained on the original
path. Moreover, it can be shown [39] that

$\mathcal{U}(0,-x)=\mathcal{U}(x,0)^{T}.$ (84)

The consequences of Eq.s (83, 84) are twofold: first, for a path starting at
$x=x_{i}$, with $0\leq x_{i}<x_{1}$, crossing $n$ shells with boundaries at
trajectory coordinate $(0,x_{1},x_{2},\dots,x_{n})$, and ending at the point
$x=x_{f}$, with $x_{n-1}<x_{f}\leq x_{n}$, the full evolutor can be expressed
as

$\mathcal{U}(x_{f},x_{i})=\mathcal{U}(x_{f},x_{n-1})\mathcal{U}(x_{n-1},x_{n-2})\dots\mathcal{U}(x_{2},x_{1})\mathcal{U}(x_{1},x_{i}).$
(85)

Second, for detectors placed at surface, the Earth spherical symmetry implies
that the electron density is symmetric with respect to the trajectory midpoint
at $x=0$; thus we only need to compute the evolutor on one half-path (cf. Fig.
2)

$\mathcal{U}(x_{F},x_{I})=\mathcal{U}(x_{F},0)\mathcal{U}(0,-x_{F})=\mathcal{U}(x_{F},0)\mathcal{U}(x_{F},0)^{T}.$
(86)

Notice that the final evolutor is only a function of $\eta$, since both the
density profile and travelled distance within Earth are a function of the
nadir angle.

If the detector is placed underground at trajectory coordinate
$x_{det}<x_{F}$, we distinguish two cases. For $0\leq\eta<\pi/2$

$\mathcal{U}(x_{det},x_{I})=\mathcal{U}(x_{det},0)\mathcal{U}(0,-x_{F})=\mathcal{U}(x_{det},0)\mathcal{U}(x_{F},0)^{T},\hskip
5.69046pt\left(0\leq\eta<\frac{\pi}{2}\right)$ (87)

where $x_{det}=r_{det}\cos{\eta}$. For $\pi/2\leq\eta\leq\pi$ the electron
density can be approximated to a constant value, since the neutrino path is
never deeper than $H$ and density variations are negligible for realistic
detectors. For instance, taking the SNO reference value $H=2$ km,

$\frac{N(R_{E})-N(R_{E}-H)}{N(R_{E})+N(R_{E}-H)}=10^{-4}.$ (88)

We fix for simplicity the electron density value to the one at Earth surface,
$n_{1}=1.67$ mol/cm3. Having assumed constant density along the path, the
evolutor is simply given by

$\mathcal{U}(\eta)=e^{-iR_{E}\Delta
x(\eta)\left(\tilde{H}_{0}+\text{diag}\left(\sqrt{2}G_{F}n_{1},0,0\right)\right)},\hskip
14.22636pt\left(\frac{\pi}{2}\leq\eta<\pi\right).$ (89)

To obtain the full dynamics and thus the full evolutor $\mathcal{U}$, PEANUTS
computes the time-ordered product of “reduced” evolutors $\tilde{\mathcal{U}}$
from above, and inputs it in Eq. (63) to re-introduce the dependence on
$\theta_{23}$ and $\delta$. The PEANUTS function that implements these two
steps is

⬇

FullEvolutor(density, DeltamSq21, DeltamSq3l, pmns, E, eta, depth, antiNu)

with dependency also on the mass differences squared, PMNS matrix and neutrino
energy, but also on the full Earth density profile density, the nadir angle
eta and experiment depth depth. This full evolutor is finally used to compute
the probability of oscillation at the detector by the function

⬇

Pearth_analytical(nustate, density, pmns, DeltamSq21, DeltamSq3l, E, eta,
depth, basis=”mass”, antinu=False)

and the evolved coefficients with

⬇

evolved_state_analytical(nustate, density, pmns, DeltamSq21, DeltamSq3l, E,
eta, depth, antinu=False)

##  5 Time integration and exposure

Solar neutrino experiments typically collect data over a finite interval of
time. As such, the measured survival probability is averaged over exposure

$\left<P_{E}\right>=\frac{\int_{\tau_{d_{1}}}^{\tau_{d_{2}}}\text{d}\tau_{d}\int_{\tau_{h_{1}}(\tau_{d})}^{\tau_{h_{2}}(\tau_{d})}\text{d}\tau_{h}P_{E}\left(\eta(\tau_{d},\tau_{h})\right)}{\int_{\tau_{d_{1}}}^{\tau_{d_{2}}}\text{d}\tau_{d}\int_{\tau_{h_{1}}(\tau_{d})}^{\tau_{h_{2}}(\tau_{d})}\text{d}\tau_{h}},$
(90)

where $\tau_{d}$ is the daily time, $\tau_{h}$ the hourly time and $\eta$ the
Sun nadir angle at detector location. The integration in Eq. (90) is typically
not the most convenient choice for practical applications; a more effective
option is to transform the double integral into a single one over $\eta$ [39]

$\left<P_{E}\right>=\int_{0}^{\pi}\text{d}\eta W(\eta)P_{E}(\eta),$ (91)

where $W(\eta)$ is a normalized weight function representing the fraction of
time in which the experiment collected data at nadir angle $\eta$. For real
experiments, $W(\eta)$ must be provided by the collaboration, taking into
account the actual times at which the detector collected data or has been
offline. It is nevertheless possible to compute $W(\eta)$ analytically, for
the ideal case of an experiment continuously taking data between days
$\tau_{d_{1}}$ and $\tau_{d_{2}}$ [39]. This is done by changing integration
variables

$\displaystyle\int_{\tau_{d_{1}}}^{\tau_{d_{2}}}\text{d}\tau_{d}\int_{\tau_{h_{1}}(\tau_{d})}^{\tau_{h_{2}}(\tau_{d})}\text{d}\tau_{h}P_{E}\left(\eta(\tau_{d},\tau_{h})\right)$
$\displaystyle=$
$\displaystyle\int_{\tau_{d_{1}}}^{\tau_{d_{2}}}\text{d}\tau_{d}\int_{0}^{\pi}\text{d}\eta\frac{\text{d}\tau_{h}(\tau_{d},\eta)}{\text{d}\eta}P_{E}\left(\eta\right)$
(92) $\displaystyle=\int_{0}^{\pi}\text{d}\eta
P_{E}\left(\eta\right)\int_{\tau_{d_{1}}}^{\tau_{d_{2}}}\text{d}\tau_{d}\frac{\text{d}\tau_{h}(\tau_{d},\eta)}{\text{d}\eta}$
$\displaystyle=$ $\displaystyle\int_{0}^{\pi}\text{d}\eta P_{E}(\eta)W(\eta).$
(93)

By normalising the daily and hourly times to the interval $[0,2\pi]$

$\tau_{d}=\frac{\text{day}}{365}2\pi,\hskip
56.9055pt\tau_{h}=\frac{\text{hour}}{24}2\pi,$ (94)

with $\tau_{d}=0$ at winter solstice and $\tau_{h}=0$ at the middle of the
night, it is possible to express

$\tau_{h}=\arccos\left(\frac{\sin(\lambda)\sin(\delta_{S})+\cos(\eta)}{\cos(\lambda)\cos(\delta_{S})}\right),$
(95)

with $\lambda$ the detector latitude and $\delta_{S}$ the Sun declination,
given by

$\delta_{S}=\arcsin\left(-\sin(i)\cos(\tau_{d})\right),$ (96)

with $i=0.4091$ rad being the Earth inclination. With these definitions it is
possible to perform the integral defining $W(\eta)$ in Eq. (93); it is
convenient to restrict $\tau_{d}$ within the interval $[0,\pi]$ (the
alternative case can be easily derived from this one by using the symmetry of
the orbit) and to change the integration variable from $\tau_{d}$ to
$T=\cos(\tau_{d})$. The resulting indefinite integral is expressed in terms of
elementary functions and of the incomplete elliptic integral of the first
kind; its analytic expression is not particularly illuminating but can be
easily evaluated numerically. Some care must be taken in defining the range of
integration for the definite integral, as this is given by the intersection of
three distinct intervals: i) $T\in[-1,1]$ is the interval where
$T=\cos(\tau_{d})$ is defined, ii)
$T\in[\sin(\lambda-\eta)/\sin(i),\sin(\lambda+\eta)/\sin(i)]$ is the range
where $T$ can take values for fixed values of $\eta,\lambda,i$, iii)
$T\in[\cos(\tau_{d_{2}}),\cos(\tau_{d_{1}})]$ is the observation time. If the
intersection of the three intervals is null then $W(\eta)$ will vanish for
that given combination of $\lambda,\eta$ values.

The exposure function $W(\eta)$ is computed in PEANUTS by the function

⬇

NadirExposure(lam=-1, d1=0, d2=365, ns=1000, normalized=False, from_file=None,
angle=”Nadir”)

which has no required arguments, but either the latitude of the experiment,
lam, or an exposure file from_file, must be provided. It returns tabulated
values of the function $W(\eta)$ for ns samples (default 1000) in nadir angle
$\eta$, assuming an exposure from day d1 to day d2, where d1=0 corresponds to
the northern hemisphere winter solstice. The exposure may be selected to be
normalized with the option normalized (defaults to False).

Under default conditions, the function NadirExposure computes analytically the
ideal exposure function performing the integral in eq. (93), which is
implemented in PEANUTS by the function

⬇

IntegralDay(eta, lam, d1=0, d2=365)

which depends on the nadir angle eta, latitude lam and day interval {d1,d2}.
We plot in Fig. 5 the exposure function $W(\eta)$ for one full year of
exposure for three ideal detectors located at latitudes
$\lambda=0,45^{\circ},89^{\circ}$.

Figure 5: Weight of the nadir angle exposure for an ideal experiment located
at latitude $\lambda$, taking data continuously over a full year. The coloured
regions represent the nadir angles subtending the Earth internal shells as
parametrised in Table 1, for a detector located at the Earth surface.

In a realistic case, the exposure must be provided by the experiment. For this
purpose, one can provide an exposure file via the option from_file of the
function NadirExposure, which shall contain the tabulated values of the
exposure. By default it is assumed that the exposure is tabulated in values of
the nadir angle $\eta$, but in some cases, the experiments provide the
exposure tabulated in either the zenith angle $\theta$ or $\cos\theta$. In
those cases, one must specify which angle is used for the tabulation with the
options angle=”Zenith” or angle=”CosZenith”, respectively. If an exposure file
is provided, the latitude is no longer required, and if the argument lam is
provided it will be ignored. Note that irrespective of the original tabulated
values, the function NadirExposure only returns values tabulated in $\eta$.

Lastly, PEANUTS provides the averaged probability of oscillation at detector
location taking into account the finite time exposure, as in eq. (91), with
the function

⬇

Pearth_integrated(nustate, density, pmns, DeltamSq21, DeltamSq3l, E, depth,
mode=”analytical”, full_oscillation=False, antinu=False, lam=-1, d1=0, d2=365,
ns=1000, normalized=False, from_file=None, angle=”Nadir”, daynight=None)

which has the same arguments as the Pearth function, but without the nadir
angle eta, and adds the same optional arguments as NadirExposure, with the
final addition of an optional argument daynight to select integrating only the
nadir angles corresponding to the night period $\eta<\pi/2$, with
daynight=”night”, or the day period $\eta\geq\pi/2$, with daynight=”day”. Note
that this function requires the input neutrino state to be on the mass basis,
since the integration is performed for an incoming incoherent flux of mass
eigenstates, and thus it is not possible to provide a state on the flavour
basis.

##  6 Quick start guide

PEANUTS is an open-source software written in Python, and as such we expect it
to run seemelessly in every Python 3 environment, irrespective of the
architecture. Nevertheless it has been tested thoroughly on Linux and Mac OS X
systems. PEANUTS is optimised to run fast with minimal I/O operations, in
order to allow effortless interface with other frameworks. Hence, it uses
extensively the pre-compile numpy package as well as just-in-time compilation
from numba. The full list of required and optional packages is as follows

  * -

numpy: fast pre-compiled array operations

  * -

numba: just-in-time compilation of functions and classes

  * -

os: file paths

  * -

copy: shallow and deep copies of objects

  * -

time: timing information

  * -

math, cmath, mpmath: mathematical operations

  * -

scipy: integration and interpolation routines

  * -

pyinterval: define integration intervals

  * -

decimal: output formatting

  * -

pandas: reading of csv files

  * -

pyyaml: reading of yaml files Optional

  * -

pyslha: reading of slha files Optional

There are two operational modes of PEANUTS, which we will call the simple and
expert modes. The simple mode allows PEANUTS to be run just from the command
line, appending the input parameters as arguments to the command. Naturally
this is a more limited mode as only provides a fraction of the
functionalities, but it is a fast way to run PEANUTS just from the command
line. The expert mode uses all of the functionalities in PEANUTS, and thus
requires a configuration file, in YAML format, to be written.

In the simple mode one can run one of two provided scripts run_prob_sun.py and
run_prob_earth.py. The first script computes the probability of neutrinos at
the surface of the Sun, and has a signature

⬇

run_prob_sun.py [options] <energy> <fraction> [ <th12> <th13> <th23> <delta>
<dm21> <dm3l>]

where <energy> and <fracion> are mandatory arguments that refer to the
neutrino energy and fraction respectively. With no options provided, the
remaining arguments are required to populate the PMNS matrix and the mass
splittings. However, if the option -i/–in_slha <slha_file> is provided, a file
location is expected after the flag, corresponding to a SLHA file where the
neutrino parameters are defined, and thus only two arguments are required.
Note that this method only works if the package pyslha is installed. Other
options are -s <solar_model> for a different solar model file, -v/–verbose for
verbose output and -h/–help to print usage information.

The second script run_prob_earth.py computes the probability of neutrinos at a
given experimental location below the Earth’s surface. It is used in the
following way

⬇

run_prob_earth.py [options] -f/-m <state> <energy> <eta> <depth> [<th12>
<th13> <th23> <delta> <dm21> <dm3l>]

where, as before, one can provide the full PMNS parameters and mass splittings
explicitly, or as an SLHA file (with the -i/–in_slha <slha_file> option). In
addition one must provide the initial neutrino state with the option
-f/–flavour <state> in the flavour basis, or -m/–mass <state> in the mass
basis for an incoherent incoming flux, as well as the neutrino energy,
<energy>, the nadir angle, <eta> and the depth of the experiment <depth>. Note
that the neutrino state must be given as a sequence of comma-separated values,
without white spaces in between. Additional options include -d/–density
<density> to provide a different Earth density file, –antinu to perform the
computations for antineutrinos, –analytical or –numerical to select either
analytical or numerical evolution, -v/–verbose for verbose output and -h/–help
to print usage information.

The expert mode allows the user to exploit all the functionalities of PEANUTS.
It requires writing a configuration file in YAML format, hence this mode is
only available if the optional pyyaml module is installed. For a given YAML
file, the expert mode can then be used by using the run_peanuts.py executable
in the following way

⬇

run_peanuts.py -f <my_yaml_file>

Since most of the options are provided in the YAML file, this command only
takes the options -v/–verbose for verbose output and -h/–help to print usage
info. An example YAML file to compute the probability at the surface of the
Sun can be seen below

⬇

Energy: 15

Neutrinos:

dm21: 7.42e-05

dm3l: 2.51e-03

theta12: 5.83638e-01

theta23: 8.5521e-01

theta13: 1.49575e-01

delta: 3.40339

Solar:

fraction: ”hep”

flux: true

spectrum: ”distorted”

The YAML file must contain the Energy and Neutrinos nodes, as well as either
the Solar or Earth nodes, or both of them. The Energy node must simply contain
a real number representing the neutrino energy in MeV. The Neutrinos node can
either contain a map of the neutrino parameters, as in the example above, or a
single string with the location of a SLHA file (provided pyslha is installed),
always relative to the location of the run_peanuts.py executable, in the
following way

⬇

Neutrinos: ”examples/example_slha.slha2”

The Solar node must contain a map with at least an entry for the neutrino
fraction fraction, as above. By default, this will compute and print the
probabilities of oscillation at Sun exit for all flavour eigenstates, but it
can be disabled with the entry probabilities: false. Additionally, one can
choose to print the total flux with flux:true (defaults to false) and the
distorted or undistorted spectrum with spectrum:”distorted” or
spectrum:”undistorted”, both of which are disabled by default. Lastly one can
select the specific solar model to use with the entry solar_model, the flux
file with flux_file and the spectrum files with the entry spectra. Note that
if the solar model is not known to PEANUTS, in addition to the path to the
relevant files, one must provide the location within those files where the
information can be found, that is the entries fluxrows, fluxcols, fluxscale,
distrow, radiuscol, densitycol and fractioncols, in the format described above
in Section 2.

The Earth node must contain a map with necessary entries to compute the
probability at some location on Earth. Consequently, it requires an entry for
the depth under the Earth’s surface, depth, and either a nadir angle eta or a
latitude. In addition, if the Solar node is not present in the YAML file, a
neutrino state is required as the state entry, as well as an entry to specify
in which basis it is, basis: ”flavour” (coherent) or basis: ”mass”
(incoherent). If the entry antinu: True is select, the computations will be
performed for an antineutrino. One can select a custom Earth density profile
by providing a density file with the entry density, as well as choose to use
either the numerical or analytical computation of the evolutor with evolution:
”numerical” or evolution: ”analytical” (the latter being the default). Lastly,
if the entry latitude is present, the probabilities will be computed
integrated over exposure, and thus one provides entries to modify the
normalization, days interval, number of samples, exposure file and exposure
angle (see Section 5 for the meaning of these quantities) with the options
exposure_normalized (True or False), exposure_time ([d1,d2]), exposure_samples
(number of samples), exposure_file (path) and exposure_angle (”Nadir”,
”Zenith” or ”CosZenith”), respectively. An example YAML file with only the
earth node can be seen below

⬇

Energy: 20

Neutrinos:

dm21: 7.42e-05

dm3l: 2.51e-03

theta12: 5.83638e-01

theta23: 8.5521e-01

theta13: 1.49575e-01

delta: 3.40339

Earth:

state: [0.4, 0.2, 0.6]

basis: ”flavour”

eta: 0.8

depth: 3000

By default, the results of a PEANUTS run will be printed to screen, but one
can redirect the ouput to a file by adding the Output node to the YAML file.
This node must contain a single string corresponding to the output file
location, which will be created if it does not exist, or ”stdout” to print to
screen (default). As an example, in order to redirect output to a file called
out.dat in the same directory as the executable, the following node must be
added to the YAML file

⬇

Output: ”out.dat”

Finally, PEANUTS allows the possibility to run simple grid scans of the input
parameters provided in the YAML file. This can be achieved easily by providing
a range of values instead of a single real value for the parameters. This
range can be provided as [min,max,step], where a specific step size is
selected, or [min,max], where the step will be computed so that a total of 10
samples for the parameter are produced. For instance, to scan over the energy
between 20 and 100 MeV, with a step of 10 MeV, one could add the following to
the YAML file.

⬇

Energy: [20,100,10]

PEANUTS can perform such simplistic grid scans on the neutrino energy (as seen
above), as well as on any neutrino parameter ($\Delta m_{21}^{2}$, $\Delta
m_{3l}^{2}$, $\theta_{12}$, $\theta_{13}$, $\theta_{23}$ and
$\delta$)101010Although it is technically possible to specify a range for
$\Delta m_{3l}^{2}$ that includes negative and positive values, we do not
recommend doing so, as the meaning of parameter changes from negative to
positive (as the ordering of mass eigenstates change), which is statistically
inconsistent. Hence we recommend scanning negative and positive values
separately., and on the nadir angle eta. As these are crude grid scans, we do
not recommend using this functionality for thorough scans of the full neutrino
parameter space. For that purpose one should use smart sampling algorithms
(see [57]), which can be easily interfaced with PEANUTS 111111This was used,
for instance, in an upcoming global study of neutrino oscillations by the
GAMBIT neutrino working group[58]..

###  6.1 Validation

We have validated the results from PEANUTS with those published by the SNO
experiment [35, 36], both against their predictions, e.g. fluxes and spectra,
and their measured data, e.g. exposure and observed events. Here we document
some of the validation tests that were performed. These can be reproduced by
running the run_SNO_test.py file shipped with PEANUTS, in the following way

⬇

run_SNO_test.py [options]

which allows for additional options to be supplied, such as -s/–solar <file>
to add a custom solar model, -d/–density <file> to add a custom earth density
profile, -f/–file <file> to modify the default neutrino parameters with an
SLHA file, and the ever present options -v/–verbose for debug output and
-h/–help for usage info.

In order to compare the probability results with those from SNO, we use the
same input values for the neutrino parameters, $\theta_{12}$, $\theta_{13}$,
$\Delta m_{21}^{2}$ and $\Delta m_{3l}^{2}=\Delta m_{31}^{2}$. The values of
$\theta_{23}$ and $\delta$ are, for the most part, irrelevant for this
comparison, but required by PEANUTS to build the full PMNS matrix, so for
those we take the values found in [18]. Thus, the values used here are

$\displaystyle\tan^{2}\theta_{12}=0.469,\quad\sin^{2}\theta_{13}=0.01,\quad\theta_{23}=0.85521,\quad\delta=3.4034$
$\displaystyle\Delta m_{21}^{2}=7.9\times 10^{-5},\quad\Delta
m_{31}^{2}=2.46\times 10^{-3}.$ (97)

Figure 6: Survival probability at the surface of the Sun for ${}^{8}B$ (left)
and $hep$ (right) neutrinos. The solid lines are the PEANUTS predictions, the
dashed lines are the digitised curves for the SNO experiment, from Fig. 6.3 in
[35].

With these values, we can compare in Fig. 6 the prediction from our code with
the SNO survival probabilities at the surface of the Sun for the ${}^{8}B$ and
$hep$ neutrinos, digitised from Fig. 6.3 in [35]. Naturally, for the same
Solar model, the predictions of PEANUTS for each neutrino fraction match very
well that reported by the SNO experiment.

Figure 7: Neutrino spectrum for the ${}^{8}B$ (left) and $hep$ (right)
neutrino fractions at the surface of the Sun. Solid lines ares are the PEANUTS
predictions for the distorted spectrum (with oscillations). The dotted lines
are the undistorted spectrum (no oscillations). Dashed lines are computing
distorting the spectrum with the survival probability from the digitised
curves for the SNO experiment, from Fig. 6.3 in [35].

In addition to comparing the survival probability, we can compare the
distorted (i.e. including oscillations) energy spectrum of neutrinos at the
surface of the Sun. Figure 7 shows the effect of the oscillation distortion on
the electron neutrino spectrum, as well as the comparison with the spectrum
distorted by the SNO survival probability. For this comparison we have used
the spectra in [52] for the ${}^{8}B$ fraction and [53] for the $hep$
fraction.

Figure 8: Oscillations of neutrino flavour eigenstates for $\eta=0$ (left) and
$\eta=\pi/2$ (right). The error quoted is the difference between the
analytical and numerical results.

In addition to validating the results against those of the SNO experiment, we
have also validated our analytical approximation from Section 4.1 against the
full numerical calculation described at the beginning of Section 4. Figure 8
shows the oscillation pattern of the neutrino flavour eigenstates for a
neutrino crossing all Earth shells, $\eta=0$ (left), and for a neutrino coming
from the horizon, $\eta=\pi/2$ (right), starting from a pure mass neutrino
eigenstate $\nu_{\alpha}=U_{\alpha 2}$. As expected, only in the first case
the oscillations are significant, since the neutrinos traverse the whole Earth
before reaching the detector, whereas in the second case the neutrino path
only crosses part of Earth’s crust between the surface and the detector, and
thus only minimal oscillation occurs. Since only the numerical computation can
provide the full trajectory, the oscillations shown in Fig. 8 correspond
exclusively to the numerical mode. However, for the same choice of parameters
we have also used the analytical computation for the final probabilities, and
the error reported in the figures precisely corresponds to the relative
difference between the numerical and analytical models121212Given
$(c_{e},c_{\mu},c_{\tau})$ the complex coefficients defining the final
(evolved) state (cf. eq. 19), the relative error is defined here as the norm
of the difference between the numerical and analytical values of it, divided
by the norm of their sum. This ensures that errors in both real and imaginary
parts of the solution are correctly taken into account.. The difference is
almost negligible, of the order of $10^{-4}$ for the night period (left), and
effectively zero during the day period. The left panel in Figure 9 confirms
this by showing the relative error as a function of the energy for the worse
case of the ones above (i.e. during the night) for $\eta=0$, where it can
clearly be noticed that the error is extremely small for all values of the
energy, and only approaches $\sim 10^{-3}$ at the worst for $E_{\nu}\sim
10^{2}$ MeV. This comparison serves as a validation that the approximated
analytical solution is a very good approximation and, since it is much faster,
can be used in place of the full numerical evaluation131313Notice that, at
energies above the TeV scale, neutrino inelastic scattering becomes relevant
for Earth-type densities [59]. PEANUTS does not currently include such
effects. Fig. 9 only compares the coherent forward scattering computed
numerically and analytically.. To quantify the speed increase gained with the
analytical implementation, we show in the right panel of Figure 9 the
computational time of the numerical and analytical computations as a function
of number of evaluations, for two values of the neutrino energy $E_{\nu}=10$
MeV and $E_{\nu}=100$ MeV. For small number of evaluations the computational
time is very similar between the analytical and numerical methods, but for a
number of evaluations $N\gtrsim 10$, the computational time for the numerical
method increases drastically and soon becomes computationally unfeasible. In
contrast the total CPU time of the analytical method remains constant for
increasing number of evaluations, and it is mostly dominated by the overhead
of the initialisation step. Only for $N\gtrsim 10^{5}$ the computation time
starts to increase noticeably, but still remains manageable up to large $N$.
To emphasise further the increase of speed with the analytical implementation,
the relative error shown on the left-hand panel of Figure 9 required over
$9.2\times 10^{3}$ seconds of CPU time to perform the numerical computations
for all energies on a 2,3 GHz Intel Core i7 quad-core, compared to 5 seconds
with the analytical method. Lastly, the right panel of Figure 9 also shows a
comparison with the state-of-the-art tool for neutrino oscillations, nuSQuIDS
[60]. For low number of evaluations, nuSQuIDS clearly outperforms PEANUTS, by
virtue of being written in C++, but for large number of evaluations, as early
as $N\gtrsim 5$ for low energies and $N\gtrsim 50$ for high energies, the
analytical implementation of PEANUTS becomes significantly faster and, again,
is the only feasible option for these large number of evaluations. It is worth
noting that this comparison is done in the worst case scenario, where $\eta=0$
and the neutrino path crosses almost the whole Earth. In less extreme
scenarios, the computational speed of the numerical implementation of PEANUTS
and the computations by nuSQuIDS are somewhat faster than shown.

Figure 9: Relative error between the numerical and analytical computations of
the probability after Earth regeneration as a function of the energy for a
chosen value of $\eta=0$ (left). Speed comparison of the analytical and
numerical computations by PEANUTS and the computations by nuSQuIDS for two
values of the neutrino energy (right).

Section 5 showed how to compute the ideal exposure time for a hypothetical
experiment. For specific experiments, however, the exposure is often provided
tabulated in bins of either the nadir angle $\eta$, the zenith angle $\theta$
or $\cos\theta$. In the case of the SNO experiment, the exposure is provided
in bins of $\cos\theta$, hence we convert it into bins of $\eta$ in order to
match fit our computations of the probability. We then show in Figure 10 the
exposure of the SNO experiment compared to the ideal case. As SNO is located
at a latitude of $46.475^{\circ}$, we can see that it matches very well the
ideal exposure at $46^{\circ}$, with a slight under-exposure during the day
and a slight over-exposure, during the night, which is consistent with the
livetime of the SNO experiment.141414This pattern is due to most of the
maintenance and calibration operations taking place during daytime [35].

Figure 10: Exposure of the SNO experiment in nadir angle $\eta$ (red),
compared to the ideal exposure for hypothetical experiments at various
latitudes. Coloured bands are as in Figure 5.

Finally, in order to match the computations of PEANUTS with that of the SNO
experiment, we reproduce the statistical fit of the oscillation parameters
$\theta_{12}$ and $\Delta m_{21}^{2}$ performed by the SNO experiment. Figure
11 shows the results of the fit as a profile likelihood ratio. The dark purple
star corresponds to the best fit point as found using PEANUTS, while the red
star is that reported by the SNO collaboration. Similarly the purple and blue
shared contours correspond to the 68% and 96% confidence intervals around the
best fit point for our results, whereas the red dashed contours are the same
from the SNO results. Note that this comparison, both our results and those
from the SNO experiment, corresponds only to Phase I of the experiment [36],
and for normal ordering of neutrino masses151515We have repeated the fit for
inverted ordering and found the results to be almost identical to that of
normal ordering, which is expected given that the oscillation of solar
neutrinos is largely independent of $\Delta m^{2}_{3l}$.. More details about
this comparison will appear in a global fit by the GAMBIT Neutrino working
group [58]. The results show a good match with those reported by the SNO
collaboration, with the best fit points laying very close to each other, at
almost identical locations on the $\tan^{2}\theta_{12}$ and only slightly
displaced to smaller values for $\Delta m_{21}^{2}$. The shape and reach of
the contours is also slightly different, which can be attributed to a minor
difference on the treatment of systematic uncertainties. It is crucial here to
emphasize that a parameter scan of this magnitude is only feasible with the
analytical implementation of PEANUTS, due to the large number of evaluations
required. For reference, the scan sampled around 380k parameter points, each
of which performed, on average, around 5k evaluations of the probability.

Figure 11: Profile likelihood ratio for the results of a statistical fit of
the oscillation parameters $\theta_{12}$ and $\Delta m_{21}^{2}$ using PEANUTS
(purple-blue) compared to the results reported by the Phase I of the SNO
experiment (red). The fit was performed using GAMBIT [61] and the figure was
generated with pippi [62].

##  7 Conclusions and Outlook

We have presented in this paper PEANUTS, a fast and flexible software to
automatically compute the energy spectra of solar active neutrinos, for
arbitrary solar models and custom Earth density profiles. PEANUTS assumes
adiabatic propagation of neutrinos within the Sun, and provides analytic
computation for the coherent evolution of active neutrinos while crossing the
Earth, thus completely avoiding any time-consuming numerical integration.
This, together with extensive use of pre-compiled and just-in-time compilation
optimisations, makes the software extremely fast and optimised for large-scale
parameter scans.

PEANUTS provides algorithms to automatically perform the full chain of
computations to simulate a solar neutrino experiment, as well as easy
individual access to the modules and functions for specific computations.
These include, for instance: mixing parameters in matter, incoherent flux at
Sun surface, evolved neutrino state after Earth crossing, Earth density
profile for given nadir angle, evolutor operator for given neutrino energy and
nadir angle, nadir exposure for an experiment between two arbitrary days of
the year, integrated probability of oscillation over a finite observation
time. In the present version of PEANUTS we focused on providing automatisation
for solar neutrinos, but the modularity of PEANUTS also allows a user to
employ, for example, the function Pearth_analytical to compute the evolution
of an atmospheric neutrino, or simulate the evolved spectra for hypothetical
neutrinos produced in solar flares.

PEANUTS can be run in _simple_ mode, for quick computations directly from the
command line, as well as in _expert_ mode, in which case the user provides a
set of comprehensive instructions in the form of a YAML file. The _expert_
mode is thus ideal for scripting, as it exploits all the possible
functionalities of PEANUTS, and to perform simple explorations of the
parameter space, as it natively allows the scanning of various input
parameters on a grid.

We extensively validated PEANUTS against the results of the SNO experiment,
and provide ready-to-run scripts to reproduce our results. We observed
excellent agreement for the survival probability at Sun surface, distorted
neutrino spectra for the ${}^{8}B$ and $hep$ channels, annual nadir exposure
and ability to reproduce the profile likelihood of the Phase I of the
experiment for the parameters $\theta_{12}$ and $\Delta m^{2}_{21}$. We also
validated our analytic computation of the evolved neutrino state with a
numerical solution, and find excellent agreement on a wide range of energies.
We performed a quick estimation of the CPU time scaling of the analytical and
numerical computations with the number of evaluations, and conclusively found
our analytical implementation to significantly outperform the numerical one
for as few as $\mathcal{O}(10)$ evaluations. We also compared our results with
those of nuSQuIDS and found good agreement as well. As with the numerical
implementation, we found that our analytical solution vastly outperforms
nuSQuIDS in terms of the computational time. PEANUTS thus presents a very fast
alternative to existing tools without appreciable loss of precision.
Therefore, due to its speed and user-friendly interface, we argue that PEANUTS
is a very capable and useful tool for the computation of the propagation of
neutrinos in the Sun and through Earth and a crucial addition to the software
toolbox of the neutrino physics community.

Concerning the limitations of the software, the current PEANUTS version
assumes adiabatic evolution within the Sun, which provides an excellent
approximation for solar neutrino energies given the currently allowed range of
oscillation parameters; the user is thus encouraged to check the validity of
this regime if working in non-standard scenarios. We also stress that the code
assumes coherent forward scattering for the propagation of neutrinos in
matter: for Earth-type densities, inelastic scattering becomes important at
energies above the TeV scale, and such effects must thus be taken into
account. We plan to incorporate a more general routine for the arbitrary
evolution of neutrinos in the Sun in a future version of PEANUTS. For
inelastic scattering, other softwares are currently available, cf. e.g. [63,
64, 65, 66, 60, 67]. Though PEANUTS v1.0 does not yet provide all the
functionalities that other software tools have, it vastly outperforms them in
speed, which can be a critical factor in e.g. global studies, where the number
of evaluations is very large.

Further improvements that will extend the scope of PEANUTS are under
investigation, most notably the implementation of algorithms for the automatic
computation of the atmospheric neutrino flux at a given location, as well as
the simulation of accelerator neutrinos beams.

## Acknowledgements

We would like to thank Aaron Vincent for helpful comments about this work,
Roberto Ruiz de Austri for testing and suggestions, and the rest of the GAMBIT
community for useful discussions. We also thank Will Handley for access to HPC
resources, as for this work we used the Cambridge Service for Data Driven
Discovery(CSD3), part of which is operated by the University of Cambridge
Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).
The DiRAC component of CSD3 was funded by BEIS capital funding via STFC
capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant
ST/R00689X/1. DiRAC is part of the National e-Infrastructure. ML acknowledges
financial support from the Alexander von Humboldt Foundation during the early
stages of the work, and thanks the UCLouvain Centre for Cosmology, Particle
Physics and Phenomenology (CP3) for hosting him during the final stage. TEG is
funded by the Deutsche Forschungsgemeinschaft (DFG) through the Emmy Noether
Grant No. KA 4662/1-1.

## References

  * [1] J.N. Bahcall, _NEUTRINO ASTROPHYSICS_ (1989).
  * [2] J.N. Bahcall, A.M. Serenelli and S. Basu, _New solar opacities, abundances, helioseismology, and neutrino fluxes_ , _Astrophys. J. Lett._ 621 (2005) L85 [astro-ph/0412440].
  * [3] J.N. Bahcall and R. Davis, Jr., _AN ACCOUNT OF THE DEVELOPMENT OF THE SOLAR NEUTRINO PROBLEM_ , .
  * [4] R. Davis, Jr., D.S. Harmer and K.C. Hoffman, _Search for neutrinos from the sun_ , _Phys. Rev. Lett._ 20 (1968) 1205.
  * [5] J.N. Bahcall, W.A. Fowler, I. Iben, Jr. and R.L. Sears, _Solar neutrino flux_ , _Astrophys. J._ 137 (1963) 344.
  * [6] R.L. Sears, _Helium Content and Neutrino Fluxes in Solar Models_ , _Astrophys. J._ 140 (1964) 477.
  * [7] P. Pochoda and H. Reeves, _A revised solar model with a solar neutrino spectrum_ , _Planet. Space Sci._ 12 (1964) 119.
  * [8] J.N. Bahcall, _Solar neutrinos. I: Theoretical_ , _Phys. Rev. Lett._ 12 (1964) 300.
  * [9] D. Ezer and A.G.W. Cameron, _A study of solar evolution_ , _Can. J. Phys._ 43 (1965) 1497.
  * [10] D. Ezer and A.G.W. Cameron, _Solar evolution with varying g_ , _Can. J. Phys._ 44 (1966) 593.
  * [11] J.N. Bahcall, M. Cooper and P. Demarque, _Dependence of the 8̂{B} Solar Neutrino Flux on Heavy Element Composition_, _ApJ_ 150 (1967) 723.
  * [12] G. Shaviv, J.N. Bahcall and W.A. Fowler, _Dependence of the 8̂{B} Solar Neutrino Flux on the Rate of the Reaction 3̂He(3̂He,2p)4He_, _ApJ_ 150 (1967) 725.
  * [13] J.N. Bahcall, N.A. Bahcall, W.A. Fowler and G. Shaviv, _Solar neutrinos and low-energy nuclear cross sections_ , _Phys. Lett. B_ 26 (1968) 359.
  * [14] SNO collaboration, _Measurement of the rate of $\nu_{e}+d\to p+p+e^{-}$ interactions produced by 8B solar neutrinos at the Sudbury Neutrino Observatory_, _Phys. Rev. Lett._ 87 (2001) 071301 [nucl-ex/0106015].
  * [15] Z. Maki, M. Nakagawa and S. Sakata, _Remarks on the unified model of elementary particles_ , _Prog. Theor. Phys._ 28 (1962) 870.
  * [16] B. Pontecorvo, _Neutrino Experiments and the Problem of Conservation of Leptonic Charge_ , _Zh. Eksp. Teor. Fiz._ 53 (1967) 1717.
  * [17] V.N. Gribov and B. Pontecorvo, _Neutrino astronomy and lepton charge_ , _Phys. Lett. B_ 28 (1969) 493.
  * [18] I. Esteban, M.C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, _The fate of hints: updated global analysis of three-flavor neutrino oscillations_ , _JHEP_ 09 (2020) 178 [2007.14792].
  * [19] P.F. de Salas, D.V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C.A. Ternes et al., _2020 global reassessment of the neutrino oscillation picture_ , _JHEP_ 02 (2021) 071 [2006.11237].
  * [20] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri and A. Palazzo, _Global constraints on absolute neutrino masses and their ordering_ , _Phys. Rev. D_ 95 (2017) 096014 [2003.08511].
  * [21] BOREXINO collaboration, _Experimental evidence of neutrinos produced in the CNO fusion cycle in the Sun_ , _Nature_ 587 (2020) 577 [2006.15115].
  * [22] Borexino collaboration, _First Simultaneous Precision Spectroscopy of $pp$, 7Be, and $pep$ Solar Neutrinos with Borexino Phase-II_, _Phys. Rev. D_ 100 (2019) 082004 [1707.09279].
  * [23] BOREXINO collaboration, _Comprehensive measurement of $pp$-chain solar neutrinos_, _Nature_ 562 (2018) 505.
  * [24] Hyper-Kamiokande Proto collaboration, _Solar neutrino physics at Hyper-Kamiokande_ , _PoS_ ICRC2019 (2020) 1037.
  * [25] IceCube collaboration, _Search for GeV neutrino emission during intense gamma-ray solar flares with the IceCube Neutrino Observatory_ , _Phys. Rev. D_ 103 (2021) 102001 [2101.00610].
  * [26] Super-Kamiokande collaboration, _Searching for neutrinos from solar flares across solar cycles 23 and 24 with the Super-Kamiokande detector_ , 2210.12948.
  * [27] IceCube collaboration, _Recent Progress in Solar Atmospheric Neutrino Searches with IceCube_ , _PoS_ ICRC2021 (2021) 1174 [2107.13696].
  * [28] IceCube collaboration, _Latest results and sensitivities for solar dark matter searches with IceCube_ , _PoS_ ICRC2017 (2018) 912.
  * [29] Super-Kamiokande collaboration, _Determination of solar neutrino oscillation parameters using 1496 days of Super-Kamiokande I data_ , _Phys. Lett. B_ 539 (2002) 179 [hep-ex/0205075].
  * [30] J.N. Bahcall and C. Pena-Garay, _Global analyses as a road map to solar neutrino fluxes and oscillation parameters_ , _JHEP_ 11 (2003) 004 [hep-ph/0305159].
  * [31] Super-Kamiokande collaboration, _Testing Non-Standard Interactions Between Solar Neutrinos and Quarks with Super-Kamiokande_ , 2203.11772.
  * [32] Super-Kamiokande collaboration, _Search for solar electron anti-neutrinos due to spin-flavor precession in the Sun with Super-Kamiokande-IV_ , _Astropart. Phys._ 139 (2022) 102702 [2012.03807].
  * [33] L. Wolfenstein, _Neutrino Oscillations in Matter_ , _Phys. Rev. D_ 17 (1978) 2369.
  * [34] S.P. Mikheyev and A.Y. Smirnov, _Resonance Amplification of Oscillations in Matter and Spectroscopy of Solar Neutrinos_ , _Sov. J. Nucl. Phys._ 42 (1985) 913.
  * [35] N.F. Fiúza de Barros, _Precision Measurement of Neutrino Oscillation Parameters: Combined Three-phase Results of the Sudbury Neutrino Observatory_ , Ph.D. thesis, Lisbon U., 2011.
  * [36] SNO collaboration, _Determination of the $\nu_{e}$ and total 8B solar neutrino fluxes with the Sudbury neutrino observatory phase I data set_, _Phys. Rev. C_ 75 (2007) 045502 [nucl-ex/0610020].
  * [37] J.N. Bahcall, M.H. Pinsonneault and S. Basu, _Solar models: Current epoch and time dependences, neutrinos, and helioseismological properties_ , _Astrophys. J._ 555 (2001) 990 [astro-ph/0010346].
  * [38] S.P. Mikheev and A.Y. Smirnov, _NEUTRINO OSCILLATIONS IN MATTER_ , in _International Symposium on Weak and Electromagnetic Interactions in Nuclei_ , pp. 405–415, 1987.
  * [39] E. Lisi and D. Montanino, _Earth regeneration effect in solar neutrino oscillations: An Analytic approach_ , _Phys. Rev. D_ 56 (1997) 1792 [hep-ph/9702343].
  * [40] D. Bruss and L.M. Sehgal, _DISTINGUISHING A COHERENT FROM AN INCOHERENT MIXTURE OF NEUTRINO FLAVORS_ , _Phys. Lett. B_ 216 (1989) 426.
  * [41] S.P. Mikheev and A.Y. Smirnov, _Neutrino Oscillations in a Variable Density Medium and Neutrino Bursts Due to the Gravitational Collapse of Stars_ , _Sov. Phys. JETP_ 64 (1986) 4 [0706.0454].
  * [42] P.B. Denton, H. Minakata and S.J. Parke, _Compact Perturbative Expressions For Neutrino Oscillations in Matter_ , _JHEP_ 06 (2016) 051 [1604.08167].
  * [43] P.B. Denton and S.J. Parke, _Addendum to “Compact perturbative expressions for neutrino oscillations in matter”_ , 1801.06514.
  * [44] A. Ioannisian and S. Pokorski, _Three Neutrino Oscillations in Matter_ , _Phys. Lett. B_ 782 (2018) 641 [1801.10488].
  * [45] P.B. Denton, H. Minakata and S.J. Parke, _Comment on 1801.10488v3_ , .
  * [46] G. Fantini, A. Gallo Rosso, F. Vissani and V. Zema, _Introduction to the Formalism of Neutrino Oscillations_ , _Adv. Ser. Direct. High Energy Phys._ 28 (2018) 37 [1802.05781].
  * [47] T.-K. Kuo and J.T. Pantaleone, _Neutrino Oscillations in Matter_ , _Rev. Mod. Phys._ 61 (1989) 937.
  * [48] H.A. Bethe and J.N. Bahcall, _Solar neutrinos and the Mikheev-Smirnov-Wolfenstein theory_ , _Phys. Rev. D_ 44 (1991) 2962.
  * [49] D. Harley, T.-K. Kuo and J.T. Pantaleone, _Solar neutrinos with three flavor mixings_ , _Phys. Rev. D_ 47 (1993) 4059 [hep-ph/9208241].
  * [50] E.K. Akhmedov, M.A. Tortola and J.W.F. Valle, _A Simple analytic three flavor description of the day night effect in the solar neutrino flux_ , _JHEP_ 05 (2004) 057 [hep-ph/0404083].
  * [51] N. Vinyoles, A.M. Serenelli, F.L. Villante, S. Basu, J. Bergström, M.C. Gonzalez-Garcia et al., _A new Generation of Standard Solar Models_ , _Astrophys. J._ 835 (2017) 202 [1611.09867].
  * [52] C.E. Ortiz, A. Garcia, R.A. Waltz, M. Bhattacharya and A.K. Komives, _Shape of the B-8 alpha and neutrino spectra_ , _Phys. Rev. Lett._ 85 (2000) 2909 [nucl-ex/0003006].
  * [53] J.N. Bahcall, _Gallium solar neutrino experiments: Absorption cross-sections, neutrino spectra, and predicted event rates_ , _Phys. Rev. C_ 56 (1997) 3391 [hep-ph/9710491].
  * [54] W.T. Winter, S.J. Freedman, K.E. Rehm and J.P. Schiffer, _The B-8 neutrino spectrum_ , _Phys. Rev. C_ 73 (2006) 025503 [nucl-ex/0406019].
  * [55] F.J. Dyson, _The Radiation theories of Tomonaga, Schwinger, and Feynman_ , _Phys. Rev._ 75 (1949) 486.
  * [56] T. Ohlsson and H. Snellman, _Three flavor neutrino oscillations in matter_ , _J. Math. Phys._ 41 (2000) 2768 [hep-ph/9910546].
  * [57] S.S. AbdusSalam et al., _Simple and statistically sound recommendations for analysing physical theories_ , _Rept. Prog. Phys._ 85 (2022) 052201 [2012.09874].
  * [58] GAMBIT Neutrino Workgroup collaboration, _Global analysis of neutrino oscillations with gambit_ , _in preparation_ (2023) .
  * [59] IceCube collaboration, _Measurement of the multi-TeV neutrino cross section with IceCube using Earth absorption_ , _Nature_ 551 (2017) 596 [1711.08119].
  * [60] C.A. Argüelles, J. Salvado and C.N. Weaver, _nuSQuIDS: A toolbox for neutrino propagation_ , _Comput. Phys. Commun._ 277 (2022) 108346 [2112.13804].
  * [61] GAMBIT collaboration, _GAMBIT: The Global and Modular Beyond-the-Standard-Model Inference Tool_ , _Eur. Phys. J. C_ 77 (2017) 784 [1705.07908].
  * [62] P. Scott, _Pippi - painless parsing, post-processing and plotting of posterior and likelihood samples_ , _Eur. Phys. J. Plus_ 127 (2012) 138 [1206.2245].
  * [63] A. Gazizov and M.P. Kowalski, _ANIS: High energy neutrino generator for neutrino telescopes_ , _Comput. Phys. Commun._ 172 (2005) 203 [astro-ph/0406439].
  * [64] A.C. Vincent, C.A. Argüelles and A. Kheirandish, _High-energy neutrino attenuation in the Earth and its associated uncertainties_ , _JCAP_ 11 (2017) 012 [1706.09895].
  * [65] A. Garcia, R. Gauld, A. Heijboer and J. Rojo, _Complete predictions for high-energy neutrino propagation in matter_ , _JCAP_ 09 (2020) 025 [2004.04756].
  * [66] I. Safa, J. Lazar, A. Pizzuto, O. Vasquez, C.A. Argüelles and J. Vandenbroucke, _TauRunner: A public Python program to propagate neutral and charged leptons_ , _Comput. Phys. Commun._ 278 (2022) 108422 [2110.14662].
  * [67] D. Garg et al., _Neutrino propagation in the Earth and emerging charged leptons with nuPyProp_ , _JCAP_ 01 (2023) 041 [2209.15581].