Sound crossing a layered star

Slow down one part of the journey

Keep the constant-density spherical cavity and its pressure-release boundary. Let the core extend from the center to radius $a$ , with sound speed $c_1$ . The envelope extends to $R$ , with sound speed $c_2$ . Both speeds are positive. Changing sound speed here changes the local compressibility; density remains fixed.

Use dimensionless radius $x=r/R$ , interface position $b=a/R$ and speed ratio $v=c_1/c_2$ . Time is measured in $R/c_2$ , and ordinary frequency in $c_2/R$ . These choices separate the structure from an arbitrary overall size and speed. A dimensionless frequency of one does not mean one hertz unless the chosen physical scales make it so.

The center-to-boundary acoustic travel time is the sum of the times spent in the two layers. This suggests an approximate high-order frequency spacing. It is a starting estimate, not a requirement that every neighboring low-order pair has exactly that separation.

\tau=\frac{a}{c_1}+\frac{R-a}{c_2},\qquad \frac{\tau c_2}{R}=\frac{b}{v}+1-b

Each term is distance divided by local wave speed. The dimensionless travel time changes continuously as the interface moves.

\frac{\Delta\nu_{\rm asym}R}{c_2}\simeq\frac{1}{2(b/v+1-b)}

The subscript marks the high-order acoustic estimate. Reflection at an internal discontinuity can disturb individual spacings.

Reading: MESA School 2026, Lab 3: Asteroseismology

Join the waves at the interface

For harmonic motion, pressure must be continuous across an interface with no imposed surface force. Normal gas velocity must also be continuous. The momentum equation therefore makes the pressure gradient divided by density continuous. Since our density is the same on both sides, the pressure gradient itself is continuous.

Within either layer, the transformed variable $\psi=rP$ still obeys a sinusoidal equation. At the interface, both $\psi$ and its radial derivative are continuous. We may choose any nonzero overall normalization, propagate the regular central solution to the interface, and then propagate through the envelope.

If density were allowed to jump as well, the pressure-gradient condition would change. The simple continuity of the derivative of $\psi$ would no longer be the correct general rule. The pressure equation below retains the required density factors.

[P]=0,\qquad\left[\frac{P^\prime}{\rho}\right]=0

Brackets denote the difference between the two sides of the interface. They express continuity of pressure and normal velocity.

-\frac{d}{dr}\left(\frac{r^2}{\rho}P^\prime\right)=\omega^2\frac{r^2}{\rho c^2}P

This pressure Sturm–Liouville equation retains the density factors needed when the density varies. The browser’s controllable model uses constant density.

Reading: MESA School 2026, Lab 3: Asteroseismology

An exact two-layer frequency calculation

Set $R=1$ in the dimensionless equation and define $s=\omega R/c_2$ . The core wavenumber is $s/v$ and the envelope wavenumber is $s$ . Start with a central solution whose transformed derivative is one. At the interface its value and derivative are known analytically.

In a uniform layer, a value and derivative at one end determine the value at the other. The envelope propagation formula below therefore gives the transformed pressure at the outer boundary. Its zeros select the modes. It contains only sines and cosines, so a tangent pole cannot be mistaken for a root.

The travel-time estimate predicts where roots roughly lie. To find several modes reliably, scan a sufficiently fine range, bracket every sign change, refine each bracket, and verify the expected node count. Merely asking a nonlinear solver for several starting guesses can return the same root repeatedly.

\psi(b)=\frac{\sin(sb/v)}{s/v},\qquad\psi^\prime(b)=\cos(sb/v)

These values follow from the regular core solution. The arbitrary normalization has no effect on the allowed frequencies.

F(s)=\psi(b)\cos[s(1-b)]+\frac{\psi^\prime(b)}{s}\sin[s(1-b)]=0

The boundary residual vanishes for a permitted mode. Convert the root to dimensionless ordinary frequency by dividing by

2\pi

Reading: MESA School 2026, Lab 3: Asteroseismology

Calculate one boundary and inspect the gaps

Choose $b=0.37$ and $v=0.4$ . The dimensionless acoustic travel time is 1.555, and its high-order spacing estimate is about 0.322. The exact first three ordinary frequencies are 0.298308, 0.686279 and 0.919733. Their first two gaps are about 0.388 and 0.233, visibly different from one another and from the travel-time estimate.

To put those numbers into a physical example, take a one-meter cavity with envelope sound speed 300 meters per second. Core speed is then 120 meters per second. Multiply each dimensionless frequency by 300 per second: the first three frequencies are about 89.5, 206 and 276 hertz. The difference between the first two is about 116 hertz, while the next difference is about 70.0 hertz.

