Part I: Solar Interior | Chapter 4

Solar Neutrinos

The MSW effect, neutrino oscillations, and the resolution of the solar neutrino problem

4.1 The Solar Neutrino Problem

Derivation 1: Expected vs Observed Neutrino Fluxes

The solar neutrino problem was one of the most important puzzles in 20th-century physics. The Standard Solar Model predicts specific neutrino fluxes from each branch of the pp chain and CNO cycle. The total $\nu_e$ flux at Earth is:

Step 1. Each pp chain termination produces neutrinos. The pp neutrino flux is directly tied to the solar luminosity:

$$\Phi_{pp} = \frac{2 L_\odot}{4\pi d^2 \langle Q_{pp} \rangle} \approx 6.0 \times 10^{10} \text{ cm}^{-2}\text{s}^{-1}$$

where $\langle Q_{pp}\rangle \approx 13.4$ MeV is the average thermal energy release per neutrino (accounting for the two neutrinos per ${}^4\text{He}$ produced).

Step 2. The Davis Homestake experiment (1968-1998) used the reaction$\nu_e + {}^{37}\text{Cl} \to {}^{37}\text{Ar} + e^-$ with threshold 0.814 MeV:

$$\text{Observed: } 2.56 \pm 0.23 \text{ SNU} \quad \text{vs} \quad \text{SSM prediction: } 7.6 \pm 1.2 \text{ SNU}$$

where 1 SNU (Solar Neutrino Unit) = $10^{-36}$ captures per target atom per second.

Step 3. The deficit of a factor of ~3 persisted across all experiments (Homestake, Kamiokande, SAGE, GALLEX/GNO) for three decades. The resolution required new physics: neutrino oscillations.

4.2 Neutrino Oscillations in Vacuum

Derivation 2: Two-Flavor Vacuum Oscillation Probability

If neutrinos have mass, the flavor eigenstates ($\nu_e, \nu_\mu$) differ from the mass eigenstates ($\nu_1, \nu_2$), connected by a unitary mixing matrix.

Step 1. In the two-flavor approximation, the mixing is parameterized by a single angle $\theta$:

$$\begin{pmatrix} \nu_e \\ \nu_\mu \end{pmatrix} = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} \nu_1 \\ \nu_2 \end{pmatrix}$$

Step 2. A neutrino produced as $\nu_e$ at $t=0$ evolves as:

$$|\nu(t)\rangle = \cos\theta \, e^{-iE_1 t} |\nu_1\rangle + \sin\theta \, e^{-iE_2 t} |\nu_2\rangle$$

Step 3. For ultra-relativistic neutrinos, $E_i \approx p + m_i^2/(2p) \approx E + m_i^2/(2E)$. The survival probability is:

$$\boxed{P(\nu_e \to \nu_e) = 1 - \sin^2 2\theta \sin^2\!\left(\frac{\Delta m^2 L}{4E}\right)}$$

where $\Delta m^2 = m_2^2 - m_1^2$ and $L$ is the propagation distance. The oscillation length is:

$$L_{\text{osc}} = \frac{4\pi E}{\Delta m^2} \approx 2.48 \frac{E \text{ [MeV]}}{\Delta m^2 \text{ [eV}^2\text{]}} \text{ m}$$

4.3 The MSW Effect

Derivation 3: Matter-Enhanced Oscillations

The Mikheyev-Smirnov-Wolfenstein (MSW) effect describes how the effective mixing angle is modified by coherent forward scattering of $\nu_e$ off electrons in matter. This is the dominant effect for solar neutrinos above about 2 MeV.

Step 1. In matter with electron density $n_e$, the $\nu_e$acquires an additional potential energy from charged-current interactions with electrons:

$$V = \sqrt{2} G_F n_e$$

Step 2. The effective Hamiltonian in the flavor basis becomes:

$$H = \frac{\Delta m^2}{4E} \begin{pmatrix} -\cos 2\theta + A & \sin 2\theta \\ \sin 2\theta & \cos 2\theta - A \end{pmatrix}$$

where the matter parameter is $A = 2\sqrt{2} G_F n_e E / \Delta m^2$.

Step 3. The effective mixing angle in matter:

