Part I: Solar Interior | Chapter 3

Helioseismology

Probing the solar interior with acoustic and gravity modes

3.1 Pressure Modes (p-modes)

The Sun acts as a resonant cavity for acoustic waves. Turbulent convection near the surface stochastically excites standing sound waves that penetrate to different depths depending on their angular degree $\ell$. These pressure modes, or p-modes, have periods of roughly 5 minutes and have been observed for millions of individual oscillation frequencies.

Derivation 1: Adiabatic Wave Equation

Step 1. Start from the linearized equations of hydrodynamics. Perturb all quantities about the equilibrium state: $\rho = \rho_0 + \rho', \; P = P_0 + P', \; \mathbf{v} = \mathbf{v}'$. The linearized momentum equation is:

$$\rho_0 \frac{\partial \mathbf{v}'}{\partial t} = -\nabla P' + \rho' \mathbf{g}$$

Step 2. The linearized continuity equation:

$$\frac{\partial \rho'}{\partial t} + \nabla \cdot (\rho_0 \mathbf{v}') = 0$$

Step 3. For adiabatic perturbations, $P' = c_s^2 \rho'$ where$c_s^2 = \Gamma_1 P_0 / \rho_0$ is the adiabatic sound speed squared. Introducing the displacement vector $\boldsymbol{\xi}$ via $\mathbf{v}' = \partial\boldsymbol{\xi}/\partial t$, we obtain the oscillation equation:

$$\boxed{\rho_0 \frac{\partial^2 \boldsymbol{\xi}}{\partial t^2} = \nabla(c_s^2 \rho_0 \nabla \cdot \boldsymbol{\xi} + \boldsymbol{\xi} \cdot \nabla P_0) - \mathbf{g} \nabla \cdot (\rho_0 \boldsymbol{\xi})}$$

Step 4. Decompose $\boldsymbol{\xi}$ in spherical harmonics:$\xi_r(r, \theta, \phi, t) = \xi_r(r) Y_\ell^m(\theta, \phi) e^{-i\omega t}$. Each mode is characterized by three quantum numbers: radial order $n$, angular degree $\ell$, and azimuthal order $m$.

For p-modes, the restoring force is pressure. The turning point radius where a mode of degree $\ell$ is reflected inward is determined by$c_s(r_t) / r_t = \omega / L$ where $L = \sqrt{\ell(\ell+1)}$. Low-$\ell$ modes penetrate deep into the core; high-$\ell$ modes are confined near the surface.

3.2 Gravity Modes (g-modes)

Derivation 2: Brunt-Vaisala Frequency

Gravity modes are driven by buoyancy as the restoring force. They exist in stably stratified regions (the radiative zone) and are evanescent in convection zones.

Step 1. Consider a fluid element displaced vertically by $\delta r$in a stratified atmosphere. The buoyancy force per unit volume is:

$$f = -g(\rho_{\text{element}} - \rho_{\text{surroundings}})$$

Step 2. If the element moves adiabatically and the surroundings have a different temperature gradient, expanding to first order:

$$f = \rho_0 N^2 \delta r$$

where the Brunt-Vaisala (buoyancy) frequency is:

$$\boxed{N^2 = g\left(\frac{1}{\Gamma_1}\frac{d\ln P}{dr} - \frac{d\ln\rho}{dr}\right) = \frac{g^2\rho}{P}\left(\nabla_{\text{ad}} - \nabla + \nabla_\mu\right)}$$

Step 3. In the radiative zone, $N^2 > 0$ (stable stratification), and g-modes oscillate with frequencies $\omega < N$. In the convection zone,$N^2 < 0$, and g-modes become evanescent.

Solar g-modes would provide direct information about the core rotation and composition gradients. Despite decades of effort, their detection remains controversial, with claimed detections by the GOLF instrument on SOHO showing period spacings consistent with $\Delta P \approx 24$ minutes for dipole ($\ell=1$) modes.

3.3 Asymptotic Theory (Tassoul)

Derivation 3: Tassoul's Asymptotic Formula

For high-order, low-degree p-modes, the JWKB (Jeffreys-Wentzel-Kramers-Brillouin) approximation yields a remarkably simple expression for the frequencies.

Step 1. In the asymptotic limit ($n \gg \ell$), the radial eigenfunction oscillates many times between the turning points. The JWKB quantization condition gives:

$$\int_{r_t}^{R} \left(\frac{\omega^2}{c_s^2} - \frac{L^2}{r^2}\right)^{1/2} dr = (n + \alpha)\pi$$

where $\alpha$ is a phase correction related to the surface reflection.

Step 2. For low-degree modes where $r_t \to 0$, this simplifies using the Duvall law. Expanding to leading order, Tassoul (1980) showed:

$$\boxed{\nu_{n,\ell} \approx \Delta\nu\left(n + \frac{\ell}{2} + \epsilon\right) - D_0 \ell(\ell+1)}$$

where the large frequency separation is:

$$\Delta\nu = \left(2\int_0^R \frac{dr}{c_s}\right)^{-1} \approx 135 \; \mu\text{Hz for the Sun}$$

Step 3. The small frequency separation is:

$$\delta\nu_{02} = \nu_{n,0} - \nu_{n-1,2} \approx 6D_0 \approx 9 \; \mu\text{Hz}$$

The large separation $\Delta\nu$ depends on the mean density (sound travel time across the full diameter), while the small separation $\delta\nu$ is sensitive to conditions in the core, making it an excellent probe of the central hydrogen abundance and hence stellar age.

3.4 Rotational Splitting

Derivation 4: Rotational Frequency Splitting

Solar rotation breaks the degeneracy of modes with different azimuthal orders $m$. This splitting provides a direct measurement of the internal rotation profile.

Step 1. In a rotating frame, each mode frequency is shifted by an amount proportional to $m$. To first order in the rotation rate $\Omega$:

$$\nu_{n\ell m} = \nu_{n\ell 0} + m \int_0^R \int_0^\pi K_{n\ell}(r, \theta) \, \Omega(r, \theta) \, r \, dr \, d\theta$$

Step 2. The rotational kernel $K_{n\ell}(r,\theta)$ is determined by the mode eigenfunctions. For a purely radially-dependent rotation $\Omega(r)$:

$$\boxed{\delta\nu_{n\ell m} = m \int_0^R K_{n\ell}(r) \, \Omega(r) \, dr}$$

Step 3. The splitting coefficients are commonly expressed as a polynomial expansion:$\delta\nu = m\sum_{j} a_j(n,\ell) P_j(\ell, m)$ where $a_1$gives the mean rotation and $a_3, a_5$ probe the latitude dependence.

Helioseismic inversions have revealed that: (1) the convection zone rotates differentially (equator ~25 days, poles ~35 days), (2) the radiative zone rotates nearly uniformly at ~27 days, and (3) a thin shear layer called the tachocline exists at the base of the convection zone ($r \approx 0.693 R_\odot$), believed to be crucial for the solar dynamo.

3.5 Sound Speed Inversions and the Duvall Law

Derivation 5: The Duvall Law

The Duvall law provides a one-dimensional function that all p-mode frequencies must satisfy, and it forms the basis for inverting the observed frequencies to obtain the sound speed profile.

Step 1. For a spherically symmetric model, define the acoustic radius$\tau = \int_0^r dr'/c_s(r')$. The ray-theory turning point condition is:

