Standard solar models
The Sun entry points — osc_prob_{2,3,4,5}nu_sun[_nsi|_liv] and
osc_prob_sun() — describe the Sun by an exponential fit to its
electron density by default. density_profile puts a published standard
solar model in its place, from tables that ship with the package.
Added in version 1.2.0.
|
Model |
Tabulated range [\(R_\odot\)] |
Reference |
|---|---|---|---|
|
Bahcall, Pinsonneault & Basu |
0.0065–0.949 |
[BPB01] |
|
Bahcall & Pinsonneault |
0.0065–0.947 |
[BP04] |
|
Bahcall, Serenelli & Basu, GS98 composition |
0.0016–0.983 |
[BSB05] |
|
Bahcall, Serenelli & Basu, AGS05 composition |
0.0016–0.983 |
[BSB05] |
|
Vinyoles et al., GS98 composition |
0.0005–1 |
[VSV+17] |
|
Vinyoles et al., AGSS09met composition |
0.0005–1 |
[VSV+17] |
|
Herrera & Serenelli, one per solar composition |
0–1 |
[HS23] |
Names are case-insensitive. The default is 'exp'.
import magnus.globaldefs as gd
from magnus import oscprob
R_sun = gd.SUN_RADIUS*gd.UNIT_KM
for profile in ('exp', 'BP04', 'B16-GS98', 'B23-AAG21'):
P_ee = oscprob.osc_prob_3nu_sun(energy=10*gd.UNIT_MEV, L=R_sun, L0=0.0,
nu_i=gd.NUE, nu_f=gd.NUE, average=True,
density_profile=profile)
print('%-10s <P_ee> = %.4f' % (profile, P_ee))
# exp 0.2971, BP04 0.3183, B16-GS98 0.3187, B23-AAG21 0.3191
exp <P_ee> = 0.2971
BP04 <P_ee> = 0.3183
B16-GS98 <P_ee> = 0.3187
B23-AAG21 <P_ee> = 0.3191
L0 is the production radius and must be a nonnegative scalar. L may be an array
of end points, but L0 may not; for several production points, make one call per point.
From the command line, magnus --environment sun takes the same names through
--density-profile (see Command-line calculator). Notebook 13
compares all twelve on the averaged observable: they agree to \(2.4\times10^{-3}\),
and the exponential fit is off by up to 0.1, at 2.5 MeV (0.02 at 10 MeV, as above).
Top: the tabulated 'BS05-AGS-OP' and 'B16-GS98' models and the exponential fit
the Sun functions use by default, with the ratio of each to 'BS05-AGS-OP'. Bottom:
\(\langle P_{\nu_e \to \nu_e}\rangle\) at three flavors on each profile, and
its difference from the 'BS05-AGS-OP' result. The two tables agree to
\(1.3 \times 10^{-3}\); the fit is low by up to 0.1. From the Magνs paper.
What is tabulated, and what is computed from it
Each table holds three columns of the authors’ file, copied as written: the radius, the mass density \(\rho\) and the hydrogen mass fraction \(X\). The electron density is
with \(m_{\rm H} = 1.00783\) and \(m_{\rm He} = 4.00260\) the atomic masses of
\(^1\)H and \(^4\)He in atomic mass units: everything heavier than hydrogen is
counted as helium-4. Against the full composition of 'BS05-OP' (\(^1\)H,
\(^4\)He, \(^3\)He, \(^{12}\)C, \(^{14}\)N, \(^{16}\)O) this is
within 0.002% at the center and 0.07% at worst, near \(0.3\,R_\odot\) where
\(^3\)He peaks. The textbook \(\rho N_A (1 + X)/2\), which Bahcall’s own
electron-density tables follow, is 0.4 to 0.8% higher. Notebooks 13 and 28 compute
\(n_e\) with the same formula, and inside the table they reproduce
'BS05-AGS-OP' exactly. The density \(n_e\) is interpolated linearly in its logarithm.
For the wrappers with sterile states, the neutral-current term needs the neutron-to-proton ratio, helium’s neutrons over all the protons, from the same atomic masses:
With a solar model and
ratio_number_neutrons_to_protons left at its default, None, the wrappers take it
from the same table, radius by radius; a value or a callable passed explicitly is used
instead. The exponential fit carries no composition, so with it the ratio defaults to 1.0.
osc_prob_sun() passes only the electron density to H_func; a
Hamiltonian with sterile states that wants the composition can take it from
magnus.solarmodels.neutron_to_proton_ratio_profile().
Outside the table
Only the B23 tables start at the center, and the Bahcall tables (BP2000 to BS05) stop short of the surface. Below the first row, the density holds its first value, since the core is flat. Past the last row, it continues along the logarithmic slope of the last tabulated interval, so that it keeps falling rather than stopping at a constant or dropping to zero. That is a continuation, not a model of the outer Sun. Continued from 0.95 \(R_\odot\), BP2000 and BP04 reach about \(1.6\times10^{-3}\) g cm−3 at the surface, whereas B16 and B23, which are tabulated there, give about \(1.7\times10^{-7}\) g cm−3.
stop_at_table_edge=True refuses the continuation. A baseline that ends past the last
row returns NaN, with a SolarModelRangeWarning naming the edge,
and the others are computed as usual:
import warnings
import numpy as np
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter('always')
P_ee = oscprob.osc_prob_3nu_sun(energy=10*gd.UNIT_MEV,
L=np.array([0.5, 0.9, 1.0])*R_sun, L0=0.0,
nu_i=gd.NUE, nu_f=gd.NUE, average=True,
density_profile='BP04', stop_at_table_edge=True)
print(P_ee)
print(caught[0].category.__name__)
[0.32640971 0.31842657 nan]
SolarModelRangeWarning
With stop_at_table_edge=True, a path that starts past the last row has nothing inside the
table to compute, and is refused with a ValueError; without it, the profile is continued
and a number is returned. stop_at_table_edge with the exponential fit, which has no last
row, is refused as well. magnus.solarmodels.table_edge() gives the edge in the units of
L.
Coherent probabilities through a table
The interpolated profile has a kink at every row, and the Bahcall tables, printed to four
significant figures, step through the flat core. The phase-averaged probability
(average=True) is unaffected: on all twelve models at 1, 5, 10 and 20 MeV, from the
center to 0.9 \(R_\odot\), it issues no warning.
A coherent probability over most of the Sun can be affected. The hybrid engine may then
fail to certify it and pass it to the slab ladder, which reaches its slab cap and issues
ToleranceNotAchievedWarning. On the same grid, 4 of the 48
coherent calls did so (BP04 and BS05-OP at 10 and 20 MeV). Their answers were still
accurate: on eleven cases, those four among them, every default answer was within
\(10^{-3}\) of a reference converged to \(10^{-9}\), and the four that warned were
within \(3\times10^{-5}\). The warning means that Magνs could not verify the answer,
not that the answer is wrong.
To obtain a verified coherent probability, put t_breakpoints at the rows, so that
no slab straddles a kink, and start the slab ladder fine enough to resolve the
oscillation:
from magnus import solarmodels
rows = solarmodels.load_solar_model('BS05-AGS-OP')['r_over_r_sun']*R_sun
P_ee = oscprob.osc_prob_3nu_sun(energy=10*gd.UNIT_MEV, L=0.95*R_sun, L0=0.0,
nu_i=gd.NUE, nu_f=gd.NUE,
density_profile='BS05-AGS-OP', strategy='magnus',
t_breakpoints=rows, n_slabs=200_000,
max_n_slabs=10_000_000)
print('P_ee = %.9f' % P_ee) # 0.126766302
P_ee = 0.126766302
On the same eleven cases, this came within \(4\times10^{-9}\) of the reference, with no
ToleranceNotAchievedWarning. Two issued
MagnusConvergenceWarning, which reports a slab width rather than an
error (see Accuracy and diagnostics). The measurements are in
docs/dev/adversarial_batteries/solar_model_engines.py and solar_model_coherent.py, in a
source checkout (a PyPI install does not ship them).
Provenance
magnus.solarmodels.solar_model_info() gives each model’s reference, the source it was
taken from, the date, the SHA-256 of the original file and its terms of use. Each table, in
magnus/data/solar_models/, carries the same in its header, above the authors’ original
header. The Bahcall tables come from the IAS archive, whose terms ask only for a note on how
they are used. The B16 tables come from the Internet Archive’s copy of the authors’ page,
which is offline, and the B23 tables from their Zenodo release, under CC BY 4.0. If you use
a model, cite its paper.
tools/build_solar_model_tables.py rebuilds the tables from the original files, and with
--check verifies that the committed ones are exactly what those files give. It refuses
any original whose hash has changed.