$$\boxed{\sin^2 2\theta_m = \frac{\sin^2 2\theta}{(\cos 2\theta - A)^2 + \sin^2 2\theta}}$$

Step 4. Resonance occurs when $A = \cos 2\theta$, i.e.,$n_e^{\text{res}} = \Delta m^2 \cos 2\theta / (2\sqrt{2} G_F E)$. At resonance,$\theta_m = \pi/4$ and mixing is maximal regardless of the vacuum angle.

For ${}^8\text{B}$ neutrinos ($E \sim 10$ MeV), the resonance occurs at $r \approx 0.2 R_\odot$. If the density change through the resonance is sufficiently adiabatic ($\gamma = \Delta m^2 \sin^2 2\theta / (2E \cos 2\theta |dn_e/dr|_{\text{res}}) \gg 1$), the neutrino undergoes a level crossing and emerges primarily as $\nu_2$, giving a survival probability $P_{ee} \approx \sin^2\theta \approx 0.31$.

4.4 Experimental Verification

Derivation 4: SNO Neutral Current Proof

The Sudbury Neutrino Observatory provided the definitive solution by measuring both the$\nu_e$ flux and the total neutrino flux independently.

Step 1. SNO used heavy water ($\text{D}_2\text{O}$) enabling three reactions:

$$\text{CC: } \nu_e + d \to p + p + e^- \quad \text{(sensitive only to } \nu_e\text{)}$$

$$\text{NC: } \nu_x + d \to p + n + \nu_x \quad \text{(sensitive to all flavors equally)}$$

$$\text{ES: } \nu_x + e^- \to \nu_x + e^- \quad \text{(mainly } \nu_e\text{, some } \nu_{\mu,\tau}\text{)}$$

Step 2. The results (2001-2002) showed:

• CC flux: $\Phi_{\text{CC}} = (1.76 \pm 0.11) \times 10^6$ cm$^{-2}$s$^{-1}$ (only $\nu_e$)
• NC flux: $\Phi_{\text{NC}} = (5.09 \pm 0.63) \times 10^6$ cm$^{-2}$s$^{-1}$ (all flavors)
• SSM prediction: $\Phi_{\text{SSM}} = (5.05 \pm 0.81) \times 10^6$ cm$^{-2}$s$^{-1}$

Step 3. The NC measurement matched the SSM perfectly, proving that the total neutrino flux was correct but $\nu_e$ had oscillated to other flavors during transit. The $\nu_\mu + \nu_\tau$ flux: $\Phi_{\mu\tau} = \Phi_{\text{NC}} - \Phi_{\text{CC}} = (3.41 \pm 0.64) \times 10^6$ cm$^{-2}$s$^{-1}$ was non-zero at $5.3\sigma$.

This earned the 2015 Nobel Prize in Physics for Arthur McDonald (SNO) and Takaaki Kajita (Super-K).

4.5 Energy-Dependent Survival Probability

Derivation 5: Adiabatic MSW Survival Probability

Step 1. For neutrinos produced deep in the solar core where the density is high ($A \gg 1$), the matter mixing angle $\theta_m \approx \pi/2$, so the neutrino starts as nearly pure $\nu_2^m$ (the heavier mass eigenstate in matter).

Step 2. If the density decreases adiabatically (the Landau-Zener condition is satisfied), the neutrino remains in $\nu_2$ as it exits the Sun. In vacuum, $\nu_2$has $\nu_e$ content $\sin^2\theta$, so:

$$\boxed{P_{ee}^{\text{MSW}} = \sin^2\theta_{12} \approx 0.31 \quad \text{(high-energy, adiabatic limit)}}$$

Step 3. For low-energy neutrinos ($E \lesssim 2$ MeV), the matter effect is negligible, and vacuum oscillations dominate. After averaging over the rapid oscillations:

$$P_{ee}^{\text{vac}} = 1 - \frac{1}{2}\sin^2 2\theta_{12} \approx 0.57$$

Step 4. The general expression including the crossing probability $P_c$ is:

$$P_{ee} = \frac{1}{2} + \left(\frac{1}{2} - P_c\right)\cos 2\theta_m^0 \cos 2\theta_{12}$$