$$\frac{c_s(r_t)}{r_t} = \frac{\omega}{L}$$

Step 2. Duvall (1982) showed that all p-mode frequencies can be collapsed onto a single function:

$$\boxed{\frac{n + \alpha}{\nu} = F\left(\frac{\nu}{L}\right) = \int_{r_t}^{R} \left(1 - \frac{L^2 c_s^2}{\omega^2 r^2}\right)^{1/2} \frac{dr}{c_s}}$$

Step 3. The sound speed profile is obtained by Abel inversion:

$$r = R \exp\left(-\frac{2}{\pi}\int_w^{w_{\max}} \frac{dF/dw}{\sqrt{w^2 - w_t^2}} \, dw\right)$$

where $w = \nu/L$. This integral equation can be solved for $c_s(r)$given the observed $F(w)$ function.

Modern inversions using regularized least-squares (RLS) or optimally localized averages (OLA) methods achieve relative precision of $\delta c_s/c_s \sim 10^{-4}$ throughout most of the Sun. The inversions confirm the SSM to better than 0.2% except near the base of the convection zone where the solar abundance problem manifests.

Historical Context

The 5-minute oscillations were discovered by Leighton, Noyes, and Simon (1962). Ulrich (1970) and Leibacher and Stein (1971) identified them as trapped acoustic modes. Deubner (1975) confirmed the diagnostic $\ell$-$\nu$ diagram. The field was revolutionized by the BiSON, GONG, and MDI/SOHO networks providing continuous, high-resolution frequency measurements.

