Oscillations

Reading Assignment

Chapter 13 of Stellar-Astrophysics Notes

Radial, Adiabatic Oscillations

Let's start with radial oscillations. We will use the enclosed mass $m$ as the coordinate and we can write

$P(m,t)=P_{0}(m)+P_{1}(m,t)=P_{0}(m)\left[1+p(m)e^{i\omega t}\right]$

$r(m,t)=r_{0}(m)+r_{1}(m,t)=r_{0}(m)\left[1+x(m)e^{i\omega t}\right]$

$\rho (m,t)=\rho _{0}(m)+\rho _{1}(m,t)=\rho _{0}(m)\left[1+d(m)e^{i\omega t}\right]$

and write the equation of motion of a shell of fluid (hydrostatic equilibrium if no oscillation)

${\ddot {r}}=-{\frac {Gm}{r^{2}}}-4\pi r^{2}{\frac {\partial P}{\partial m}}.$

Substituting the perturbed expressions and keeping only the first-order terms gives

$-\omega ^{2}r_{0}x=-{\frac {Gm}{r_{0}^{2}}}(-2x)-4\pi r_{0}^{2}(2x){\frac {\partial P_{0}}{\partial m}}-4\pi r_{0}^{2}\left(P_{0}{\frac {\partial p}{\partial m}}.+p{\frac {\partial P_{0}}{\partial m}}\right)$

If we use the unperturbed equation and define $g_{0}=Gm/r_{0}^{2}$

$-\omega ^{2}r_{0}x=g_{0}(2x)+(2x+p)g_{0}-4\pi r_{0}^{2}P_{0}{\frac {\partial p}{\partial m}}$

and

$\omega ^{2}r_{0}x+g_{0}(p+4x)={\frac {P_{0}}{\rho _{0}}}{\frac {\partial p}{\partial r_{0}}}$

where we used the zeroth-order continuity equation

${\frac {\partial r}{\partial m}}={\frac {1}{4\pi r^{2}\rho }}.$

Now if we linearize this equation we get

$r_{0}{\frac {\partial x}{\partial m}}+x{\frac {\partial r_{0}}{\partial m}}={\frac {1}{4\pi r^{2}\rho }}(-2x-d)$

$r_{0}{\frac {\partial x}{\partial m}}={\frac {1}{4\pi r^{2}\rho }}(-3x-d)$

and

$r_{0}{\frac {\partial x}{\partial r_{0}}}=-3x-d.$

We now have two equations but three unknowns $p,d$ and $x$ . The final ingredient is the equation of state which for an adiabatic gas gives

$p=\gamma _{\mathrm {ad} }d=-3\gamma _{\mathrm {ad} }x-\gamma _{\mathrm {ad} }{\frac {\partial x}{\partial r_{0}}},$

so now we can write the first equation just in terms of $x$ :

${\frac {\partial }{\partial r_{0}}}\left(\gamma _{\mathrm {ad} }{\frac {\partial x}{\partial r_{0}}}\right)+{\frac {4}{r_{0}}}{\frac {\partial }{\partial r_{0}}}\left(\gamma _{\mathrm {ad} }x\right)-{\frac {\rho _{0}g_{0}}{P_{0}}}\gamma _{\mathrm {ad} }{\frac {\partial x}{\partial r_{0}}}+{\frac {\rho _{0}}{P_{0}}}\left[{\frac {g_{0}}{r_{0}}}\left(4-3\gamma _{\mathrm {ad} }\right)+\omega ^{2}\right]x^{2}=0$

We can rewrite this and a Sturm-Liouville equation

${\hat {\mathcal {L}}}x=-{\frac {1}{\rho _{0}r_{0}^{4}}}{\frac {\partial }{\partial r_{0}}}\left(\gamma _{\mathrm {ad} }P_{0}r_{0}^{4}{\frac {\partial x}{\partial r_{0}}}\right)-{\frac {1}{r_{0}\rho _{0}}}\left\{{\frac {\partial }{\partial r_{0}}}\left[\left(3\gamma _{\mathrm {ad} }-4\right)P_{0}\right]\right\}x=\omega ^{2}x.$

This is an eigenvalue equation for the function $x$ . There are a countably infinite number of eigenvalues, all real. We will order this such that $\omega _{n}^{2}<\omega _{n+1}^{2}$ . The smallest eigenvalue $\omega _{0}^{2}$ will correspond to an eigenfunction $x_{0}(r_{0})$ without any nodes and similarly $\omega _{n}^{2}$ will correspond to an eigenfunction with $n$ nodes. Furthermore, these eigenfunctions will be orthogonal with the weighting $\rho r_{0}^{4}$ over the interval $[0,R]$ .