where $\theta_m^0$ is the matter angle at the production point and$P_c = \exp(-\pi\gamma/2)$ is the Landau-Zener level crossing probability. For the measured parameters ($\Delta m_{21}^2 = 7.53 \times 10^{-5}$ eV$^2$,$\theta_{12} = 33.4^\circ$), the transition between vacuum and MSW regimes occurs around 2-5 MeV. Borexino has beautifully confirmed this energy-dependent survival probability by measuring pp, ${}^7\text{Be}$, pep, and ${}^8\text{B}$neutrino rates individually.

Numerical Simulation

This simulation shows the energy-dependent survival probability, matter potential in the Sun, the MSW level-crossing diagram, and a comparison of experimental results.

Solar Neutrinos: Survival Probability, MSW Effect, and Experimental Results

Python

script.py175 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# Solar Neutrino Physics
# ============================================================

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Solar Neutrino Physics', fontsize=16, fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

for ax in axes.flat:
    ax.set_facecolor('#0a0a0a')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')
    for spine in ax.spines.values():
        spine.set_color('#333')

# Parameters
theta12 = np.radians(33.44)
dm2_21 = 7.53e-5  # eV^2
G_F = 1.166e-5     # GeV^-2

# --- Panel 1: Survival probability vs energy ---
ax1 = axes[0, 0]
E_arr = np.logspace(-1, 1.3, 500)  # MeV

# Electron density at production point (center)
n_e_center = 6e25  # cm^-3 (approximate)

# Matter parameter A = 2*sqrt(2)*G_F*n_e*E / dm^2
# Convert: G_F in GeV^-2, n_e in cm^-3, E in MeV
# A = 2*sqrt(2) * 1.166e-5 * n_e * (E*1e-3) / dm2
# with proper unit conversion: A = 7.63e-14 * n_e[cm^-3] * E[MeV] / dm2[eV^2]
A = 7.63e-14 * n_e_center * E_arr / dm2_21

# Effective mixing angle at production
sin2_2theta_m = np.sin(2*theta12)**2 / ((np.cos(2*theta12) - A)**2 + np.sin(2*theta12)**2)
cos_2theta_m = (np.cos(2*theta12) - A) / np.sqrt((np.cos(2*theta12) - A)**2 + np.sin(2*theta12)**2)

# Adiabatic survival probability
P_ee = 0.5 + 0.5 * cos_2theta_m * np.cos(2*theta12)

# Vacuum average
P_vac = 1 - 0.5 * np.sin(2*theta12)**2

ax1.semilogx(E_arr, P_ee, color='#f43f5e', linewidth=2.5, label='MSW + vacuum')
ax1.axhline(y=np.sin(theta12)**2, color='#fb923c', linestyle='--', linewidth=1.5,
           label='MSW limit (sin^2 theta)')
ax1.axhline(y=P_vac, color='#60a5fa', linestyle='--', linewidth=1.5,
           label='Vacuum avg limit')

# Experimental data points
expts = [
    ('pp', 0.267, 0.57, 0.09),
    ('7Be', 0.862, 0.51, 0.07),
    ('pep', 1.44, 0.43, 0.11),
    ('8B', 8.0, 0.31, 0.04),
]
for name, E, P, err in expts:
    ax1.errorbar(E, P, yerr=err, fmt='o', color='#4ade80', markersize=8, capsize=4)
    ax1.annotate(name, xy=(E, P), xytext=(E*1.3, P+0.05), fontsize=8, color='#4ade80')

ax1.set_xlabel('Neutrino energy [MeV]')
ax1.set_ylabel('P(nu_e -> nu_e)')
ax1.set_title('Survival Probability vs Energy')
ax1.set_ylim(0, 0.8)
ax1.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: MSW level crossing diagram ---
ax2 = axes[0, 1]
A_range = np.linspace(0, 3, 500)

# Effective masses (eigenvalues of the 2x2 Hamiltonian)
# lambda_+/- = (A +/- sqrt((cos2theta - A)^2 + sin^2(2theta))) / 2
discriminant = np.sqrt((np.cos(2*theta12) - A_range)**2 + np.sin(2*theta12)**2)
lambda_plus = 0.5 * (A_range + discriminant)
lambda_minus = 0.5 * (A_range - discriminant)