Numerical Simulation

This simulation computes the p-mode power spectrum, ray paths for different angular degrees, the Duvall law collapse, and the internal rotation profile from splitting measurements.

Helioseismology: Mode Frequencies, Ray Paths, Sound Speed, and Rotation

Python

script.py162 lines

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

# ============================================================
# Helioseismology Simulations
# ============================================================

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Helioseismology: Probing the Solar Interior', 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')

# --- Panel 1: l-nu diagram (echelle-like) ---
ax1 = axes[0, 0]
Delta_nu = 135.0  # microHz
D_0 = 1.5  # microHz
epsilon = 1.5

ell_vals = [0, 1, 2, 3]
colors_l = ['#f43f5e', '#fb923c', '#fbbf24', '#4ade80']
for ell, col in zip(ell_vals, colors_l):
    n_vals = np.arange(10, 30)
    nu = Delta_nu * (n_vals + ell / 2.0 + epsilon) - D_0 * ell * (ell + 1)
    # Add small random scatter to simulate real data
    np.random.seed(42 + ell)
    nu += np.random.normal(0, 0.3, len(nu))
    ax1.scatter(ell * np.ones_like(nu) + np.random.normal(0, 0.05, len(nu)),
               nu, c=col, s=15, alpha=0.8, label='l = ' + str(ell))

ax1.set_xlabel('Angular degree l')
ax1.set_ylabel('Frequency [microHz]')
ax1.set_title('p-mode Frequencies (Tassoul Pattern)')
ax1.set_xlim(-0.5, 3.5)
ax1.set_ylim(1500, 4200)
ax1.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Ray paths ---
ax2 = axes[0, 1]

# Draw solar disk
theta_circle = np.linspace(0, 2*np.pi, 200)
ax2.plot(np.cos(theta_circle), np.sin(theta_circle), 'w-', linewidth=1.5)

# Sound speed profile (simplified)
def c_s(r):
    # Approximate sound speed in km/s, normalized
    return 0.1 + 0.5 * (1 - r**2)

# Ray tracing for different l values
for ell, col, lbl in [(1, '#f43f5e', 'l=1'), (5, '#fb923c', 'l=5'),
                       (20, '#fbbf24', 'l=20'), (75, '#4ade80', 'l=75')]:
    # Turning point
    L = np.sqrt(ell * (ell + 1))
    omega = 2 * np.pi * 3000e-6  # ~3 mHz
    # r_t where c_s/r = omega/L
    r_t_approx = max(0.05, 1.0 - ell / 100.0)

# Draw approximate ray path as an arc
    n_pts = 100
    theta_ray = np.linspace(-0.8, 0.8, n_pts)
    r_ray = r_t_approx + (1.0 - r_t_approx) * np.cos(theta_ray * np.pi / 1.6)**2

x_ray = r_ray * np.cos(theta_ray + np.pi/4)
    y_ray = r_ray * np.sin(theta_ray + np.pi/4)