Let's integrate all three terms times $\rho r_{0}^{4}$ from $0$ to $R$ for the eigenfunction $x_{0}$ to get

$\left.-\left(\gamma _{\mathrm {ad} }P_{0}r_{0}^{4}{\frac {\partial x_{0}}{\partial r_{0}}}\right)\right|_{0}^{R}-\int _{0}^{R}r_{0}^{3}{\frac {\partial }{\partial r_{0}}}\left[\left(3\gamma _{\mathrm {ad} }-4\right)P_{0}\right]x_{0}dr_{0}=\omega _{0}^{2}\int _{0}^{R}x_{0}\rho _{0}r_{0}^{4}dr_{0}.$

and

$\omega _{0}^{2}=-\int _{0}^{R}r_{0}^{3}{\frac {\partial }{\partial r_{0}}}\left[\left(3\gamma _{\mathrm {ad} }-4\right)P_{0}\right]x_{0}dr_{0}\left[\int _{0}^{R}x_{0}\rho _{0}r_{0}^{4}dr_{0}\right]^{-1}.$

If we take the adiabatic index to be constant and use the equation of hydrostatic equilibrium, we get

$\omega _{0}^{2}=\left(3\gamma _{\mathrm {ad} }-4\right)\int _{0}^{R}r_{0}^{3}\rho _{0}g_{0}x_{0}dr_{0}\left[\int _{0}^{R}x_{0}\rho _{0}r_{0}^{4}dr_{0}\right]^{-1}.$

Because $x_{0}$ has no nodes, $\omega _{0}^{2}$ has the same sign as $3\gamma _{\mathrm {ad} }-4$ ; therefore $\gamma _{\mathrm {ad} }>4/3$ means that $\omega ^{2}$ for all of the modes are positive (all stable) and for $\gamma _{\mathrm {ad} }<4/3$ at least $\omega _{0}^{2}<0$ and other might be too (at least one unstable mode).

Let's write everything in terms of dimensionless variables

$q={\frac {m}{M}},{\hat {r}}={\frac {r_{0}}{R}},{\hat {\rho }}={\frac {\rho }{\bar {\rho }}}$

to get

$\omega ^{2}=\left(3\gamma _{\mathrm {ad} }-4\right){\frac {GM}{R^{3}}}\int _{0}^{1}{\hat {r}}{\hat {\rho }}qxd{\hat {r}}\left[\int _{0}^{1}x{\hat {\rho }}{\hat {r}}^{4}d{\hat {r}}\right]^{-1}.$

We see that $P^{2}M/R^{3}$ is a constant for homologous stars. For polytropes one can also get the frequencies of the overtones

$\omega ^{2}={\frac {4\pi G\gamma _{\mathrm {ad} }\rho _{c}}{n+1}}\Omega _{n}^{2}$

where $\Omega _{n}^{2}$ is a dimensionless eigenvalue that depends only on $\gamma _{\mathrm {ad} }$ and $n$ .

Can we estimate the frequencies in general for modes with many nodes? Let's use a variational principle and guess the eigenfunction to be

$x_{n}=\sin \left(\pi {\frac {n+1}{2}}{\hat {r}}\right).$

If we substitute this into the homologous frequency equation and integrate the denominator by parts twice we will get

$\omega ^{2}\propto \left({\frac {n+1}{2}}\right)^{2}~\mathrm {so} ~\omega \approx (n+1)\omega _{0}.$

Radial, Non-Adiabatic Oscillations

The key idea here is that for the modes to be sustained, they must not be strictly adiabatic! We can imagine the mode is like an engine and work must be done to drive the oscillation so heat must enter the gas at high temperature. For this to work the gas must be more opaque at the high temperature, high density part of the oscillation. We have

$\kappa \propto {\frac {\rho }{T^{3.5}}}$

and

$T\propto \rho ^{\gamma _{\mathrm {ad} }-1}$

so

$\kappa \propto \rho ^{(9-7\gamma _{\mathrm {ad} })/2}.$

Therefore, the opacity will be smaller at the high density if $\gamma _{\mathrm {ad} }>9/7$ , and this is typically the case.

However, an exception is a partially ionised layer. The compression ionises the layer but the temperature does not increase. Heat is absorbed during the compression and released during the expansion. This is called the $\kappa$ mechanism.

H-ionization zone at 10,000-15,000 K: HI to HII and He I to He II.

He-ionization zone at 40,000 K: He II to He III

For $T_{\mathrm {eff} }\approx 7500\mathrm {K}$ these zone are too close to surface not much mass!