# Also show the diagonal (unmixed) states
lambda_1_unmixed = np.zeros_like(A_range)
lambda_2_unmixed = np.cos(2*theta12) * np.ones_like(A_range)
lambda_e = A_range - np.cos(2*theta12)

ax2.plot(A_range, lambda_plus, color='#f43f5e', linewidth=2.5, label='nu_2^m (heavy)')
ax2.plot(A_range, lambda_minus, color='#60a5fa', linewidth=2.5, label='nu_1^m (light)')
ax2.plot(A_range, lambda_e - np.cos(2*theta12), color='#fbbf24', linewidth=1, linestyle=':', alpha=0.5, label='nu_e (unmixed)')
ax2.plot(A_range, lambda_1_unmixed, color='#4ade80', linewidth=1, linestyle=':', alpha=0.5, label='nu_mu (unmixed)')

# Mark resonance
A_res = np.cos(2*theta12)
ax2.axvline(x=A_res, color='white', linestyle='--', alpha=0.3)
ax2.annotate('Resonance', xy=(A_res, 0.5), fontsize=9, color='white')

ax2.set_xlabel('Matter parameter A = 2sqrt(2) G_F n_e E / dm^2')
ax2.set_ylabel('Effective eigenvalue')
ax2.set_title('MSW Level-Crossing Diagram')
ax2.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: Electron density and resonance location ---
ax3 = axes[1, 0]
r_frac = np.linspace(0.01, 0.95, 200)

# Approximate electron density profile (exponential model)
n_e_profile = 6e25 * np.exp(-10.54 * r_frac)  # cm^-3

ax3.semilogy(r_frac, n_e_profile, color='#f43f5e', linewidth=2.5)

# Mark resonance locations for different energies
for E_val, col, lbl in [(1.0, '#fbbf24', '1 MeV'), (5.0, '#fb923c', '5 MeV'), (10.0, '#4ade80', '10 MeV')]:
    n_res = dm2_21 * np.cos(2*theta12) / (7.63e-14 * E_val)
    r_res_idx = np.argmin(np.abs(n_e_profile - n_res))
    if r_res_idx > 0 and r_res_idx < len(r_frac) - 1:
        ax3.axhline(y=n_res, color=col, linestyle='--', alpha=0.5)
        ax3.plot(r_frac[r_res_idx], n_e_profile[r_res_idx], 'o', color=col, markersize=10)
        ax3.annotate(lbl + ' res.', xy=(r_frac[r_res_idx], n_e_profile[r_res_idx]),
                    xytext=(r_frac[r_res_idx]+0.1, n_e_profile[r_res_idx]*2),
                    fontsize=8, color=col,
                    arrowprops=dict(arrowstyle='->', color=col))

ax3.set_xlabel('r / R_sun')
ax3.set_ylabel('Electron density [cm^-3]')
ax3.set_title('Resonance Locations in the Sun')

# --- Panel 4: Experimental results comparison ---
ax4 = axes[1, 1]

experiments = ['Homestake', 'SAGE', 'GALLEX', 'Super-K', 'SNO CC', 'SNO NC', 'Borexino']
measured = [2.56, 67.2, 69.3, 2.35, 1.76, 5.09, 3.0]
expected_ssm = [7.6, 128, 128, 5.05, 5.05, 5.05, 5.05]
# Normalize each to SSM prediction
ratios = [m/e for m, e in zip(measured, expected_ssm)]

colors_exp = ['#f43f5e', '#fb923c', '#fbbf24', '#4ade80', '#60a5fa', '#a78bfa', '#e879f9']
bars = ax4.barh(experiments, ratios, color=colors_exp, alpha=0.8, height=0.6)
ax4.axvline(x=1.0, color='white', linestyle='--', linewidth=1.5, label='SSM prediction')
ax4.set_xlabel('Measured / SSM prediction')
ax4.set_title('Solar Neutrino Experiments')
ax4.set_xlim(0, 1.2)

for i, (r, exp) in enumerate(zip(ratios, experiments)):
    ax4.text(r + 0.02, i, str(round(r, 2)), va='center', fontsize=8, color='white')

ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Plot saved to output.png")