ax2.plot(x_ray, y_ray, color=col, linewidth=1.5, alpha=0.8, label=lbl)
    # Draw turning point circle
    theta_tp = np.linspace(0, 2*np.pi, 100)
    ax2.plot(r_t_approx * np.cos(theta_tp), r_t_approx * np.sin(theta_tp),
            color=col, linewidth=0.5, linestyle=':', alpha=0.4)

ax2.set_xlim(-1.2, 1.2)
ax2.set_ylim(-1.2, 1.2)
ax2.set_aspect('equal')
ax2.set_title('Ray Paths (different l)')
ax2.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')
ax2.set_xlabel('x / R_sun')
ax2.set_ylabel('y / R_sun')

# --- Panel 3: Sound speed profile ---
ax3 = axes[1, 0]
r_arr = np.linspace(0.05, 0.99, 200)

# Model sound speed (approximate SSM values in km/s)
cs_model = 510 * (1 - 0.8 * r_arr**1.3) * np.exp(-0.1 * r_arr)
cs_model = np.clip(cs_model, 10, 520)

# "Observed" from inversions (slightly different)
np.random.seed(123)
cs_obs = cs_model * (1 + 0.001 * np.sin(5 * np.pi * r_arr) + 0.0005 * np.random.randn(len(r_arr)))

ax3.plot(r_arr, cs_model, color='#f43f5e', linewidth=2.5, label='SSM (model)')
ax3.plot(r_arr, cs_obs, color='#60a5fa', linewidth=2, linestyle='--', label='Helioseismic inversion')
ax3.axvline(x=0.25, color='#fbbf24', linestyle=':', alpha=0.5, label='Core boundary')
ax3.axvline(x=0.713, color='#4ade80', linestyle=':', alpha=0.5, label='CZ base')
ax3.set_xlabel('r / R_sun')
ax3.set_ylabel('Sound speed [km/s]')
ax3.set_title('Sound Speed Profile')
ax3.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 4: Internal rotation ---
ax4 = axes[1, 1]
r_rot = np.linspace(0.1, 0.95, 200)

# Latitude-dependent rotation (nHz)
lats = [0, 30, 60]  # degrees
colors_lat = ['#f43f5e', '#fb923c', '#fbbf24']

for lat, col in zip(lats, colors_lat):
    lat_rad = np.radians(lat)
    # Surface differential rotation
    Omega_surf = 460 - 65 * np.sin(lat_rad)**2 - 55 * np.sin(lat_rad)**4  # nHz
    # Radiative interior: solid body
    Omega_rad = 430  # nHz

# Tachocline transition
    r_tach = 0.693
    w_tach = 0.02
    Omega = Omega_rad + (Omega_surf - Omega_rad) * 0.5 * (1 + np.tanh((r_rot - r_tach) / w_tach))

ax4.plot(r_rot, Omega, color=col, linewidth=2.5, label=str(lat) + ' deg latitude')

ax4.axvline(x=0.693, color='white', linestyle=':', alpha=0.3, label='Tachocline')
ax4.set_xlabel('r / R_sun')
ax4.set_ylabel('Rotation rate [nHz]')
ax4.set_title('Internal Rotation Profile')
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("HELIOSEISMOLOGY RESULTS")
print("=" * 60)
print("\n--- Asymptotic Parameters ---")
print("  Large separation Delta_nu = " + str(Delta_nu) + " microHz")
print("  Small separation delta_nu_02 ~ " + str(round(6 * D_0, 1)) + " microHz")
print("  Phase epsilon = " + str(epsilon))
print("\n--- Sound Speed ---")
print("  Central c_s ~ 510 km/s")
print("  Surface c_s ~ 10 km/s")
print("  Sound crossing time ~ 2 x R_sun / <c_s> ~ " + str(round(2 * 6.957e8 / (200e3) / 3600, 1)) + " hours")
print("\n--- Internal Rotation ---")
print("  Equatorial surface: ~460 nHz (25.4 days)")
print("  Polar surface: ~340 nHz (34.0 days)")
print("  Radiative interior: ~430 nHz (26.9 days)")
print("  Tachocline: r ~ 0.693 R_sun, width ~ 0.04 R_sun")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Full Derivation: Tassoul's Asymptotic Formula

