Lecture IV: Modeling Stellar Pulsation A pulsating star is not in hydrostatic equilbrium. For example ρd2r/dt2 = -GMrρ/r2 – dP/dr. Mass continuity equation still holds. Energy equation: dE/dt + PdV/dt + dL/dm = 0, where L(r) = -4πr24σ/3κ . dT4/dm ρ(r) = 1/V(r), P = P(ρ,T), E=E(ρ,T), κ=κ(ρ,T).

Modeling Stellar Pulsation Boundary Cnditions: L0=Lcons., dr/dt)0 = 0. Psurface = 0. Tsurface = f(Tef) ie. a grey solution to the equationof radiative transfer. 1D radiative codes. Now there are “numerical recipes” to model time dependent turbulent convection.

Linear Models Assume displacement from equilbrium, δr, are small. Write variables as P = P0 + δP, r = r0 + δr, ρ0 + δρ etc. Expand pulsation equations and drop second order terms. This is linear stellar pulsation. Assume δr = |δr|eiωt, solve resulting eigenvalue problem. Leads to linear periods and growth rates ie. Whether a given perturbation is stable or will continue to grow. Can investigate boundaries of “instability strip” with such a technique.

Non-Linear Models Write diferential equations as diference equations over a computational grid covering the star. Zones 1,……,N, with interfaces 0,1,….N+1. Extensive variables r, velocity, vr, luminosity, Lr, defined at zone interfaces. Intensive variables defined at zone centers, T, ρ, P, κ etc. Sometimes may need to extrapolate intensive/extensive variables to zone interface/centers. Time mesh: tn+1 = tn + Δtn+1/2,tn+1/2 – tn-1/2 = Δtn, Δtn = ½(Δtn-1/2 + Δtn+1/2).

Non-Linear Models Momentum equation: vn+1/2(I) = vn-1/2(I) – Δtn(GM(I)/rn(I)2 + 4π(rn(I))2/ΔM(I)[Pn(I) – Pn(I-1) + Qn-1/2(I) – Qn+1/2(I-1)]) Leads to a matrix equation Ax=d to be solved for the increments to the physical variables at each time step. Q: Artifical vsicosity. Field in its own right.

Pulsation Modeling Linear model to find set of L,M, X,Z,Tef. Also get eigenvector showing ampltide of rafial displacement. Non-linear model with an initial “kick” scaled by linear eigenvector for that model Continue pulsation until amplitude increase levels of: several hundred cycles, maybe 1-2 hours on a modern fast PC. Need opacity tables, equation of state (usually Saha). Result is a nonlinear full amplitude variation of L with T. Stellar atmosphere converts this to magnitude and color. Compare with observations via Fourier analysis. This is for radial oscillations. No time dependent code to model non-radial oscillations exists.

Non-Radial Oscillations Expand perturbatin δr in terms of spherical harmonics, specified by 3 numerbs, n, l, m. δr = R(r)Y(θ,φ): n is for the radial part, l, m the angular part. l=m=0, pulsation purely radial. l=0,1,2,,,n-1 and m=-l+1,-l+2,….l-1 With l,m non-zero need to worry about Poisson’s equation as well. n: number of nodes radially outward from Sun’s center. m: number of nodes found around the equator. l: number of nodes found around the azimuth (great circle through the poles) Hard mathematical/numerical problem. P-modes: pressure is the restoring force, G modes: gravity is the restoring force.

Helioseismology Sun is a non-radial oscillator. Modes with periods between 3 an d8 minutes – five minute oscillations are p modes: l going from 0 to 1000. Modes with longer periods – about 160 minutes could be g modes: l ~1-4. Comparison of observed and theoretical frequencies can be used to calibrate solar models: helioseismology. Can reveal the depth of the solar convection zone, plus rotation and composition of the outer layers of the Sun.

One Zone Models Central point mass of mass M. At a radius R is a thin spherical shell, mass m. There is a pressure P in this shell which provides support against gravity. Newton’s second law: md2R/dt2 = -GMm/R2 + 4πR2P In equilbrium, GMm/R02 = 4πR02P0 Linearize: R = R0+δR, P = P0+δP Insert into momentum equation, linearize, keep only first powers of δs and use d2R0/dt2 = 0 to give

One Zone Models md2(δR)/dt2 = 2GMm(δR)/R03 + 8πR0P0(δR) + 4πR02δP Adiabatic oscillations:PVγ = const. Linearized version: δP/P0 = -3γδR/R0 Hydrostatic equilbrium means 8πR0P0 = 2GMm/R03. The the linearized equation for δR is d2(δR)/dt2 = -(3γ – 4)GM(δR)/R03 Simple Harmonic Motion, δR = Asin(ωt) with ω2=(3γ-4)GM/R03 Since, the pulsation period, Π = 2π/ω, we have Π = 2π/(√[4πGρ0(3γ-4)]), the period mean density theorem.