print("\n" + "=" * 60)
print("SOLAR NEUTRINO PHYSICS")
print("=" * 60)
print("\n--- Oscillation Parameters ---")
print("  theta_12 = " + str(round(np.degrees(theta12), 2)) + " degrees")
print("  dm^2_21 = " + str(dm2_21) + " eV^2")
print("  sin^2(theta_12) = " + str(round(np.sin(theta12)**2, 4)))
print("\n--- Survival Probabilities ---")
P_vac_val = 1 - 0.5 * np.sin(2*theta12)**2
print("  Low-E vacuum avg: P_ee = " + str(round(P_vac_val, 4)))
print("  High-E MSW adiab: P_ee = sin^2(theta) = " + str(round(np.sin(theta12)**2, 4)))
print("\n--- Neutrino Fluxes (SSM, Vinyoles 2017) ---")
print("  pp:  5.98 x 10^10 cm^-2 s^-1")
print("  7Be: 4.93 x 10^9  cm^-2 s^-1")
print("  8B:  5.46 x 10^6  cm^-2 s^-1")
print("  CNO: 4.88 x 10^8  cm^-2 s^-1")
print("\n--- Nobel Prizes ---")
print("  2002: Ray Davis Jr. (Homestake), Masatoshi Koshiba (Kamiokande)")
print("  2015: Arthur McDonald (SNO), Takaaki Kajita (Super-K)")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Full Derivation: The MSW Effect

We derive the matter potential, the resonance condition, and the Landau-Zener crossing probability that governs the flavor conversion of solar neutrinos propagating through the electron-dense solar interior.

Step 1: Origin of the Matter Potential

In matter, electron neutrinos $\nu_e$ interact with electrons via both charged-current (W exchange) and neutral-current (Z exchange) interactions. Muon and tau neutrinos interact only via neutral current. The neutral-current contribution is the same for all flavors and cancels in the oscillation Hamiltonian. The remaining charged-current forward scattering gives an effective potential:

$$V_{\text{CC}} = \sqrt{2} \, G_F \, n_e$$

where $G_F = 1.166 \times 10^{-5}$ GeV$^{-2}$ is the Fermi constant and$n_e$ is the electron number density. At the solar center,$n_e \approx 6 \times 10^{25}$ cm$^{-3}$, giving:

$$\boxed{V = \sqrt{2} \, G_F \, n_e \approx 7.6 \times 10^{-12} \text{ eV at solar center}}$$

Step 2: Hamiltonian in Matter

In the flavor basis $(\nu_e, \nu_\mu)$, the effective Hamiltonian for propagation is:

$$H = \frac{\Delta m^2}{4E} \begin{pmatrix} -\cos 2\theta + A & \sin 2\theta \\ \sin 2\theta & \cos 2\theta - A \end{pmatrix}$$

where the dimensionless matter parameter is:

$$A = \frac{2\sqrt{2} \, G_F \, n_e \, E}{\Delta m^2} = \frac{2EV}{\Delta m^2}$$

Step 3: Resonance Condition

The effective mixing angle in matter $\theta_m$ is determined by diagonalizing $H$:

$$\tan 2\theta_m = \frac{\sin 2\theta}{\cos 2\theta - A}$$

Resonance (maximal mixing, $\theta_m = \pi/4$) occurs when the denominator vanishes:

$$\boxed{A_{\text{res}} = \cos 2\theta \implies n_e^{\text{res}} = \frac{\Delta m^2 \cos 2\theta}{2\sqrt{2} \, G_F \, E}}$$

The MSW resonance condition: mixing becomes maximal at a specific density

Step 4: Adiabaticity and the Landau-Zener Formula

As the neutrino propagates outward through decreasing density, it passes through the resonance. The key question is whether this passage is adiabatic (the neutrino stays in the same mass eigenstate). The adiabaticity parameter is:

$$\gamma = \frac{\Delta m^2 \sin^2 2\theta}{2E \cos 2\theta} \left|\frac{1}{n_e} \frac{dn_e}{dr}\right|_{\text{res}}^{-1}$$