The remarkably regular spacing of solar p-mode frequencies can be derived from the JWKB approximation applied to the radial oscillation equation. We present the full derivation leading to the large and small frequency separations.

Step 1: The Radial Wave Equation

In the Cowling approximation (neglecting the perturbation to the gravitational potential), the radial displacement satisfies:

$$\frac{d^2\Psi}{dr^2} + \left[\frac{\omega^2}{c_s^2} - \frac{L^2}{r^2} - \omega_c^2\right]\Psi = 0$$

where $\Psi$ is related to the radial displacement,$L = \sqrt{\ell(\ell+1)}$ is the angular degree factor, and $\omega_c$ is the acoustic cutoff frequency.

Step 2: JWKB Quantization Condition

For high-order modes ($n \gg \ell$), the JWKB approximation gives standing wave solutions between the inner turning point $r_t$ (where $c_s/r = \omega/L$) and the surface. The quantization condition requires an integer number of half-wavelengths:

$$\int_{r_t}^{R} k_r \, dr = \int_{r_t}^{R} \left(\frac{\omega^2}{c_s^2} - \frac{L^2}{r^2}\right)^{1/2} dr = \left(n + \alpha\right)\pi$$

Step 3: Low-Degree Limit

For low-degree modes ($\ell = 0, 1, 2, 3$), the turning point approaches the center ($r_t \to 0$). In this limit, the integral becomes:

$$\int_0^{R} \frac{\omega}{c_s} dr - \frac{L}{\omega} \int_0^R \frac{c_s}{r^2} \, dr + \cdots = (n + \alpha)\pi$$

The leading term gives:

$$\omega \int_0^R \frac{dr}{c_s} \approx (n + \alpha)\pi$$

Step 4: Large Frequency Separation

Converting to cyclic frequency $\nu = \omega/(2\pi)$:

$$\nu \approx \frac{n + \alpha}{2\int_0^R dr/c_s}$$

The large frequency separation is the frequency difference between consecutive radial orders:

$$\boxed{\Delta\nu = \nu_{n+1,\ell} - \nu_{n,\ell} = \left(2\int_0^R \frac{dr}{c_s}\right)^{-1} \approx 135 \; \mu\text{Hz}}$$

$\Delta\nu$ is the inverse of the sound travel time across the full solar diameter

Step 5: Including the $\ell$-Dependence

Tassoul (1980) carried the expansion to higher order, including the second integral which depends on $\ell$. The result is:

$$\boxed{\nu_{n,\ell} \approx \Delta\nu\left(n + \frac{\ell}{2} + \epsilon\right) - D_0\,\ell(\ell+1)}$$

where $\epsilon \approx 1.5$ is a phase offset and the second-order term is:

$$D_0 = -\frac{\Delta\nu}{4\pi^2\nu_{n,\ell}} \int_0^R \frac{dc_s}{dr}\frac{dr}{r}$$

Step 6: Small Frequency Separation

The key consequence of the $\ell/2$ term is a near-degeneracy:$\nu_{n,\ell} \approx \nu_{n-1,\ell+2}$. The small frequency separation measures their difference:

$$\boxed{\delta\nu_{02} = \nu_{n,0} - \nu_{n-1,2} \approx 6D_0 \approx 9 \; \mu\text{Hz}}$$

$\delta\nu_{02}$ probes the sound speed gradient in the core, sensitive to the central H abundance (age)

Full Derivation: Rotational Frequency Splitting

Step 1: Rotation as a Perturbation

In a non-rotating star, modes with the same $n$ and $\ell$ but different azimuthal order $m$ are degenerate (same frequency). Rotation breaks this degeneracy. Treating rotation as a first-order perturbation:

$$\delta\omega_{n\ell m} = m \int_0^R \int_0^\pi K_{n\ell}(r,\theta) \, \Omega(r,\theta) \, r \, dr \, d\theta$$