For $T_{\mathrm {eff} }\approx 5500\mathrm {K}$ these zone are deeper and the fundamental can be excited. There is in instability strip around effective temperatures between 5,500 and 7,500 K.

Non-Radial, Adiabatic Oscillations

Now things get more complicated, but we will make some simplifications to make our lives easier. First, the equations

$\nabla ^{2}\Phi =4\pi G\rho$

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

$\rho \left({\frac {\partial }{\partial t}}+\mathbf {v} \cdot \nabla \right)=-\nabla P-\rho \nabla \Phi .$

Remember that the unperturbed solution is just a function of radius $r$ . We can follow the changes in the fluid in two ways, the Eulerian (I'm standing at a point) and Lagrangian (I'm moving with the fluid). We will use the prime to denote the Eulerian perturbation and $\delta$ to denote the Lagrangian perturbation. We can relate them by

$\delta \rho =\rho '+\mathbf {x} \cdot \nabla \rho$

where $\mathbf {x} (\mathbf {r} ,t)$ is the perturbation of the fluid from its equilibrium position. We can also think of changes in time as measured by someone in the fluid and someone standing still

${\frac {d}{dt}}={\frac {\partial }{\partial t}}+\mathbf {v} \cdot \nabla .$

We will connect the perturbation of the pressure to that of the density

${\frac {\delta P}{P}}=\gamma _{\mathrm {ad} }{\frac {\delta \rho }{\rho }}.$

We can now write perturbed Euler equation

$\rho {\frac {d^{2}\mathbf {x} }{\partial t^{2}}}=-\nabla P'-\rho \nabla \Phi '-\rho '\nabla \Phi .$

and we can connect the density perturbation to the displacement field through the continuity equation

${\frac {\delta \rho }{\rho }}=-\nabla \cdot \mathbf {x} .$

We would also need to calculate the perturbation to the gravitational potential through the Poisson equation, but we will make the Cowling approximation in which we take $\Phi '=0$ and we write

$\mathbf {x} (\mathbf {r} ,t)=\mathbf {x} \exp(i\sigma t)$

and get equations for the three components of the displacement

