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

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.


Introduction
Solar neutrinos represent the most abundant source of neutrinos on Earth, with a flux of the order of 6 × 10 10 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 8 B 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-changing 1 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 opensource, fast and flexible package to compute the neutrino oscillation dynamics in all the abovementioned 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 com- 1 Slow or fast here refers to the importance of variation of matter density over a neutrino oscillation length scale.
putations 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.PEANUTS2 is provided as an open-source Python package, under the GPL-3.0license.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 vacuum 3 .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 (θ 12 , θ 13 , θ 23 , δ) 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 θ 12 , θ 13 with their matter-rotated ones [42,43] where θ 12 , θ 13 are the PMNS mixing angles in vacuum, θ 12 , θ 13 the corresponding ones in matter, 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, is the θ 13 -modified matter potential and ∆m The matrix T is then simply defined as 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]) where the wavenumber k is useful in the relation The probability of producing a neutrino mass eigenstate i in matter is thus |T αi (E, n e )| 2 , where α 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 θ 12 , θ 13 , ∆m 2  21 (eV 2 ), ∆m 2 3ℓ (eV 2 ), E (MeV) and n e (mol/cm 3 ).Note that in PEANUTS the input variable ∆m 3ℓ has different meanings according to the ordering of neutrino mass eigenstates.In normal ordering (NO) l = 1, and thus ∆m .PEANUTS also provides accessible functions for the mixing angles in matter θ 12 and θ 13 , as well as useful quantities such as ∆m 2  ee 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 ≡ R/R ⊙ , where R ⊙ 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 with the normalization The computation of f (r) assumes a specific solar model.In our validation of PEANUTS we use the BS05(AGS,OP) model [2] 5 , 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 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 fluxes.By default, the following are provided and taken from: pp and hep [52], 8 B [53], 13 N , 15 O and 17 F [54], 7 Be [55].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 [56], 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)

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 with α = e, µ, τ 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 states6 implying that the probability of observing a neutrino of flavour α from a mass eigenstate i is given by Thus, the probability for a solar neutrino to manifest as flavour α is given by 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 θ ij and CP phase δ CP as 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 ∆, which allow the PMNS mixing matrix U to be expressed as and are defined as 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

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, a generic neutrino state at time7 t can be expressed in terms of the state at time t 0 by evolving it with an appropriate evolutor operator Û(t, t 0 ) The generic state |ν, t⟩ can be expressed as a linear superposition of pure flavour eigenstates, where c α (t) are complex numbers, implying that the probability of it to interact as a neutrino with flavour α at time t will be given by |c α (t)| 2 .From Eq. 17 it follows that the evolved probability amplitudes are given by where U αβ (t, t 0 ) are the matrix elements of the evolutor operator in flavour basis.The determination of the evolutor operator 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 mass eigenstate expressed as linear combination of flavour eigenstates is which implies a transition amplitude from (evolved) mass to flavour eigenstate Putting everything together, the final probability for a solar neutrino to manifest as α flavour while crossing the Earth is given by where t 0 is defined at the time of neutrino crossing the Earth surface.The interpretation of Eq. ( 22) is the following: U βi are the coefficients of the mass eigenstate i expressed as linear combination of flavour eigenstates, |ν i ⟩ = U βi |ν β ⟩, and U αβ (t, t 0 )U βi = ⟨ν α |ν i , t⟩ is the transition amplitude from the evolved mass eigenstate i to interaction eigenstate α.Finally, each probability | ⟨ν α |ν i , t⟩ | 2 is multiplied by the weight of the mass eigenstate i in the incoherent solar flux, P ⊙ νe→ν i (E).The probability of oscillation for each flavour eigenstate in eq. ( 22) is implemented in PEANUTS by the function Pearth, with signature Pearth(nustate, density, pmns, DeltamSq21, DeltamSq3l, E, eta, depth, mode="analytical", massbasis=True, 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 massbasis defines the basis of the neutrino state.If massbasis=False the initial state is assumed to be a coherent one, expressed in flavour basis, with nustate defining the complex coefficients c α in eq. ( 18), and final probabilities computed by squaring the coefficients evolved as in eq. ( 19).If massbasis=True, the initial state is assumed to be an incoherent superposition of mass eigenstates, with nustate defining the real weights P ⊙ νe→ν i (E) and final probabilities computed as in eq. ( 22).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 α (t) from a coherent (flavour basis) initial neutrino state, from eq. 19, 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.

Evolutor operator
The evolutor operator Û is defined by the equations The equation of motion for Û(t, t 0 ) can be derived from the Schrödinger equation implying The Schrödinger equation in flavour basis takes the form The explicit expression for H αβ (t) can be readily derived in vacuum (where there is no time dependence) where E i is the energy of the mass eigenstate |ν i ⟩.In the presence of matter the Hamiltonian matrix elements receive an additional (time dependent) term, cf.Section 4; Eq. ( 27) determines its vacuum CP-structure, following the definition of the neutrino fields/states in Eq.s (12,13).
In PEANUTS the probability of oscillation through vacuum, and its evolved state, are computed with the functions Pvacuum(nustate, pmns, DeltamSq21, DeltamSq3l, E, L, antinu=False, massbasis=True) vacuum_evolved_state(nustate, pmns, DeltamSq21, DeltamSq3l, E, L, antinu=False) with similar arguments as Pearth above, with the notable difference of the oscillation length or baseline L, to be provided in km.
Once the Hamiltonian matrix elements are known, they can be used to derive the evolutor ones.
The formal solution is where T is the time-order operator.Eq. ( 28) 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 [57] which allows for an approximate solution obtained by truncating Eq. ( 29) 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.

Earth matter regeneration
Due to the incoherent nature of the solar neutrino flux, its relative weights in mass eigenstates components, P ⊙ νe→ν i (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 ν 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] 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.1: Values of the parameters for the electron density expressed as N j (r) = α j + β j r 2 + γ j r 4 with [N ] = mol/cm 3 , for each of the Earth internal shells, as derived in [39].The radial distance r is normalised to the radius of Earth.
The parametrisation in Eq. ( 30) is valid for radial trajectories, i.e. paths crossing the center of the Earth.For a path forming a nadir angle η with the radial trajectory (cf.Fig. 2), the parametrisation is functionally invariant, with modified coefficients where the trajectory coordinate x is defined as the distance from the trajectory mid-point, i.e. x = r 2 − sin 2 η.The Earth density profiles for some example values of the nadir angle η are reported in Fig. 3.
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 where R E is the (not rescaled) Earth radius, R E = 6.371 the scenario of an underground detector with respect to the treatment in [39]: the first is that the measured nadir angle η differs from the angle η ′ = arcsin(r det sin η) that one would measure at the Earth surface, for same neutrino trajectory, cf. Figure (4).This implies that it is the η ′ angle that should be used in Eq. ( 31) to compute the value of the electron density profile along the neutrino trajectory, and not the value η measured by the experiment.The second modification is that matter effects are present even for values of η ≥ π/2, if the detector is underground.The contribution to the trajectory from the outer layer (between Earth surface and detector shell) is given by ∆x reduces to h for η = 0, π, and attains the maximum value h(2 − h) for η = π/2 (this is approximately 160 km for SNO).
In general, not every shell is crossed by solar neutrinos: a shell i is crossed by a neutrino trajectory with nadir angle η if r i > r det sin η.The value of the trajectory coordinate at each shell crossing is  given by PEANUTS implements the Earth density as a Python class called EarthDensity, with signature EarthDensity(density_file=None, tabulated_density=False, custom_density=False) where the optional arguments density_file, tabulated_density and custom_density allow the user to supply any Earth density profile.If supplied by file, using the density_file argument, the density profile is either expected as a table with columns {r j , α j , β j , γ j }, following eq.( 30) (additional optional columns δ 2n j , with n ≥ 3, are allowed, corresponding to higher orders in the polynomial expansion of the density); or if tabulated_density=True it is expected as a two-column table of radii and density, where each entry will be treated as a new layer of constant density.Consequently PEANUTS can work with an arbitrary number of shells and arbitrary densities.In addition to supplying the density profile by file, it is also possible to provide an analytical expression for the density.Upon selecting the optional argument custom_density=True in the EarthDensity constructor, PEANUTS will use the earth density computed by the member function custom_density(r), where a user can implement their own analytical density profile.However, custom density profiles and higher orders in the polynomial expansion are only fully used when computing the oscillation probabilities numerically; the analytical computation, described below in Section 4.1, relies on the density described as in eq. ( 30), and thus any density provided will be Taylor expanded and truncated to fit that form, when used for the analytical computation of the oscillation probability.Regardless of the source and form of the density profile, the EarthDensity class provides methods to compute the value of the density at given coordinate x, as well as, when appropriate, the modified coefficients α ′ j , β ′ j and γ ′ j (and δ ′ 2n j if needed), corresponding to each Earth shell, from the radii and nadir angle η using Eq. ( 31) (or modified to accommodate higher orders).

Neutrino propagation Hamiltonian
The propagation Hamiltonian for an ultrarelativistic neutrino propagating in a medium with electron density n e (x) is, in the flavour basis [46] with For antineutrinos, the same Eq.( 37) holds, with the replacements To streamline computations and statistical analysis, it is more convenient to redefine the Hamiltonian by subtracting a constant term with j = 1 for normal ordering (NO) and j = 2 for inverted ordering (IO), such that so that we can use as free parameters ∆m 2 21 > 0 and ∆m 2 3ℓ , with ℓ = 1 for NO (∆m 2 31 > 0) and ℓ = 2 for IO (∆m 2 32 < 0).In the following we keep using the notation 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 (43,44).
With the parametrisation in Eq. ( 15) and using [∆, Hamiltonian can be rewritten as Notice that H does not depend on θ 23 , δ.
Given that R 23 and ∆ do not depend on position, they can be factorised in the time-ordered definition of the evolutor operator + . . .( other terms of the Dyson series for n → ∞) A numerical solution for the evolutor can be obtained by solving eq. ( 50) or by resolving the differential equation in eq. ( 25).PEANUTS offers a numerical evaluation of the evolutor, via the function Pearth_numerical(nustate, density, pmns, DeltamSq21, DeltamSq3l, E, eta, depth, massbasis=True, 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.

Perturbative expansion of the neutrino propagation Hamiltonian
We are interested in an expression for the operator in Eq. ( 50) 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 H can be divided in a kinetic and a matter dependent terms where Hk does not depend on x.
To work out a perturbative expression for U8 it is convenient to express the electron density as a perturbation along its mean value along the path [39] n e (x) = ne + δn(x), ne = 1 from which it follows We can analogously divide the Hamiltonian into a zeroth order term and a perturbation where again H0 does not depend on x.
The evolutor for a constant Hamiltonian H can generally be expressed in a closed form [58] e −i Hx = ϕ where T = H − Tr( H)1/3 is a traceless matrix and λ a are the roots of the characteristic equation with Finally By noticing that the full dependence on x in e −i Hx is now contained within the scalar functions e −iλax , the first order correction in Eq. ( 56) can be computed as where λa = λ a + Tr( H)/3 and we defined For a path fully contained within one shell we can parametrise where α′ = α ′ − ne , implying that I ab (x 2 , x 1 ) can be expressed analytically in closed form.
Summarising, we can perturbatively expand the evolutor operator as In PEANUTS we implement these perturbative expressions to compute U in Eq. ( 51) at first order in perturbation theory, using a "reduced" mixing matrix Ũ = R 13 R 12 , with he function 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 = α, b = β and c = γ.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 δ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 where x i is a generic point x 1 < x i < x 2 contained on the original path.It can be shown that [39] The consequences of Eq.s (70, 71) are twofold: first, for a path starting at x = x i , with 0 ≤ x i < x 1 , crossing n shells with boundaries at trajectory coordinate (0, x 1 , x 2 , . . ., x n ), and ending at the point x = x f , with x n−1 < x f ≤ x n , the full evolutor can be expressed as 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) Notice that the final evolutor is only a function of η, 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 ≤ η < π/2 where x det = r det cos η.For π/2 ≤ η ≤ π 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, We fix for simplicity the electron density value to the one at Earth surface, n 1 = 1.67 mol/cm 3 .Having assumed constant density along the path, the evolutor is simply given by To obtain the full dynamics and thus the full evolutor U, PEANUTS computes the time-ordered product of "reduced" evolutors Ũ from above, and inputs it in Eq. ( 50) to re-introduce the dependence on θ 23 and δ.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, massbasis=True, antinu=False) and the evolved coefficients with evolved_state_analytical(nustate, density, pmns, DeltamSq21, DeltamSq3l, E, eta, depth, antinu= False)

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 where τ d is the daily time, τ h the hourly time and η the Sun nadir angle at detector location.The integration in Eq. ( 77) 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 η [39] where W (η) is a normalized weight function representing the fraction of time in which the experiment collected data at nadir angle η.For real experiments, W (η) 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 (η) analytically, for the ideal case of an experiment continuously taking data between days τ d 1 and τ d 2 [39].This is done by changing integration variables By normalising the daily and hourly times to the interval [0, 2π] with τ d = 0 at winter solstice and τ h = 0 at the middle of the night, it is possible to express with λ the detector latitude and δ S the Sun declination, given by with i = 0.4091 rad being the Earth inclination.With these definitions it is possible to perform the integral defining W (η) in Eq. ( 80); it is convenient to restrict τ d within the interval [0, π] (the alternative case can be easily derived from this one by using the symmetry of the orbit) and to change the integration variable from τ d to T = cos(τ 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 ∈ [−1, 1] is the interval where T = cos(τ d ) is defined, ii) T ∈ [sin(λ − η)/ sin(i), sin(λ + η)/ sin(i)] is the range where T can take values for fixed values of η, λ, i, iii) T ∈ [cos(τ d 2 ), cos(τ d 1 )] is the observation time.If the intersection of the three intervals is null then W (η) will vanish for that given combination of λ, η values.
The exposure function W (η) 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 (η) for ns samples (default 1000) in nadir angle η, 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. ( 80), 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 (η) for one full year of exposure for three ideal detectors located at latitudes λ = 0, 45 Annual nadir exposure for an experiment at various latitudes = 0° = 45° = 89°F igure 5: Weight of the nadir angle exposure for an ideal experiment located at latitude λ, 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 η, but in some cases, the experiments provide the exposure tabulated in either the zenith angle θ or cos θ.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 η.
Lastly, PEANUTS provides the averaged probability of oscillation at detector location taking into account the finite time exposure, as in eq. ( 78), 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 η < π/2, with daynight ="night", or the day period η ≥ π/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.