The Landau-Zener crossing probability (probability of a non-adiabatic jump) is:

$$\boxed{P_c = \exp\left(-\frac{\pi}{2}\gamma\right) = \exp\left(-\frac{\pi \Delta m^2 \sin^2 2\theta}{4E \cos 2\theta \left|\frac{d\ln n_e}{dr}\right|_{\text{res}}}\right)}$$

For solar ${}^8\text{B}$ neutrinos with $E \sim 10$ MeV and the measured oscillation parameters, $\gamma \gg 1$ and $P_c \approx 0$: the transition is fully adiabatic.

Step 5: Survival Probability Formula

The general survival probability for $\nu_e$ produced at matter mixing angle$\theta_m^0$ in the core is:

$$P(\nu_e \to \nu_e) = \frac{1}{2} + \left(\frac{1}{2} - P_c\right) \cos 2\theta_m^0 \cos 2\theta$$

In the two limiting cases:

High energy (adiabatic MSW): $A \gg 1, \; \theta_m^0 \approx \pi/2, \; P_c = 0$

$$P_{ee} \approx \frac{1}{2}(1 + \cos 2\theta_m^0 \cos 2\theta) \to \sin^2\theta \approx 0.31$$

Low energy (vacuum oscillations): $A \ll 1, \; \theta_m^0 \approx \theta$

$$P_{ee} \approx 1 - \frac{1}{2}\sin^2 2\theta \approx 0.57$$

MSW Resonance Level-Crossing Diagram

The effective mass eigenvalues as a function of electron density. At resonance, the gap between the levels is minimized. An adiabatic neutrino follows the upper curve from high density (core) to vacuum.

A $\nu_e$ produced at high density in the core starts on the lower branch ($\nu_1^m \approx \nu_e$ at high density). As it propagates outward, an adiabatic transition follows the curve to emerge as $\nu_2$ in vacuum, with survival probability $P_{ee} = |\langle\nu_e|\nu_2\rangle|^2 = \sin^2\theta$.

Extended Simulation: MSW Survival Probability

Detailed computation of the neutrino survival probability comparing vacuum oscillations with the full MSW matter effect across the solar neutrino energy spectrum.

MSW Survival Probability, Matter Mixing, Level Crossing, and Landau-Zener Adiabaticity

Python