$\sigma ^{2}x_{r}={\frac {\partial }{\partial r}}\left({\frac {P'}{\rho }}\right)-A{\frac {\gamma _{\mathrm {ad} }P}{\rho }}\nabla \cdot \mathbf {x}$

$\sigma ^{2}x_{\theta }={\frac {\partial }{\partial \theta }}\left({\frac {1}{r}}{\frac {P'}{\rho }}\right)$

$\sigma ^{2}x_{\phi }={\frac {1}{\sin \theta }}{\frac {\partial }{\partial \phi }}\left({\frac {1}{r}}{\frac {P'}{\rho }}\right)$

where

$A={\frac {d\ln \rho }{dr}}-{\frac {1}{\gamma _{\mathrm {ad} }}}{\frac {d\ln P}{dr}}.$

Notice that $A=0$ if the unperturbed star is adiabatic.

We can solve the two angular equations by using spherical harmonics

$\mathbf {x} (r,\theta ,\phi )=\left[x_{r}(r)\mathbf {e} _{r}+x_{t}(r)\mathbf {e} _{\theta }{\frac {\partial }{\partial \theta }}+x_{t}(r)\mathbf {e} _{\phi }{\frac {1}{\sin \phi }}{\frac {\partial }{\partial \phi }}\right]Y_{lm}(\theta ,\phi )$

and

${\frac {P'(\mathbf {r} )}{\rho }}={\frac {P'(r)}{\rho }}Y_{lm}(\theta ,\phi )$

where the transverse component $x_{t}(r)$ is related to $P'/\rho$ by

$x_{t}(r)={\frac {1}{\sigma ^{2}}}{\frac {1}{r}}{\frac {P'(r)}{\rho }}$

so we can focus on $x_{r}$ and $x_{t}$ yielding the following two equations

$r{\frac {dx_{r}}{dr}}=\left({\frac {k_{t}^{2}gr}{S_{l}^{2}}}-2\right)x_{r}+r^{2}k_{t}^{2}\left(1-{\frac {\sigma ^{2}}{S_{l}^{2}}}\right)x_{t}$

and

$r{\frac {dx_{t}}{dr}}=\left(1-{\frac {N^{2}}{\sigma ^{2}}}\right)x_{r}+\left({\frac {r}{g}}N^{2}-1\right)x_{t}$

where we have defined the Brunt-Vaisala frequency $N$ , the Lamb frequency $S_{l}$ and the transverse wavenumber $k_{t}$ as

$N^{2}\equiv -Ag=-g\left({\frac {d\ln \rho }{dr}}-{\frac {1}{\gamma _{\mathrm {ad} }}}{\frac {d\ln P}{dr}}\right)$

$S_{l}^{2}\equiv {\frac {l(l+1)}{r^{2}}}{\frac {\gamma _{\mathrm {ad} }P}{\rho }}={\frac {l(l+1)}{r^{2}}}v_{s}^{2}$

and

$k_{t}^{2}=k_{l}^{2}={\frac {l(l+1)}{r^{2}}}={\frac {S_{l}^{2}}{v_{s}^{2}}}.$

The Brunt-Vaisala frequency is simply the oscillation frequency of a blob displaced vertically

$N^{2}=g\left[\left({\frac {d\rho }{dP}}\right)_{*}-\left({\frac {d\rho }{dP}}\right)_{S}\right]$

and the Lamb frequency is the frequency of sound waves travelling transversely.

Unlike the equation for purely radial modes, this is not a Sturm-Liouville equation, so we have no guarantee that the value of $\sigma ^{2}$ is bounded from below. All that we do know is that $\sigma ^{2}$ is real and that the solutions will be orthogonal.

We can talk about modes with $m=0$ which will be constant on lines of constant latitude ( $\theta$ ). These are called zonal modes. We can also have modes with $m=|l|$ . These are constant on lines of longitude and are called sectoral modes. The other modes are called tesseral modes.

To make some further progress we will take a WKB approximation in which the radial variation is proportional to $\exp \left[ik_{r}(r)r\right]$ in the limit where $k_{r}r\gg 1$ . This yields a dispersion relation for the radial wavenumber

$k_{r}^{2}={\frac {k_{t}^{2}}{\sigma ^{2}S_{l}^{2}}}\left(\sigma ^{2}-N^{2}\right)\left(\sigma ^{2}-S_{l}^{2}\right).$

For a wave that is oscillatory in time we have $\sigma ^{2}>0$ . Furthermore for it to be oscillatory in the radial direction as well, we must have $k_{r}^{2}>0$ ; therefore, we can have waves for values of $\sigma$ outside the range in between $N$ and $S_{l}$ . We have two types of waves those with $\sigma$ greater than both $N$ and $S_{l}$ and those with $\sigma$ less than both $N$ and $S_{l}$ . Usually we have $N<S_{l}$ . In this case we call the first $P$ modes or pressure modes and the second we call $g$ modes or gravity modes.

If we assume that $\sigma ^{2}\gg N^{2},S_{l}^{2}$ we have

$\sigma _{P}^{2}\approx {\frac {k_{r}^{2}+k_{t}^{2}}{k_{t}^{2}}}S_{l}^{2}=\left(k_{r}^{2}+k_{t}^{2}\right)v_{s}^{2}$

and $\sigma ^{2}\ll N^{2},S_{l}^{2}$

$\sigma _{g}^{2}\approx {\frac {k_{t}^{2}}{k_{r}^{2}+k_{t}^{2}}}N^{2}.$

The size of the perturbations for the two types of modes are given by

$\left|{\frac {x_{r}}{x_{t}}}\right|\sim \left\{{\begin{array}{ll}rk_{r}&(P~\mathrm {modes} )\\l(l+1)/(rk_{r})&(g~\mathrm {modes} )\end{array}}\right.$

We can approximate the number of radial nodes as

$n\approx {\frac {1}{\pi }}\int _{0}^{R}k_{r}(r)dr$

and the frequencies as

$\sigma _{P}\approx n\pi \left(\int _{0}^{R}{\frac {dr}{v_{s}}}\right)^{-1}$

and

$P_{g}={\frac {2\pi }{\sigma _{g}}}\approx n{\frac {2\pi ^{2}}{\left[l(l+1)\right]^{1/2}}}\left(\int _{0}^{R}{\frac {N}{r}}dr\right)^{-1}.$

Check out the latest Doppler images of the Sun at http://jsoc.stanford.edu/data/hmi/images/latest/

Assignment

You will use the inlist called astero_adipls for this assignment. This model evolves a 1.2 solar mass star from the ZAMS until the hydrogen is nearly exhausted in the centre. I would like you to plot the sound speed (profile column csound) and Brunt frequency (brunt_N2) as a function of radius for the initial stage (ZAMS) , halfway through the evolution and at the final stage. Identify the convectively unstable regions on the plots. Plot the Lamb frequency for l=1,2,3. Identify the turning points of p-mode with ω=3000μHz for different values of angular momentum. What about 10,000μHz? Similarly what are the turning points of a 2000μHz g-mode?

Plot the structure of the first several modes (l=0, n=0 to 6) for the final model. The results are in the file ttt.adipls.prt.