The Sun
The averaged solar survival probability. Top: the electron number density along a radial ray of the BS2005-AGS,OP solar model [BSB05]. The marked densities are those of the MSW resonance, \(\sqrt{2} G_F n_e = \Delta m^2_{21}\cos 2\theta_{12}/2E\), at three energies: a 20-MeV neutrino meets it a third of the way out; a 1-MeV neutrino never does. Bottom: the phase-averaged survival probability, the averaged limit on a varying profile, computed with average=True for four Hamiltonians on the same profile. The two sterile cases lie lowest, \(3+2\) below \(3+1\), because each extra state takes a share of the flux that does not return. Listing Averaged solar probabilities generates the data. See notebooks #10, #13, and #28.
Figure The averaged solar survival probability shows the electron density of the BS2005-AGS,OP solar model [BSB05], and the averaged survival probability it produces for four Hamiltonians. Along a radial ray, the density falls smoothly by five orders of magnitude, with no discontinuity. On the way out, a 5-MeV neutrino accumulates a phase of \(3 \cdot 10^4\) radians between the first two eigenstates, and of \(9 \cdot 10^5\) radians between the first and the third. So the survival probability at the surface runs through a full oscillation each time the production point moves by about 160 km, or the energy by a few parts in \(10^4\). Yet, neutrinos are produced over tens of thousands of kilometers, and no detector resolves energy that finely. Therefore, what a solar experiment measures is the phase-averaged probability of Phase-averaged probabilities.
The energy dependence of the probability in Figure The averaged solar survival probability is set by where the resonance density sits relative to the density at production. The resonance density falls with energy as \(1/E\). At 20 MeV, the profile falls to the resonance density a third of the way out. A neutrino born at the center, at ten times that density, is mostly the matter eigenstate that turns into \(\nu_2\) on the way out, so it leaves the Sun as \(\nu_2\), and its average is \(\lvert \mathbb{R}_{e2} \rvert^2 \approx 0.30\). At 1 MeV, the resonance density lies above the central density, so the neutrino is born below it and never crosses it. Matter then only lowers the average a little, to 0.51, from the vacuum value, 0.55, which the curve reaches at 0.1 MeV.
The probability in Figure The averaged solar survival probability shows the passage between the two regimes, over an energy range that spans the spectrum of neutrinos from \(^8\)B decay. Along the ray, the density changes slowly compared with the local oscillation length, so a neutrino stays in the eigenstate it was born in: \(P^{\rm cross}\) in the averaged limit on a varying profile is the identity, at every energy. What remains are the two eigenbases at the ends of the ray. The one at the surface is the vacuum mixing matrix, the same for any profile. The one at the other end is set by the density where the neutrino is born.
Listing Averaged solar probabilities computes the four curves of Figure The averaged solar survival probability. The BS2005-AGS,OP table is one of the twelve standard solar models that ship with Magνs (Comparing solar models); every Sun wrapper takes it by name through density_profile. Without that keyword, a wrapper uses its default, 'exp', an exponential fit to the solar density (Comparing solar models). For the sterile states, it also reads the neutron-to-proton ratio from the table (Comparing solar models).
The only difference from an ordinary probability call is average=True (Asking for the averaged probability). With it, Magνs does not propagate the phase along the ray. It searches the ray for crossings that are not adiabatic, which it would bridge with the Magnus patch of Adiabatic + Magnus hybrid strategy, and combines the eigenbases at the two ends. On the Sun, it finds no such crossing for any of the four Hamiltonians, so the result depends only on those two eigenbases. The neutrino starts in a flavor state, the default of average_initial_state, so the result is the averaged limit on a varying profile. Its interference terms carry the phases above, which damp them to nothing, so the result equals the averaged limit on a varying profile to within \(5 \cdot 10^{-15}\).
In Listing Averaged solar probabilities, only the wrapper changes between the four probability calls, and with it the Hamiltonian. The cost of a call follows the number of pairs of levels that the engine tests for adiabaticity along the ray, not the content of the Hamiltonian: three pairs at three flavors, six at \(3+1\), and ten at \(3+2\). The non-standard interactions change the matter term but add no levels, so they cost the same as the standard case. Sterile states add levels, and so do pseudo-Dirac partners, one per split pair (Shipped Hamiltonians). The four curves take about 0.5 s, 0.5 s, 0.8 s, and 1.2 s, respectively, over ninety energies.
Averaged solar probabilities. The four averaged solar probabilities of Figure The averaged solar survival probability, complete and runnable as it stands: the BS2005-AGS,OP table [BSB05] ships with Magνs, and every Sun wrapper takes it by name through density_profile. The same keywords serve the four wrappers. The keyword average=True is the only difference from an instantaneous call: it returns the averaged limit on a varying profile from the eigenbases at the two ends of the ray, without propagating the phase. The four curves take about three seconds. Naming 'B16-GS98' instead gives the B16-GS98 curve of Figure Three solar models compared. See notebooks #10, #13, and #28.
import numpy as np
import magnus.oscprob as oscprob
import magnus.globaldefs as gd
import magnus.solarmodels as solarmodels
# The table ships with Magnus. Its columns are read
# only to plot the density (top panel of the figure)
model = 'BS05-AGS-OP'
r = solarmodels.load_solar_model(model)['r_over_r_sun']
ne_sun = solarmodels.electron_density_profile(model)
ne = ne_sun(r*gd.SUN_RADIUS*gd.UNIT_KM)
# The solar spectrum, 0.1 to 20 MeV
E = np.logspace(-1.0, np.log10(20.0), 90)*gd.UNIT_MEV
osc = gd.load_nufit_params('NuFIT 6.1')
# Shared by the four calls: the model by name, from
# the center to its last row, nu_e to nu_e, the average
R = solarmodels.table_edge(model)
run = dict(nu_i=gd.NUE, nu_f=gd.NUE, average=True,
density_profile=model)
# Standard three-flavor oscillations
P3 = oscprob.osc_prob_3nu_sun(E, R, 0.0, **osc, **run)
# The same, with non-standard interactions
eps = dict(eps_ee=0.10, eps_em=0.05, eps_mt=0.03)
P3_nsi = oscprob.osc_prob_3nu_sun_nsi(
E, R, 0.0, **osc, **eps, **run)
# One and two sterile states. n_n/n_p is read from
# the table; passing ratio_number_neutrons_to_protons
# overrides it
ster = dict(s14=0.10**0.5, s24=0.10**0.5, D41=1.0)
P4 = oscprob.osc_prob_4nu_sun(
E, R, 0.0, **osc, **ster, **run)
P5 = oscprob.osc_prob_5nu_sun(
E, R, 0.0, **osc, **ster, s15=0.06**0.5,
s25=0.06**0.5, D51=1.7, **run)
# P3 = 0.5449 at 0.1 MeV, 0.2993 at 20 MeV
Without the average keyword, the same call returns the instantaneous probability at the surface. That answer is correct, even if it is not what Figure The averaged solar survival probability shows. With an oscillation length of about 160 km, the probability swings across most of the interval [0,1] within a hundred kilometers of the surface:
dL = np.array([0, 50, 100])*gd.UNIT_KM
P = oscprob.osc_prob_3nu_sun(
5.0*gd.UNIT_MEV, R - dL, 0.0, **osc,
nu_i=gd.NUE, nu_f=gd.NUE,
density_profile=model)
# P = 0.666, 0.092, 0.259 in about 7 s
These values are not numerical noise. The hybrid engine of Adiabatic + Magnus hybrid strategy computes and certifies them: two successive passes, on finer grids, agree to the requested tolerance, \(10^{-3}\) by default.
Varying the production point
Dependence on the production point in the Sun. Phase-averaged survival probability of a solar neutrino as a function of where it is produced, along the radial ray of Figure The averaged solar survival probability, at three energies, with and without the non-standard interactions (NSI) of Figure The averaged solar survival probability. The density is from the BS2005-AGS,OP solar model [BSB05]. The band marks where \(90\%\) of the \(^8\)B neutrinos are made, in the same model. See notebooks #10, #13, and #28.
The Sun placed the neutrino at the center of the Sun. Figure Dependence on the production point in the Sun moves its production point instead, \(r_0\), at fixed energy. At 20 MeV, a neutrino born at the center leaves as \(\nu_2\), with an average of 0.30. The farther out it is born, the lower the density at production, and the closer the average comes to the vacuum value, 0.55, which it nearly reaches by \(R_\odot/2\). The rise passes through the radius where the resonance density of Figure The averaged solar survival probability meets the profile, about \(R_\odot/3\); a neutrino born there has an average of \(\cos^4\theta_{13}/2 \approx 0.48\). As the energy drops, the rise moves inward, to about \(R_\odot/6\) at 5 MeV. At 1 MeV, the resonance density lies above the central density, and the curve rises only mildly, from 0.51 to 0.55. Each point of the figure is one call with average=True, as in Listing Averaged solar probabilities, with a different L0:
P = oscprob.osc_prob_3nu_sun(
5.0*gd.UNIT_MEV, R,
0.05*gd.SUN_RADIUS*gd.UNIT_KM,
**osc, **run)
# P = 0.389, against 0.375 from the center
At 5 MeV, moving the production point from the center to \(0.05\,R_\odot\), in the middle of where the \(^8\)B neutrinos are made, raises the average by 0.014.
The \(^8\)B neutrinos are made between \(0.02\,R_\odot\) and \(0.09\,R_\odot\) (the band in Figure Dependence on the production point in the Sun). The density is lowest at the outer edge of the band, and even there it exceeds the resonance density of a neutrino with more than 2.6 MeV. So a \(^8\)B neutrino with more than 2.6 MeV is born at a density above its resonance density, wherever in the band it is made, just like one born at the center. That is why Figure The averaged solar survival probability, computed for a neutrino born at the center, holds for the \(^8\)B neutrinos at high energy: at 20 MeV, the average is 0.301 from \(0.05\,R_\odot\), against 0.299 from the center. The non-standard interactions shift the curves by up to 0.012.
Probability averaging
In principle, the averaged probability could be estimated from the instantaneous one, as its mean over a window of baselines at the end of the ray. On the Sun, that estimate is close, but biased. The window must span many oscillation lengths for the mean to wash out the oscillation, yet across such a window the density changes, and so does the averaged probability, which depends on the density at the end of the ray. In notebook #13, for a 5-MeV neutrino at two flavors, on a ray from the center through the inner tenth of the solar radius, the mean over the last 6, 12, 24, and 48 oscillation lengths lies between 0.599 and 0.601, while the true averaged probability at the end of the ray is 0.595. The bias is small, but it grows as the window widens, from below \(10^{-6}\) over 6 oscillation lengths to 0.009 over 48.
A scan of the instantaneous probability along the ray has a second problem: aliasing. A scan needs at least two points per cycle of the oscillation it samples. With fewer, each value is still correct, but a curve drawn through them is not the probability. On the Sun, the fastest oscillation of a 5-MeV neutrino has a length of 4.9 km, so a scan of the last 100 km in steps of 50 km is aliased. Passing strategy_info reports it:
dL = np.array([0, 50, 100])*gd.UNIT_KM
info = {}
P = oscprob.osc_prob_3nu_sun(
5.0*gd.UNIT_MEV, R - dL, 0.0, **osc,
nu_i=gd.NUE, nu_f=gd.NUE,
density_profile=model,
strategy_info=info)
s = info['sampling']
s['cycles_per_step'] # 10.1
s['aliased'] # True
s['nyquist_points'] # 277028
Ten cycles fit between neighboring points. To resolve every cycle from the center to the surface, a scan would need \(2.8 \cdot 10^5\) points. Magνs reports aliasing, but does not warn about it, since nearly every realistic scan of a solar ray is aliased.
The keyword average=True avoids both problems: it takes no window and runs no scan. It evaluates the limit directly from the eigenbases at the two ends of the ray (Phase-averaged probabilities). At two flavors, on an adiabatic passage, the averaged survival probability of a neutrino started decohered has the textbook form
with \(\theta_m\) the mixing angle in matter at the two ends of the ray. Notebook #13, with average_initial_state='decohered', reproduces this expression to within \(3 \cdot 10^{-16}\) across 1–20 MeV.
Comparing solar models
|
Solar density profile |
|---|---|
|
Exponential fit [GK07b]; the default |
|
BP2000 [BPB01]; to \(0.95\,R_\odot\) |
|
BP04 [BP04]; to \(0.95\,R_\odot\) |
|
BS2005, GS98 composition [BSB05]; to \(0.98\,R_\odot\) |
|
BS2005, AGS05 composition [BSB05]; to \(0.98\,R_\odot\) |
|
B16, GS98 composition [VSV+17] |
|
B16, AGSS09met composition [VSV+17] |
|
B23, GS98 composition [HS23] |
|
B23, AGSS09 composition [HS23] |
|
B23, C11 composition [HS23] |
|
B23, AAG21 composition [HS23] |
|
B23, MB22 meteoritic composition [HS23] |
|
B23, MB22 photospheric composition [HS23] |
The solar density profiles that every Sun wrapper takes through density_profile. The tabulated standard solar models ship with Magνs. The B16 and B23 tables reach the surface. The solar compositions are GS98 [GS98], AGS05 [AGS05], AGSS09 [AGSS09], C11 [CLS+11], AAG21 [AAG21], and MB22 [M+22]. See notebook #13.
Three solar models compared. Top: the tabulated BS2005-AGS,OP [BSB05] and B16-GS98 [VSV+17] standard solar models, and the exponential fit that the Sun wrappers use by default, \(n_e = 245\,N_A\,e^{-10.54\,r/R_\odot}\) cm\(^{-3}\) [GK07b], with \(N_A\) the Avogadro number; underneath, the ratio of each to BS2005-AGS,OP. Bottom: \(\langle P_{\nu_e \to \nu_e}\rangle\) at three flavors on each profile; underneath, the difference from the BS2005-AGS,OP result. See notebooks #10, #13, and #28.
Table Solar density profiles lists the solar density profiles that every Sun wrapper accepts through density_profile. Besides an exponential profile, which is the default, Magνs ships twelve tabulated standard solar models from five series published between 2001 and 2023, several of them in more than one solar composition. Each is selected by name, which is not case-sensitive, e.g., density_profile='B16-GS98'.
From the tabulated mass density, \(\rho\), and hydrogen mass fraction, \(X\), the wrapper builds the electron number density, \(n_e = \rho N_A\,[X/m_{\rm H} + 2(1 - X)/m_{\rm He}]\), with \(m_{\rm H}\) and \(m_{\rm He}\) the atomic masses of hydrogen and helium-4 (Standard solar models), and interpolates it log-linearly in radius. Below the first tabulated radius, it holds the density at its first value, since the core is flat. Past the last one, which lies at \(0.95\,R_\odot\) for BP2000 and BP04 and at \(0.98\,R_\odot\) for BS2005, it continues the density along the logarithmic slope of the last tabulated interval. Passing stop_at_table_edge=True returns NaN there instead, with a warning. With sterile states, the wrapper reads the neutron-to-proton ratio automatically from the same table, with the same atomic masses, unless ratio_number_neutrons_to_protons is passed. In notebook #13, we compare all twelve models: they give the same averaged \(\nu_e\) survival probability to within \(2.4 \cdot 10^{-3}\). Most of that spread is between series of models; the six compositions of the B23 series differ by at most \(2 \cdot 10^{-4}\).
Figure Three solar models compared compares the probabilities computed with three choices of solar density profile. Two are tabulated standard solar models: BS2005-AGS,OP [BSB05], used throughout this section, and the newer B16-GS98 [VSV+17], with updated nuclear rates, equation of state, and opacities. The third is the exponential fit that textbooks quote [GK07b] and that the Sun wrappers use by default. The densities of the two solar models are within 2% of each other over the inner half of the Sun and 9% over the outer half. The exponential fit is 2.4 times too dense at the center, where the real profile flattens, and several times too dense again near the surface, where the real profile plunges.
The two tabulated models give the same averaged probability to within \(1.3 \cdot 10^{-3}\) at every energy, far below the precision of any solar measurement. The exponential fit gives a probability lower by up to 0.1, at 2.5 MeV. The reason is in The Sun: on an adiabatic passage, the averaged probability depends only on the eigenbases at the two ends of the ray. At the surface, the eigenbasis is that of vacuum, whatever the profile. At the center, it is set by the central density, where the fit is 2.4 times too dense. The probability falls from 0.55 at low energy to 0.30 at high energy, and the fall happens near the resonance energy, which is inversely proportional to the density. Therefore, with the fit, the fall happens at energies about 2.4 times lower. The values at low and high energy are nearly the same for every profile, so the difference is largest in between, at a few MeV.
Comparing against approximations
Error of the two-flavor solar approximation. Relative error of the two-flavor approximation, \(\langle P^{3\nu}_{\nu_e \to \nu_e} \rangle_{\rm approx}\) of Eq. (2), to the averaged solar survival probability \(\langle P^{3\nu}_{\nu_e \to \nu_e} \rangle\) computed with Magνs, on the three solar density profiles of Figure Three solar models compared. The approximation overestimates the probability at every energy, by an amount that grows in proportion to the energy. On the two tabulated models, it reaches 0.6% at 20 MeV. On the exponential fit, whose central density is 2.4 times higher, it is 2.3 to 3.4 times larger, and its curve leaves the panel above 12 MeV. See notebook #28.
Solar analyses rarely solve the three-flavor problem. The standard approximation [KP89] reduces it to two flavors in two steps. First, the 1–2 sector is solved on its own, as a two-flavor problem with \(\theta_{12}\) and \(\Delta m^2_{21}\), on the electron density multiplied by \(\cos^2\theta_{13}\). Second, \(\theta_{13}\) is folded back in through
Here, \(\langle P^{2\nu}_{\nu_e \to \nu_e} \rangle_{\rm approx}\) is Eq. (1), with \(\theta_m\) computed from \(\theta_{12}\), \(\Delta m^2_{21}\), and the electron density multiplied by \(\cos^2\theta_{13}\). The approximation rests on \(\Delta m^2_{31}\) being large compared with \(\Delta m^2_{21}\) and with the matter potential, so that the third mass state decouples. In vacuum, the approximation is exact. In matter, it drops terms of relative size \(2 E V_{\rm CC}/\Delta m^2_{31}\), about 0.1 at the center of the Sun for a 20-MeV neutrino, multiplied by the admixture of the decoupled state, \(\sin^2\theta_{13} \approx 0.02\).
Figure Error of the two-flavor solar approximation evaluates them on the three density profiles of Figure Three solar models compared. The three-flavor result, \(\langle P^{3\nu}_{\nu_e \to \nu_e} \rangle\), is computed with Magνs: P3 of Listing Averaged solar probabilities for the BS2005-AGS,OP table, and the same call on the other two profiles. The two-flavor result needs the density multiplied by \(\cos^2\theta_{13}\), which a wrapper cannot provide, since it takes the solar model by name. So it goes through the scenario function, which accepts any profile, here ne_sun of Listing Averaged solar probabilities multiplied by \(\cos^2\theta_{13}\):
# The 1-2 sector on cos^2 th13 * n_e
c13sq = 1.0 - osc['s13']**2 # cos^2 th13
P2 = oscprob.osc_prob_matter_std_potential(
2, lambda l: c13sq*ne_sun(l), E, R,
dict(sth=osc['s12'], Dm2=osc['D21']),
L0=0.0, nu_i=gd.NUE, nu_f=gd.NUE, average=True,
density_is_of_number_of_electrons=True)
# Fold theta_13 back in
P_approx = osc['s13']**4 + c13sq**2*P2
(P_approx - P3)/P3 # 6e-3 at 20 MeV
The P2 equals Eq. (1) to within \(10^{-16}\) at every energy.
The approximation overestimates the probability at every energy, on every profile, and its relative error grows in proportion to the energy, as the size of the dropped terms does. On the two tabulated models, which give nearly the same probability (Comparing solar models), the error is about \(3 \cdot 10^{-5}\) at 0.1 MeV, \(3 \cdot 10^{-4}\) at 1 MeV, and 0.6% at 20 MeV. On the exponential fit, it is 2.3 to 3.4 times larger, 2% at 20 MeV, because the fit’s central density is 2.4 times higher and the dropped terms grow with it.