Step 2: Rotational Kernel

The rotational kernel $K_{n\ell}(r,\theta)$ is constructed from the mode eigenfunctions$\xi_r(r)$ (radial) and $\xi_h(r)$ (horizontal):

$$K_{n\ell}(r) = \frac{\rho_0 r^2 \left[\xi_r^2 + L^2\xi_h^2 - 2\xi_r\xi_h - \xi_h^2\right]}{\int_0^R \rho_0 r^2 \left[\xi_r^2 + L^2\xi_h^2\right] dr}$$

Step 3: 1D Inversion for Radial Rotation

For a latitude-independent rotation profile $\Omega(r)$, the splitting simplifies to:

$$\boxed{\delta\nu_{n\ell m} = m \int_0^R K_{n\ell}(r) \, \Omega(r) \, dr}$$

Step 4: Latitude-Dependent Splitting

Real solar rotation depends on latitude. The splitting is expanded in odd-order a-coefficients:

$$\delta\nu_{n\ell m} = \sum_{j=1,3,5,...} a_j(n,\ell) \, \mathcal{P}_j^{(\ell)}(m)$$

where $a_1$ gives the mean rotation rate, $a_3$ measures the equator-to-pole differential rotation, and $a_5$ captures finer latitudinal structure. Helioseismic inversions of the full 2D rotation profile $\Omega(r,\theta)$ use thousands of measured splittings with regularization techniques.

Key Results from Rotational Inversions

• The convection zone exhibits strong differential rotation: equator rotates ~30% faster than poles
• The radiative interior rotates nearly as a solid body at ~430 nHz (~27 day period)
• The tachocline at $r \approx 0.693 R_\odot$ is a thin shear layer (width ~0.04 $R_\odot$)
• Near-surface shear layer (NSSL) shows rotation rate decreasing just below the photosphere

$\ell$-$\nu$ Diagram (Echelle Diagram)

In an echelle diagram, frequencies are plotted modulo the large separation $\Delta\nu$. Modes of the same $\ell$ form vertical ridges, offset by $\ell/2 \cdot \Delta\nu$from Tassoul's formula:

The l=0 and l=2 ridges are close together (separated by $\delta\nu_{02} \approx 9$ $\mu$Hz), as are l=1 and l=3. This near-degeneracy is a direct consequence of Tassoul's $\ell/2$ term.

Extended Simulation: Echelle Diagram and Frequency Separations

Compute the large separation $\Delta\nu$, small separation $\delta\nu_{02}$, and construct a synthetic echelle diagram from asymptotic p-mode theory.

Echelle Diagram, Frequency Separations, C-D Diagram, and Synthetic Power Spectrum

Python

script.py169 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('Helioseismology: Echelle Diagram and Frequency Analysis', fontsize=16,
             fontweight='bold', color='white')
fig.patch.set_facecolor('#0a0a0a')

# Asymptotic parameters (solar values)
Delta_nu = 134.9    # microHz
D_0 = 1.5           # microHz
epsilon = 1.55

# --- Panel 1: Echelle diagram ---
ax1 = axes[0, 0]
colors_ell = ['#f43f5e', '#fb923c', '#fbbf24', '#4ade80']
labels_ell = ['l=0', 'l=1', 'l=2', 'l=3']
sizes_ell = [40, 40, 25, 20]

for ell, col, lbl, sz in zip([0,1,2,3], colors_ell, labels_ell, sizes_ell):
    n_vals = np.arange(8, 28)
    nu = Delta_nu * (n_vals + ell / 2.0 + epsilon) - D_0 * ell * (ell + 1)
    # Add small scatter to simulate measurement uncertainty
    np.random.seed(100 + ell)
    nu += np.random.normal(0, 0.2, len(nu))
    nu_mod = nu % Delta_nu
    ax1.scatter(nu_mod, nu, c=col, s=sz, alpha=0.9, label=lbl, zorder=3)