The wave reflects at the internal boundary as well as interacting with the outer boundary. Some mode shapes place strong pressure variation in the slow core; others give that layer a different weight. One average travel time cannot describe every low-order phase shift. The changing gaps record those extra phase shifts.

import numpy as np
from scipy.optimize import brentq
b, v = 0.37, 0.4
def boundary(s):
    value = np.sin(s*b/v)/(s/v)
    derivative = np.cos(s*b/v)
    return value*np.cos(s*(1-b)) + derivative*np.sin(s*(1-b))/s
first_root = brentq(boundary, 1.0, 2.0)
print(first_root/(2*np.pi))  # 0.298307644...

Reading: MESA School 2026, Lab 3: Asteroseismology

Turn the differential equation into a checked matrix

An arbitrary positive speed profile will not usually have a two-layer analytic solution. Approximate the transformed shape using small linear hat functions on a radial mesh. Multiplying the differential equation by each hat function and integrating gives a stiffness matrix and a weighted mass matrix.

The stiffness matrix measures spatial gradients. The mass matrix weights each layer by inverse sound speed squared. Here $v(x)$ is local sound speed divided by envelope speed: it equals $v$ in the core and one in the envelope. A generalized eigenvalue problem asks which shape can balance these two contributions. The eigenvalue is dimensionless angular frequency squared. Positivity of both quadratic forms keeps squared frequencies positive.

At an interface crossing a mesh cell, split the integral at the actual boundary. Rounding the boundary to a grid point would make the controls jump in steps. The continuum problem has no such snapping, so the numerical model should not introduce it.

For the example above, the independent consistent finite-element reference gives maximum first-eight frequency errors of about 1.10 percent at 64 cells,0.272 percent at 128 cells and 0.0682 percent at 256 cells, compared with the exact transfer roots. The browser uses a separately checked positive lumped discretization. Its 289-point structure test gives a largest first-eight error of 0.179 percent at 192 cells. These are two distinct numerical schemes and evidence sets.

K_{ij}=\int_0^1\phi_i^\prime\phi_j^\prime\,dx,\qquad M_{ij}=\int_0^1\frac{\phi_i\phi_j}{v(x)^2}\,dx,\qquad K\mathbf a=\lambda M\mathbf a

The hat functions

\phi_i

form a local basis.

\mathbf a

contains shape coefficients and

\lambda=(\omega R/c_2)^2

Reading: MESA School 2026, Lab 3: Asteroseismology

Check shapes as well as frequencies

A small frequency error does not prove that a mode was identified correctly. Check that frequencies are ordered, that the nth radial pressure mode has $n-1$ internal nodes, and that the outer pressure and central regularity conditions hold. The discrete equation should also have a small residual.

Different modes are orthogonal under the wave equation’s mass weight. This means their weighted overlap vanishes, not that their curves never cross. The property follows by multiplying each eigenvalue equation by the other mode, subtracting and integrating: the boundary terms vanish, leaving the frequency difference times the weighted overlap.

Checking frequencies, shapes and boundaries gives the displayed nodes a physical meaning. The drawn node rings come from the same eigenfunctions that supply the model frequencies.

\int_0^1\frac{\psi_m(x)\psi_n(x)}{v(x)^2}\,dx=0\qquad(m\ne n)

The weight follows from the constant-density transformed wave equation and is applied to the eigenfunctions.

Reading: MESA School 2026, Lab 3: Asteroseismology

One pitch still leaves many possible interiors

Consider two exact models. The first has interface 0.37 and core speed ratio 0.65. The second has interface 0.60 and core speed ratio about 0.797655. Both have the same fundamental dimensionless frequency, about 0.421445. Their interiors are clearly different.

The third-mode frequencies are about 1.26389 and 1.29771. An additional well-resolved mode could separate these two cases. It still would not prove that only one arbitrary continuous sound-speed profile is possible. Every measured mode supplies a finite constraint, while a freely varying profile contains many degrees of freedom.

Real stellar modeling adds mass, composition, mixing and atmospheric physics. That is why an interior comparison must name its allowed family. A single fixed evolutionary track answers a conditional question about that track. A broad age or composition inference needs an appropriate ensemble and its uncertainties.

Reading: MESA School 2026, Lab 3: Asteroseismology · Lindsay et al. 2025, Precise Asteroseismic Ages for the Helmi Streams

Sound crossing a layered star

Slow down one part of the journey

Join the waves at the interface

An exact two-layer frequency calculation

Calculate one boundary and inspect the gaps

Turn the differential equation into a checked matrix

Check shapes as well as frequencies

One pitch still leaves many possible interiors

Exercises and explanations

Try the prediction

Change a layer and read its modes

Sources and calculation scope