The Hamiltonians
A direct call to osc_prob takes the Hamiltonian as a matrix, or as a function of position that returns one (Using wrappers vs. scenario calls vs. direct calls). Magνs ships builders for standard and non-standard terms (Shipped Hamiltonians), and any other term can be written by hand (Building a new Hamiltonian from scratch).
Returns |
|
|---|---|
Vacuum term |
|
Same, without the \(1/E\) |
|
Vacuum term at a position |
|
Both of the above |
|
\(V_{\rm CC}\) times the projector |
|
Same, from a potential function |
|
Non-standard term alone |
|
Same, from a potential function |
|
Lorentz-violating term alone |
|
Same, without the energy factor |
|
Pseudo-Dirac vacuum term |
|
Same, without the \(1/E\) |
|
Its matter term |
|
The \(43\) Hamiltonian builders of magnus.hamiltonians, each named with the prefix hamiltonian_. The ten variants above the table break exist at two, three, four, and five flavors. The three below the break build a pseudo-Dirac spectrum from the active one, and take the pairs of states to split, so their names carry no number of flavors.
Shipped Hamiltonians
Table The Hamiltonian builders lists the 43 Hamiltonian builders shipped in magnus.hamiltonians: vacuum, matter, non-standard interactions, and Lorentz violation at two to five flavors, plus three for a pseudo-Dirac spectrum.
A Hamiltonian is a sum of terms, as in the matter Hamiltonian. Only the vacuum builders return a complete Hamiltonian; every other builder in Table The Hamiltonian builders returns its own term alone. For instance, at three flavors, the matter builder returns \(V_{\rm CC}\) times the projector, which is zero except in the electron entry, [0][0]. Added to a vacuum term, it gives the matter Hamiltonian. Passed to osc_prob on its own, it returns the probabilities as the identity matrix, since a diagonal Hamiltonian does not mix flavors.
Every builder returns a \(d \times d\) NumPy array, with \(d\) the number of flavors, so terms add with + and scale with *. The builders take their parameters under the same names as the wrappers, so a three-flavor set loaded with load_nufit_params can be passed to them with **. Every physical parameter must be given, since a builder has no default set to fall back on. The one exception is the neutron-to-proton ratio of the four- and five-flavor matter builders, which is 1 by default.
The snippets on this page reuse E, L, Es, and osc from Listing Vacuum and constant density, three ways, and import the builders as ham.
Vacuum.—At three flavors, the vacuum Hamiltonian is
import magnus.hamiltonians as ham
H = ham.hamiltonian_3nu_vacuum(E, **osc)
H.shape # (3, 3)
The call returns the vacuum Hamiltonian at three flavors, i.e.,
with \(\mathbb{R}\) the PMNS matrix (The Hamiltonian and the probability), built from \(\theta_{12}\), \(\theta_{23}\), \(\theta_{13}\), and \(\delta_{\rm CP}\). At two flavors, the builder takes sth and Dm2, and returns the vacuum Hamiltonian with \(\mathbb{R} = \mathbb{R}_{12}\) and \(\mathbb{M}^2 = {\rm diag}(0, \Delta m^2)\).
At four and five flavors, the builder also takes the mixing angles, phases, and mass-squared splitting of each sterile state, under the same names as in the wrappers:
st4 = dict(s14=0.1, s24=0.1, s34=0.1,
d14=0.0, d24=0.0, D41=1.0)
H = ham.hamiltonian_4nu_vacuum(
E, **osc, **st4)
H.shape # (4, 4)
st5 = dict(st4, s15=0.1, s25=0.1, s35=0.1,
d15=0.0, d35=0.0, D51=1.7)
H = ham.hamiltonian_5nu_vacuum(
E, **osc, **st5)
H.shape # (5, 5)
The two calls return the vacuum Hamiltonian at four and five flavors, i.e.,
with \(\mathbb{R}\) the mixing matrix at \(n = 4\) and 5, respectively.
A builder marked _energy_independent returns the same matrix without its \(1/E\) factor. At one energy, the two builders are interchangeable. Across a scan in energy, the energy-independent one is faster, because the mixing rotation does not depend on the energy and is then computed only once. Over the 200 energies Es of Listing Vacuum and constant density, three ways, rebuilding the Hamiltonian at each energy takes about 4 ms, and building it once and dividing it by each energy takes 0.2 ms:
# Interchangeable at one energy
H1 = ham.hamiltonian_3nu_vacuum(E, **osc)
H0 = ham.\
hamiltonian_3nu_vacuum_energy_independent(
**osc)
H1 - H0/E # All zeros
# Over 200 energies: 4.1 ms
Hs = [ham.hamiltonian_3nu_vacuum(e, **osc)
for e in Es]
# ... and 0.16 ms
Hs = [H0/e for e in Es]
Matter.—At three flavors, the matter term is
V = gd.VCC_EARTH_CRUST # 1.145e-13 eV
H = ham.hamiltonian_3nu_matter(V)
H.shape # (3, 3)
The call returns the matter term of the matter Hamiltonian at three flavors, i.e.,
the charged-current potential acting on the electron flavor alone. (VCC_EARTH_CRUST is \(V_{\rm CC}\) at density 3 g cm\(^{-3}\) with \(Y_e = 0.5\).) At two flavors, the call is the same; the projector is \({\rm diag}(1, 0)\).
At four and five flavors, the projector also has an entry for each sterile flavor, \(r/2\), where \(r\) is the neutron-to-proton ratio of the medium (Conventions). The builder takes \(r\) as a second parameter, 1 by default. In neutral matter, \(r = (1 - Y_e)/Y_e\) (Constant density), so \(r\) follows from the electron fraction; for instance, from that of the core of the Earth, which Magνs stores in magnus.earth:
import magnus.earth as earth
Ye = earth.Y_E_CORE_PREM # 0.4656
r = (1.0 - Ye)/Ye # 1.15
H = ham.hamiltonian_4nu_matter(V, r)
H.shape # (4, 4)
H = ham.hamiltonian_5nu_matter(V, r)
H.shape # (5, 5)
The builder takes the potential and the ratio independently. The two calls return
Four of the ten variants in Table The Hamiltonian builders end in _td (for time-dependent): the two vacuum forms, matter, and nsi. Each takes the position along the path, \(l\), as its first argument, just as the function of position passed to osc_prob does. The vacuum ones ignore it and return the same matrix at every position. matter_td and nsi_td take the potential as a function of position, f, rather than as a number, and evaluate it at \(l\); e.g., matter_td(l, f) returns the same matrix as matter(f(l)).
Non-standard interactions.—At three flavors, the NSI term is built as below. This block and the other NSI and LIV blocks show the signatures only: each argument stands for a number.
H = ham.hamiltonian_3nu_nsi(
V, eps_ee, eps_em, eps_et,
eps_mm, eps_mt, eps_tt)
H.shape # (3, 3)
The call returns the charged-current potential times the matrix of NSI couplings, i.e.,
The six couplings passed are the diagonal entries and those above it. The builder fills the entries below the diagonal with their complex conjugates, so the matrix is Hermitian even when the off-diagonal couplings are complex. At two flavors, the builder takes two couplings, eps_aa and eps_ab, and returns
where \(a\) and \(b\) are the two flavors, \(e\) and \(\mu\) by default (Two flavors). The lower diagonal entry is zero because only the difference between the two diagonal entries affects probabilities. Adding the same amount to both shifts every energy level of the Hamiltonian equally, which changes only an overall phase of the evolution, so \(\varepsilon_{aa}\) stands for the difference \(\varepsilon_{aa} - \varepsilon_{bb}\).
At four and five flavors, the matrix of couplings also has rows and columns for the sterile flavors, and the builder takes the entries on and above its diagonal, ten and fifteen of them:
H = ham.hamiltonian_4nu_nsi(
V, eps_ee, eps_em, eps_et, eps_es,
eps_mm, eps_mt, eps_ms, eps_tt,
eps_ts, eps_ss)
H.shape # (4, 4)
H = ham.hamiltonian_5nu_nsi(
V, eps_ee, eps_em, eps_et, eps_es1,
eps_es2, eps_mm, eps_mt, eps_ms1,
eps_ms2, eps_tt, eps_ts1, eps_ts2,
eps_s1s1, eps_s1s2, eps_s2s2)
H.shape # (5, 5)
The two calls return the matrix whose entries are
with \(\alpha, \beta \in \{e, \mu, \tau, s\}\) at four flavors and \(\alpha, \beta \in \{e, \mu, \tau, s_1, s_2\}\) at five; Eq. (5) is the same expression at three flavors.
Lorentz-invariance violation.—At three flavors, the Lorentz-invariance-violating (LIV) term is
H = ham.hamiltonian_3nu_liv(
E, sxi12, sxi23, sxi13, dxiCP,
b1, b2, b3, Lambda, n_liv)
H.shape # (3, 3)
The call returns
Here, \(b_1\), \(b_2\), and \(b_3\) are the eigenvalues of the LIV operator, in eV. The unitary matrix \(\mathbb{V}_\xi\) rotates them into the flavor basis; it is built from three angles, \(\xi_{12}\), \(\xi_{23}\), and \(\xi_{13}\), and one phase, \(\delta_\xi\), in the same way as the PMNS matrix (The Hamiltonian and the probability). The prefactor sets how the term grows with energy: \(\Lambda\) is the energy scale of the operator, in eV, and \(n_{\rm LIV}\) is its dimension minus three, so that \(n_{\rm LIV} = 0\) gives a term that does not depend on the energy.
The angles \(\xi_{ij}\) and the phase \(\delta_\xi\) are parameters of the LIV operator alone, independent of the PMNS angles and phase. They fix the basis in which the operator is diagonal, which need be neither the flavor basis nor the mass basis. For instance, with every \(\xi_{ij}\) and \(\delta_\xi\) set to zero, \(\mathbb{V}_\xi\) is the identity and the operator is diagonal in the flavor basis; with them set equal to the PMNS values, it is diagonal in the mass basis.
At four and five flavors, the operator has one eigenvalue per flavor, and \(\mathbb{V}_\xi\) is built as the mixing matrix \(\mathbb{R}\) at \(n = 4\) and \(n = 5\), with angles \(\xi_{ij}\) in place of \(\theta_{ij}\) and phases \(\delta_{\xi, ij}\) in place of \(\delta_{ij}\):
H = ham.hamiltonian_4nu_liv(
E, sxi12, sxi23, sxi13, dxiCP,
sxi14, dxi14, sxi24, dxi24, sxi34,
b1, b2, b3, b4, Lambda, n_liv)
H.shape # (4, 4)
H = ham.hamiltonian_5nu_liv(
E, sxi12, sxi23, sxi13, dxiCP,
sxi14, dxi14, sxi15, dxi15, sxi24,
dxi24, sxi25, sxi34, sxi35, dxi35,
b1, b2, b3, b4, b5, Lambda, n_liv)
H.shape # (5, 5)
The two calls return
At two flavors, \(\mathbb{V}_\xi\) is a single rotation, by the angle \(\xi\), and the operator has two eigenvalues.
The builders marked _energy_independent return these terms without the factor \(E^{n_{\rm LIV}}\), which you then multiply back in at each energy, as for the vacuum term.
Pseudo-Dirac.—A pseudo-Dirac neutrino is a mass state split into a nearly degenerate pair, one member active and the other sterile. The pseudo-Dirac builders take the mixing matrix and the mass-squared values of the active sector, and the mass states to split:
R = ham.pmns_mixing_matrix(
osc['s12'], osc['s23'], osc['s13'],
osc['dCP'])
H = ham.hamiltonian_pseudo_dirac_vacuum(
E, R, [0.0, osc['D21'], osc['D31']],
{1: 1.0e-18})
H.shape # (4, 4)
The last argument lists the mass states to split, each with the mass-squared splitting of its pair, \(\delta m^2_j\), in eV\(^2\). The mass states are counted from 0, so {1: 1.0e-18} splits \(\nu_2\) by \(\delta m^2_2 = 10^{-18}\) eV\(^2\). States not listed stay unsplit; with an empty mapping, the builder returns the three-flavor vacuum Hamiltonian, to round-off.
The call returns the vacuum Hamiltonian for the enlarged set of states, i.e.,
Each split state \(\nu_j\) is replaced by two mass states, \((\nu_j \pm s_j)/\sqrt{2}\), where \(s_j\) is its sterile partner, with masses squared \(m_j^2\) and \(m_j^2 + \delta m^2_j\). For the example above, which splits \(\nu_2\), the rows of \(\mathbb{R}_{\rm PD}\) are the flavors \(e\), \(\mu\), \(\tau\), and \(s_2\), and its columns are \(\nu_1\), the two members of the pair, and \(\nu_3\):
with \(\mathbb{R}\) the PMNS matrix, and
In general, there is one state per active state, plus one per split state. Both matrices are available on their own, from pseudo_dirac_mixing_matrix and pseudo_dirac_mass_squared.
The matter builder, hamiltonian_pseudo_dirac_matter, takes the potential, the number of active states, the same mapping of split states, and the neutron-to-proton ratio, \(r\), 1 by default. It returns \(V_{\rm CC}\) times the projector in the same enlarged flavor basis. For the example above, which splits \(\nu_2\),
with rows \(e\), \(\mu\), \(\tau\), and \(s_2\). Since each partner \(s_j\) is sterile, its entry is \(r/2\), as for the sterile states of Eq. (4). Split states other than \(\nu_2\) add one such entry each.
Building a new Hamiltonian from scratch
There are two ways to build a Hamiltonian that Magνs does not ship. The first is to add shipped terms, as Using wrappers vs. scenario calls vs. direct calls did with a vacuum, a matter, and a non-standard term. With a _td builder, a term can also vary along the path. For instance, a matter potential that grows linearly along the baseline is
H_vac = ham.hamiltonian_3nu_vacuum(E, **osc)
VCC_func = lambda l: V*(1.0 + l/L)
def H_func(l):
H_mat = ham.hamiltonian_3nu_matter_td(
l, VCC_func)
return H_vac + H_mat
P = oscprob.osc_prob(H_func, 0.0, L,
rtol=1e-8, atol=1e-8)
# Pme = 0.021309
The second way is to write the matrix directly. osc_prob requires only a function that returns a Hermitian matrix, of any size, so an entirely new term reaches the solver in the same way as a shipped one. For instance, an off-diagonal term between \(\nu_e\) and \(\nu_\tau\) that grows along the path is
def H_own(l):
H = H_vac.copy()
H[0, 2] += 1.0e-13*(l/L)
H[2, 0] += 1.0e-13*(l/L)
return H
P = oscprob.osc_prob(H_own, 0.0, L,
rtol=1e-8, atol=1e-8)
# Pme = 0.004865
The Hermiticity of the matrix is yours to guarantee, since it is what keeps the evolution unitary (Against other codes). Magνs neither checks nor enforces it, and a matrix that is not Hermitian returns probabilities that do not add up to one. Every entry of the matrix is in eV, and every position in eV\(^{-1}\) (Conventions). Long-range interactions in the Sun shows a complete example.
A custom Hamiltonian may have more than five flavors. The wrappers of Table The probability wrapper functions, the builders of Table The Hamiltonian builders, and the compiled exponential kernels of Methodology stop at five. Above five, you supply the matrix, either to osc_prob or, as a vacuum term, to a scenario function through h_vac_energy_indep (Using wrappers vs. scenario calls vs. direct calls), and the exponential of each slab goes through a general Hermitian eigendecomposition, numpy.linalg.eigh. Each slab then costs more: on one exponential density profile, at the default tolerance, a probability takes about twice as long at six flavors as at five, and three and a half times as long at eight. The probabilities remain unitary to within \(10^{-13}\).