ax1.set_xlabel('Frequency mod ' + str(Delta_nu) + ' [microHz]')
ax1.set_ylabel('Frequency [microHz]')
ax1.set_title('Echelle Diagram')
ax1.set_xlim(0, Delta_nu)
ax1.legend(fontsize=8, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 2: Large and small separations vs frequency ---
ax2 = axes[0, 1]
n_vals = np.arange(10, 26)
nu_0 = Delta_nu * (n_vals + 0.0 / 2.0 + epsilon)
nu_1 = Delta_nu * (n_vals + 1.0 / 2.0 + epsilon) - D_0 * 2
nu_2 = Delta_nu * (n_vals + 2.0 / 2.0 + epsilon) - D_0 * 6

# Large separation (varies slightly with frequency in real data)
np.random.seed(42)
Dnu_measured = np.diff(nu_0) + np.random.normal(0, 0.3, len(n_vals)-1)
ax2.plot(nu_0[:-1], Dnu_measured, 'o-', color='#f43f5e', markersize=5, linewidth=1.5,
         label='Large sep (l=0)')

# Small separation
delta_02 = nu_0[1:] - nu_2[:-1]
np.random.seed(43)
delta_02 += np.random.normal(0, 0.2, len(delta_02))
ax2_twin = ax2.twinx()
ax2_twin.plot(nu_0[1:], delta_02, 's-', color='#fbbf24', markersize=5, linewidth=1.5,
              label='Small sep delta_02')
ax2_twin.set_ylabel('delta_nu_02 [microHz]', color='#fbbf24')
ax2_twin.tick_params(axis='y', colors='#fbbf24')

ax2.set_xlabel('Frequency [microHz]')
ax2.set_ylabel('Delta_nu [microHz]', color='#f43f5e')
ax2.tick_params(axis='y', colors='#f43f5e')
ax2.set_title('Frequency Separations')
lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
ax2.legend(lines1 + lines2, labels1 + labels2,
          fontsize=7, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white')

# --- Panel 3: C-D diagram (asteroseismic HR diagram) ---
ax3 = axes[1, 0]
# Solar values
Dnu_sun = 134.9
dnu02_sun = 9.0

# Other stars for context (approximate)
stars = [
    ('Sun', 134.9, 9.0, '#f43f5e', 120),
    ('alpha Cen A', 105.5, 6.2, '#fb923c', 80),
    ('alpha Cen B', 161.0, 11.2, '#fbbf24', 80),
    ('KIC 8006161', 149.3, 10.5, '#4ade80', 60),
    ('16 Cyg A', 103.2, 5.8, '#60a5fa', 60),
    ('16 Cyg B', 117.0, 7.1, '#a78bfa', 60),
]

for name, dnu, d02, col, sz in stars:
    ax3.scatter(dnu, d02, c=col, s=sz, zorder=3, label=name)
    ax3.annotate(name, xy=(dnu, d02), xytext=(dnu+3, d02+0.3),
                fontsize=7, color=col)

# Age tracks (schematic)
dnu_track = np.linspace(90, 170, 50)
for age, alpha_val in [(2, 0.4), (4.6, 0.6), (8, 0.8)]:
    d02_track = 14 - alpha_val * (170 - dnu_track) / 10
    ax3.plot(dnu_track, d02_track, '--', color='white', alpha=0.2, linewidth=1)
    ax3.text(170, d02_track[-1], str(age) + ' Gyr', fontsize=7, color='white', alpha=0.4)

ax3.set_xlabel('Delta_nu [microHz]')
ax3.set_ylabel('delta_nu_02 [microHz]')
ax3.set_title('C-D Diagram (Asteroseismic HR Diagram)')
ax3.legend(fontsize=6, facecolor='#0a0a0a', edgecolor='#333', labelcolor='white',
          loc='upper left')

# --- Panel 4: Mode power spectrum ---
ax4 = axes[1, 1]
np.random.seed(200)
freq = np.linspace(1500, 4200, 20000)

# Generate synthetic power spectrum with peaks at p-mode frequencies
power = 0.5 * np.random.exponential(1, len(freq))  # noise background

# Add p-mode peaks for l=0,1,2,3
for ell in [0, 1, 2, 3]:
    for n in range(10, 28):
        nu_peak = Delta_nu * (n + ell / 2.0 + epsilon) - D_0 * ell * (ell + 1)
        # Mode power with Gaussian envelope centered at ~3000 microHz
        amplitude = 100 * np.exp(-0.5 * ((nu_peak - 3100) / 600)**2)
        # Lorentzian line profile with width ~1 microHz
        linewidth = 1.0 + 2.0 * ((nu_peak - 3000) / 2000)**2
        power += amplitude * linewidth**2 / ((freq - nu_peak)**2 + linewidth**2)

ax4.plot(freq, power, color='#f43f5e', linewidth=0.3, alpha=0.8)

# Smoothed envelope
from numpy import convolve
kernel = np.ones(100) / 100
power_smooth = convolve(power, kernel, mode='same')
ax4.plot(freq, power_smooth * 5, color='#fbbf24', linewidth=1.5, label='Smoothed envelope')

ax4.set_xlabel('Frequency [microHz]')
ax4.set_ylabel('Power spectral density')
ax4.set_title('Synthetic p-mode Power Spectrum')
ax4.set_ylim(0, 150)
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("ECHELLE DIAGRAM AND FREQUENCY ANALYSIS")
print("=" * 60)
print("\n--- Asymptotic Parameters ---")
print("  Delta_nu = " + str(Delta_nu) + " microHz")
print("  D_0 = " + str(D_0) + " microHz")
print("  delta_nu_02 = 6 * D_0 = " + str(round(6*D_0, 1)) + " microHz")
print("  epsilon = " + str(epsilon))
print("\n--- Physical Interpretation ---")
print("  Delta_nu probes mean density (sound travel time)")
print("  delta_nu_02 probes core conditions (age indicator)")
print("  Echelle diagram reveals near-degeneracy of l and l+2")
print("\n--- Solar Sound Crossing Time ---")
R_sun = 6.957e8  # m
cs_mean = 200e3   # m/s approximate
tau_cross = 2 * R_sun / cs_mean
Dnu_calc = 1e6 / tau_cross  # microHz
print("  tau_cross = 2R/c_s = " + str(round(tau_cross, 0)) + " s")
print("  Delta_nu = 1/tau_cross = " + str(round(Dnu_calc, 1)) + " microHz")
print("=" * 60)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Nuclear Fusion Next: Solar Neutrinos →

Helioseismology

3.1 Pressure Modes (p-modes)

Derivation 1: Adiabatic Wave Equation

3.2 Gravity Modes (g-modes)

Derivation 2: Brunt-Vaisala Frequency

3.3 Asymptotic Theory (Tassoul)

Derivation 3: Tassoul's Asymptotic Formula

3.4 Rotational Splitting

Derivation 4: Rotational Frequency Splitting

3.5 Sound Speed Inversions and the Duvall Law

Derivation 5: The Duvall Law

Historical Context

Numerical Simulation

Helioseismology: Mode Frequencies, Ray Paths, Sound Speed, and Rotation

Full Derivation: Tassoul's Asymptotic Formula

Step 1: The Radial Wave Equation

Step 2: JWKB Quantization Condition

Step 3: Low-Degree Limit

Step 4: Large Frequency Separation

Step 5: Including the \(\ell\)-Dependence

Step 6: Small Frequency Separation

Full Derivation: Rotational Frequency Splitting

Step 1: Rotation as a Perturbation

Step 2: Rotational Kernel

Step 3: 1D Inversion for Radial Rotation

Step 4: Latitude-Dependent Splitting

Key Results from Rotational Inversions

\(\ell\)-\(\nu\) Diagram (Echelle Diagram)

Extended Simulation: Echelle Diagram and Frequency Separations

Echelle Diagram, Frequency Separations, C-D Diagram, and Synthetic Power Spectrum