Quick start guide
PEANUTS is an open-source software written in Python, and as such we expect it to run seamlessly in every Python 3 environment, irrespective of the architecture.Nevertheless it has been tested thoroughly on Linux and Mac OS X systems.An example installation guide using the conda environment can be found in Appendix A. 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 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.pyand run_prob_earth.py.
The first script computes the probability of neutrinos at the surface of the Sun, and has a signature where <energy> and <fraction> are mandatory arguments that refer to the neutrino energy, in MeV, and fraction respectively.The list of available neutrino fractions is ["pp", "hep", "7Be", "8B", "13N" , "15O", "17F"] 10 .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.pycomputes the probability of neutrinos at a given experimental location below the Earth's surface.It is used in the following way 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.The neutrino state must be given as comma-separated list of three real or complex numbers, without spaces in between them, e.g.run_prob_earth.py-f 0.1+0.03j,0.6+0.05j,0.09-0.9j... Other necessary arguments are the neutrino energy <energy>, in MeV, the nadir angle, <eta>, in radians, and the depth of the experiment <depth>, in meters.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.pyexecutable 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 The YAML file must contain the Energy and Neutrinos nodes, as well as either the Vacuum, Solar or Earth nodes.The Solar and Earth nodes can appear simultaneously to combine the effects, whereas the Vacuum node cannot be combined with either.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.pyexecutable, 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.
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 (see below for an example), 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 provide a user-defined Earth density profile by providing a density file with the entry density, and indicate whether the density comes from tabulated data, with the entry tabulated_density or from a custom analytical expression, with the entry custom_density.It is also possible to choose 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 The Vacuum node must contain the minimal requirements to compute oscillations in vacuum.This means it requires an input neutrino state and its basis, as for the earth oscillations above.It must also contain the distance traveled as baseline, in km.Lastly, optionally one can request vacuum oscillations for antineutrinos with antinu: True, and to switch on or off the printing of the probabilities or the evolved state with the probabilities and evolved_state options, respectively.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]  The solid lines are the PEANUTS predictions, the dashed lines are the digitised curves for the SNO experiment, from Fig. 6.3 in [35].
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, η = 0 (left), and for a neutrino coming from the horizon, η = π/2 (right), starting from a pure mass neutrino eigenstate ν α = U α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 models 13 .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 η = 0, where it can clearly be noticed that the error is small for all values of the energy, and only approaches ∼ 10 −2 at the worst for E ν ∼ 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 evaluation 14 .For completeness, for the production of Fig. 9, the unitarity of the evolution matrix was confirmed, and no unphysical solutions were found.Figure 9 also shows a comparison with the results from the state-of-the-art tool for neutrino oscillations nuSQuIDS [62], where it can be seen that the agreement is also good with a relative error consistently below 10 −2 .Therefore, the massive increase in computational speed provided by PEANUTS certainly justifies the small loss in precision. 15o quantify this 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 ν = 10 MeV and E ν = 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 ≳ 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 ≳ 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 × 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 nuSQuIDS.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 ≳ 5 for low energies and N ≳ 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 η = 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.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 η, the zenith angle θ or cos θ.In the case of the SNO experiment, the exposure is provided in bins of cos θ, hence we convert it into bins of η 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 • , we can see that it matches very well the ideal exposure at 46 • , 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. 16inally, in order to match the computations of PEANUTS with that of the SNO experiment, we reproduce the statistical fit of the oscillation parameters θ  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 95% 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 masses 17 .Details about this comparison will appear in a global fit by the GAMBIT Neutrino working group [60].The results show a decent match with those reported by the SNO collaboration, with the best fit points laying very close to each other.The shape and reach of the contours is larger in our study, which can be attributed to a slight 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.

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 largescale 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 readyto-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 θ 12 and ∆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 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.[69,70,71,72,62,73].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.This command will take care of installing both the pip and conda packages with versions as defined in peanuts env.yml.

Figure 2 :
Figure 2: Earth section showing the different shells and the path of a solar neutrino having nadir angle η.The trajectory coordinate x is also schematised.

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

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

Figure 6 :
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].

Figure 7 :
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].

Figure 9 :
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 η = 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).
12 and ∆m 2 21 performed by the SNO

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

Figure 11 :
Figure11: Profile likelihood ratio for the results of a statistical fit of the oscillation parameters θ 12 and ∆m2  21 using PEANUTS (purple-blue) compared to the results reported by the Phase I of the SNO experiment (red).The fit was performed using GAMBIT[67] and the figure was generated with pippi[68].