script.py194 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('MSW Effect: Survival Probability and Matter Oscillations', fontsize=16,
             fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

# Oscillation parameters
theta12 = np.radians(33.44)
dm2 = 7.53e-5  # eV^2
G_F = 1.166e-5  # GeV^-2
n_e_center = 6e25  # cm^-3

# --- Panel 1: P_ee vacuum vs MSW ---
ax1 = axes[0, 0]
E_arr = np.logspace(-1, 1.3, 1000)  # MeV

# Matter parameter at solar center
A = 7.63e-14 * n_e_center * E_arr / dm2

# Matter mixing angle at production point
cos2theta_m = (np.cos(2*theta12) - A) / np.sqrt((np.cos(2*theta12) - A)**2 + np.sin(2*theta12)**2)

# MSW adiabatic survival probability
P_msw = 0.5 + 0.5 * cos2theta_m * np.cos(2*theta12)

# Vacuum average survival probability (distance-averaged)
P_vac = (1 - 0.5 * np.sin(2*theta12)**2) * np.ones_like(E_arr)

# Full formula with crossing probability (approximate)
gamma = dm2 * np.sin(2*theta12)**2 / (2 * E_arr * 1e-3 * np.cos(2*theta12) * 10.54 / 6.957e10)
P_c = np.exp(-np.pi * gamma / 2)
P_c = np.clip(P_c, 0, 1)
P_full = 0.5 + (0.5 - P_c) * cos2theta_m * np.cos(2*theta12)

ax1.semilogx(E_arr, P_vac, color='#60a5fa', linewidth=2, linestyle='--', label='Vacuum (averaged)')
ax1.semilogx(E_arr, P_msw, color='#f43f5e', linewidth=2.5, label='MSW adiabatic')
ax1.semilogx(E_arr, P_full, color='#fbbf24', linewidth=1.5, linestyle=':', label='Full (with P_c)')

# Experimental data
data = [
    ('pp', 0.267, 0.55, 0.07),
    ('7Be', 0.862, 0.51, 0.05),
    ('pep', 1.44, 0.43, 0.10),
    ('8B (SK)', 8.0, 0.31, 0.03),
    ('8B (SNO)', 10.0, 0.34, 0.04),
]
for name, E_val, P_val, err in data:
    ax1.errorbar(E_val, P_val, yerr=err, fmt='o', color='#4ade80', markersize=7, capsize=4)
    ax1.annotate(name, xy=(E_val, P_val), xytext=(E_val*1.4, P_val + 0.04),
                fontsize=7, color='#4ade80')

ax1.axhline(y=np.sin(theta12)**2, color='#fb923c', linestyle=':', alpha=0.5,
           label='sin^2(theta_12) = ' + str(round(np.sin(theta12)**2, 3)))
ax1.set_xlabel('Neutrino energy [MeV]')
ax1.set_ylabel('P(nu_e -> nu_e)')
ax1.set_title('Survival Probability: Vacuum vs MSW')
ax1.set_ylim(0, 0.8)
ax1.legend(fontsize=6, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Effective mixing angle vs density ---
ax2 = axes[0, 1]
A_range = np.linspace(0, 5, 500)

sin2_2theta_m = np.sin(2*theta12)**2 / ((np.cos(2*theta12) - A_range)**2 + np.sin(2*theta12)**2)
theta_m = 0.5 * np.arctan2(np.sin(2*theta12), np.cos(2*theta12) - A_range)

ax2.plot(A_range, np.degrees(theta_m), color='#f43f5e', linewidth=2.5, label='theta_m')
ax2.axhline(y=45, color='white', linestyle=':', alpha=0.3)
ax2.axhline(y=np.degrees(theta12), color='#fbbf24', linestyle='--', alpha=0.5,
           label='Vacuum theta_12 = ' + str(round(np.degrees(theta12), 1)) + ' deg')
ax2.axvline(x=np.cos(2*theta12), color='#4ade80', linestyle=':', alpha=0.5, label='Resonance (A=cos2theta)')

ax2.set_xlabel('Matter parameter A')
ax2.set_ylabel('Effective mixing angle [degrees]')
ax2.set_title('Matter Mixing Angle vs Density')
ax2.set_ylim(0, 90)
ax2.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: MSW level crossing with energy annotations ---
ax3 = axes[1, 0]

discriminant = np.sqrt((np.cos(2*theta12) - A_range)**2 + np.sin(2*theta12)**2)
lam_plus = 0.5 * (A_range + discriminant)
lam_minus = 0.5 * (A_range - discriminant)

ax3.plot(A_range, lam_plus, color='#f43f5e', linewidth=2.5, label='nu_2^m (heavy)')
ax3.plot(A_range, lam_minus, color='#60a5fa', linewidth=2.5, label='nu_1^m (light)')

# Unmixed flavor diagonal
ax3.plot(A_range, A_range - np.cos(2*theta12), color='#fbbf24', linewidth=1,
         linestyle=':', alpha=0.4, label='nu_e (unmixed)')
ax3.plot(A_range, np.zeros_like(A_range), color='#4ade80', linewidth=1,
         linestyle=':', alpha=0.4, label='nu_mu (unmixed)')

# Resonance vertical
A_res = np.cos(2*theta12)
ax3.axvline(x=A_res, color='white', linestyle='--', alpha=0.3)
ax3.annotate('Resonance', xy=(A_res, 0.4), fontsize=9, color='white')

# Mark energy-corresponding A values
for E_val, col, lbl in [(1, '#fbbf24', '1 MeV'), (5, '#fb923c', '5 MeV'), (10, '#f43f5e', '10 MeV')]:
    A_val = 7.63e-14 * n_e_center * E_val / dm2
    if A_val < 5:
        ax3.axvline(x=A_val, color=col, linestyle='--', alpha=0.3)
        ax3.text(A_val, 2.3, lbl, fontsize=7, color=col, rotation=90, va='bottom')

ax3.set_xlabel('Matter parameter A')
ax3.set_ylabel('Effective eigenvalue (arb. units)')
ax3.set_title('MSW Level-Crossing Diagram')
ax3.legend(fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')
ax3.set_ylim(-1.5, 2.5)

# --- Panel 4: Landau-Zener adiabaticity ---
ax4 = axes[1, 1]
E_lz = np.logspace(-1, 2, 500)  # MeV

# Scale height for solar density profile
H_n = 6.957e10 / 10.54  # cm (1/e scale of density)

# Adiabaticity parameter
gamma_lz = dm2 * np.sin(2*theta12)**2 * H_n / (2 * E_lz * 1e-3 * np.cos(2*theta12))
# Convert units properly
gamma_lz = gamma_lz * 5.068e9  # GeV^-1 to cm conversion factor

P_jump = np.exp(-np.pi * gamma_lz / 2)
P_jump = np.clip(P_jump, 0, 1)

ax4.loglog(E_lz, gamma_lz, color='#f43f5e', linewidth=2.5, label='Adiabaticity gamma')
ax4.axhline(y=1, color='white', linestyle='--', alpha=0.3)
ax4.text(0.15, 1.5, 'gamma = 1 (adiabatic boundary)', fontsize=8, color='white', alpha=0.5)

# Mark solar neutrino energies
for E_val, lbl, col in [(0.267, 'pp', '#4ade80'), (0.862, '7Be', '#fbbf24'),
                         (1.44, 'pep', '#fb923c'), (10, '8B', '#f43f5e')]:
    idx = np.argmin(np.abs(E_lz - E_val))
    ax4.plot(E_val, gamma_lz[idx], 'o', color=col, markersize=8)
    ax4.annotate(lbl, xy=(E_val, gamma_lz[idx]),
                xytext=(E_val*1.5, gamma_lz[idx]*2),
                fontsize=8, color=col,
                arrowprops=dict(arrowstyle='->', color=col))

ax4.fill_between(E_lz, gamma_lz, 1,
                 where=gamma_lz > 1, alpha=0.05, color='#4ade80')
ax4.fill_between(E_lz, gamma_lz, 1,
                 where=gamma_lz < 1, alpha=0.05, color='#f43f5e')
ax4.text(20, 0.1, 'Non-adiabatic', fontsize=9, color='#f43f5e', alpha=0.6)
ax4.text(0.15, 100, 'Adiabatic', fontsize=9, color='#4ade80', alpha=0.6)

ax4.set_xlabel('Neutrino energy [MeV]')
ax4.set_ylabel('Adiabaticity parameter gamma')
ax4.set_title('Landau-Zener Adiabaticity')
ax4.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Plot saved to output.png")

print("\n" + "=" * 60)
print("MSW EFFECT: DETAILED ANALYSIS")
print("=" * 60)
print("\n--- Oscillation Parameters ---")
print("  theta_12 = " + str(round(np.degrees(theta12), 2)) + " degrees")
print("  dm^2_21 = " + str(dm2) + " eV^2")
print("  sin^2(theta) = " + str(round(np.sin(theta12)**2, 4)))
print("  cos(2*theta) = " + str(round(np.cos(2*theta12), 4)))
print("\n--- Survival Probability Limits ---")
print("  Low-E vacuum: P_ee = 1 - 0.5*sin^2(2theta) = " + str(round(1 - 0.5*np.sin(2*theta12)**2, 4)))
print("  High-E MSW:   P_ee = sin^2(theta) = " + str(round(np.sin(theta12)**2, 4)))
print("  Transition region: 2-5 MeV")
print("\n--- Resonance Conditions ---")
A_res_val = np.cos(2*theta12)
n_res_10MeV = dm2 * np.cos(2*theta12) / (7.63e-14 * 10)
print("  A_res = cos(2*theta) = " + str(round(A_res_val, 4)))
print("  n_e(res) for 10 MeV = " + str(round(n_res_10MeV / 1e25, 2)) + " x 10^25 cm^-3")
print("\n--- Key Experiments ---")
print("  SNO NC flux matches SSM: neutrinos oscillate, not disappear")
print("  Borexino confirms energy-dependent transition")
print("  KamLAND confirms same dm^2 with reactor antineutrinos")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Helioseismology Next: Photosphere →