4.1: Pulsational Stability of White Dwarfs

By now, I hope the reader has been impressed by the wide range of problems which can be dealt with by the virial theorem. Some of the problems mentioned in the last chapter indicate the type of insight which can be achieved through use of virial theorem, however, the type of objects which were explicitly discussed, notably normal stars, are currently judged by the naive to be well understood. In order to illustrate the power of this remarkable theorem, I cannot resist discussing some objects about which even less is known. We shall see that at least one commonly held tenant of stellar structure, while leading to nearly the correct numerical result, is conceptually wrong.

During the 1960's advances in observational astronomy presented problems requiring theoreticians to postulate the existence of a wide range of objects previously considered only of academic interest. These terms, like supermassive stars, neutron stars, and black holes became 'household' words in the literature of astrophysics. Many of the objects were clearly so condensed as to require the application of the General Theory of Relativity or some other gravitational theory for their description. When one postulates the existence of a "new" object it is always wise to subject that object to a stability analysis. This is particularly important for highly collapsed objects as the time scale for development of the instability will be very short. Since general relativistic effects can usually be viewed as an effective increase in the gravitational force, one would expect its presence to decrease the stability of objects in which it is important. What came as a surprise is the importance of these effects where one would normally presume them to be of little or no importance.

Apparently inspired by a comment of R. P. Feynman in 1963, W. A. Fowler noted that effects of general relativity would lead to previously unexpected instabilities in supermassive stars¹.^† Noting that the conditions for this instability also exist in massive white dwarfs, Chandrasekhar and Tooper² showed by means of rather detailed calculations that a white dwarf would become unstable when its radius shrank to about 246 Schwarzschild radii or on the order of 1000 km. This corresponds to a mass about 1.5% below the well-known Chandrasekhar limiting mass for degenerate objects. During the next 15 years, this instability received a great deal of attention and I will not attempt to fully recount it here. Rather, let us examine with the aid of hindsight and the virial theorem, how this result could be anticipated without the need of detailed calculations.

One can see that the stability analysis coupled with the post-Newtonian form of the virial theorem given in Chapter II [Equation \ref{2.4.15}] would serve as the basis for investigating this effect. However the estimation or calculation of the relativistic terms on the right hand side of Equation \ref{2.4.15} is extremely difficult. Instead, by assuming spherical symmetry we may start with the spherically symmetric Equation of motion given by Meltzer and Thorne³ as did Fowler⁴ and follow the formalism of Chapter III. Thus

\[\begin{gathered}
\mathrm{y \frac{\mathrm{d}}{\mathrm{dt}}\left(y \frac{d r}{d t}\right)=-\frac{1}{\rho} \frac{\mathrm{dP}}{\mathrm{dr}}\left[\frac{1+\frac{y^2}{\mathrm{c}^2}\left(\frac{\mathrm{dr}}{\mathrm{dt}}\right)^2-\frac{2 \mathrm{Gm}(\mathrm{r})}{\mathrm{rc}^2}}{1+\frac{\mathrm{P}}{\rho \mathrm{c}^2}}\right]-\frac{\mathrm{Gm}(\mathrm{r})}{\mathrm{r}^2}-\frac{4 \pi \mathrm{G} \rho \mathrm{r}}{\mathrm{c}^2}}, \\[4pt]
\text{where}\quad y=\frac{\left(\rho+\mathrm{P} / \mathrm{c}^2\right)}{\rho_0} .
\end{gathered}\label{4.1.1}\]

If we confine our attention to objects nearly in equilibrium, no large scale radial motions can exist. Thus the term involving (dr/dt)² can be neglected and Equation \ref{4.1.1} becomes

\[\mathrm{y}^2 \frac{\mathrm{d}^2 \mathrm{r}}{\mathrm{dt}^2}=\frac{1}{\rho} \frac{\mathrm{dP}}{\mathrm{dr}}\left[\frac{1-2 \mathrm{Gm}(\mathrm{r}) / \mathrm{rc}^2}{1+\mathrm{P} / \rho \mathrm{c}^2}\right]-\frac{\mathrm{Gm}(\mathrm{r})}{\mathrm{r}^2}-\frac{4 \pi \mathrm{GPr}}{\mathrm{c}^2},\label{4.1.2}\]

or in the post-Newtonian approximation (i.e.. keeping only terms of the order 1/c²)

\[\mathrm{\rho y^2 \frac{\mathrm{d}^2 \mathrm{r}}{\mathrm{dt}^2}=\frac{\mathrm{dP}}{\mathrm{dr}}\left[1-\frac{\mathrm{P}}{\rho \mathrm{c}^2}-\frac{2 \mathrm{Gm}(\mathrm{r})}{\mathrm{rc}^2}-\cdots-\right]=\frac{\mathrm{Gm}(\mathrm{r}) \rho}{\mathrm{r}^2}-\frac{4 \pi \mathrm{GP} \rho \mathrm{r}}{\mathrm{c}^2}}.\label{4.1.3}\]

†

See Fricke⁹ who also uses a post-Newtonian virial approach to this problem.

In hydrostatic equilibrium, dP/dr is given by Oppenheimer-Volkoff as

\[\frac{\mathrm{dP}}{\mathrm{dr}}=-\frac{\mathrm{G}\left(\rho+\mathrm{P} / \mathrm{c}^2\left)\lfloor\mathrm{m}(\mathrm{r})+4 \pi \mathrm{r}^3 \mathrm{P} / \mathrm{c}^2\right\rfloor\right.}{\mathrm{r}^2\left(1-2 \mathrm{Gm}(\mathrm{r}) / \mathrm{rc}^2\right)},\label{4.1.4}\]

or

\[\frac{\mathrm{dP}}{\mathrm{dr}}=-\frac{\mathrm{Gm}(\mathrm{r}) \rho}{\mathrm{r}^2}+\frac{1}{\mathrm{c}^2}\left[\frac{2 \mathrm{G}^2 \mathrm{m}^2(\mathrm{r}) \rho}{\mathrm{r}^4}+\frac{\mathrm{Gm}(\mathrm{r}) \rho}{\mathrm{r}^2}+4 \pi \mathrm{rG} \rho^2\right]+\mathrm{O}\left(\frac{1}{\mathrm{c}^4}\right).\label{4.1.5}\]

If we retain dP/dr explicitly for the expansion of the first term of Equation \ref{4.1.3} all other relativistic terms of Equation \ref{4.1.5} will be of the order 1/c⁴ in the product in Equation \ref{4.1.3}. Thus, Equation \ref{4.1.3} becomes

\[\mathrm{\rho y^2 \frac{\mathrm{d}^2 \mathrm{r}}{\mathrm{dt}^2}=-\frac{\mathrm{dP}}{\mathrm{dr}}-\frac{\mathrm{Gm}(\mathrm{r}) \rho}{\mathrm{r}^2}\left[1+\frac{\mathrm{P}}{\rho \mathrm{c}^2}+\frac{2 \mathrm{Gm}(\mathrm{r})}{\mathrm{rc}^2}+\cdots+\right]-\frac{4 \pi \mathrm{rGP} \rho}{\mathrm{c}^2}}.\label{4.1.6}\]

Now if we again form Lagrange's identity by multiplying by r and integrating over all volume, we get

\[\mathrm{\int_v\left(\rho y^2 r \frac{d^2 r}{d t^2}\right) d V==-\int_0^R 4 \pi r^3 \frac{d P}{d r} d r-\int_v \frac{G m(r) \rho}{r} d V-\int_v \frac{G m(r) \rho}{r c^2}-2 \int_v\left(\frac{G^2 m^2(r)}{r c^2}\right) d V^2-\int_v \frac{4 \pi r G P \rho}{c^2} d V}.\label{4.1.7}\]

The last integral can be integrated by parts,^4.1 so that

\[\int_{\mathrm{v}} \frac{4 \pi \mathrm{rGP} \rho}{\mathrm{c}^2} \mathrm{dV}=\int_{\mathrm{v}} \frac{\mathrm{G}^2 \mathrm{m}^2(\mathrm{r}) \rho}{\mathrm{r}^2 \mathrm{c}^2} \mathrm{dV}.\label{4.1.8}\]

With somewhat less effort the first integral becomes

\[\int_0^{\mathrm{R}} 4 \pi \mathrm{r}^3\left(\frac{\mathrm{dP}}{\mathrm{dr}}\right) \mathrm{dr}=\int_0^{\mathrm{R}} 4 \pi \mathrm{r}^3 \mathrm{dP}=\left.4 \pi \mathrm{r}^3 \mathrm{P}\right|_0 ^{\mathrm{R}}-\int_0^{\mathrm{R}} 12 \pi \mathrm{r}^2 \mathrm{Pdr}=-3 \int_{\mathrm{v}} \mathrm{PdV}.\label{4.1.9}\]

Putting the results of Equation \ref{4.1.8} and Equation \ref{4.1.9} into Equation \ref{4.1.7}, noting that the first term on the right hand side of Equation \ref{4.1.9} is zero, and rewriting the left hand side in terms of a relativistic moment of inertia we get

\[\mathrm{\frac{1}{2} \frac{d^2 I_r}{d t^2}=3 \int_v P d V-\Omega-\frac{1}{c^2} \int_v \frac{G m(r) P}{r} d V-\frac{3}{c^2} \int_v \frac{G^2 m^2(r) \rho}{r^2} d V}.\label{4.1.10}\]

which is equivalent to Equation \ref{2.4.15} of Chapter II for spherical stars but vastly simpler. Although this approach, which is basically due to Fowler⁴, lacks the rigor of the EIH post-Newtonian approach, it does yield the same results for spherical stars nearly in hydrostatic equilibrium. It is worth noting that the relativistic correction terms are of the same mixed energy integrals as those that appear in Equation \ref{2.4.15}.

Taking the variation of Equation \ref{4.1.10} as we did in Chapter III, we have

\[\frac{1}{2} \frac{\mathrm{d}^2}{\mathrm{dt}^2}\left(\delta \mathrm{I}_{\mathrm{r}}\right)=3 \delta\left[\int_{\mathrm{v}} \mathrm{PdV}\right]-\delta \Omega-\frac{1}{\mathrm{c}^2} \delta\left[\int_{\mathrm{v}} \frac{\mathrm{Gm}(\mathrm{r}) \mathrm{P}}{\mathrm{r}} \mathrm{dV}\right]-\frac{3}{\mathrm{c}^2} \delta\left[\int_{\mathrm{v}} \frac{\mathrm{G}^2 \mathrm{m}^2(\mathrm{r}) \rho}{\mathrm{r}^2} \mathrm{dV}\right] .\label{4.1.11}\]

As before, let us suppose that the variation of these quantities results from a variation of the independent variable δr. Further assume that \(\delta \mathrm{r} / \mathrm{r}=\xi_{\mathrm{o}} \mathrm{e}^{\mathrm{i} \sigma t}\) where \(\xi_{\mathrm{o}}\) is constant and the variation is adiabatic. Since we can write the internal heat energy density as \(\left(\Gamma_1-1\right) \mathrm{u}=\mathrm{P}\), the first term becomes

\[3 \delta \int_{\mathrm{v}} \mathrm{PdV}=3<\Gamma_1-1>\delta \mathcal{U}.\label{4.1.12}\]

Equation \ref{3.2.10} (i.e. Chapter III) leads to

\[\delta \Omega=-\xi \Omega.\label{4.1.13}\]

The variation of the relativistic correction terms can be computed as follows: Let

\[\Omega_1=\frac{3}{2} \int_0^{\mathrm{M}} \frac{\mathrm{G}^2 \mathrm{m}^2(\mathrm{r})}{\mathrm{r}^2 \mathrm{c}^2} \mathrm{dm}(\mathrm{r}).\label{4.1.14}\]

so that the variation of the last term in Equation \ref{4.1.11} is

\[2 \delta \Omega_1=\frac{3}{\mathrm{c}^2} \int_0^{\mathrm{M}} \mathrm{G}^2 \mathrm{m}^2(\mathrm{r}) \delta\left(\frac{1}{\mathrm{r}^2}\right) \mathrm{dm}(\mathrm{r})=-2 \xi \Omega_1 .\label{4.1.15}\]

It is convenient (particularly for the relativistic terms) to normalize by the dimensionless quantity (2GM/Rc²). Thus,

\[\Omega_1=\frac{3}{2} \mathrm{Mc}^2 \int \frac{1}{4}\left(\frac{2 \mathrm{GM}}{\mathrm{Rc}^2}\right)^2\left(\frac{\mathrm{m}^2(\mathrm{r})}{\mathrm{M}}\right)\left(\frac{\mathrm{R}}{\mathrm{r}}\right)^2 \frac{\mathrm{dm}(\mathrm{r})}{\mathrm{M}}=\frac{3}{8} \mathrm{Mc}^2\left(\frac{2 \mathrm{GM}}{\mathrm{Rc}^2}\right)^2 \int_0^1(\mathrm{q} / \mathrm{x})^2 \mathrm{dq},\label{4.1.16}\]

where the dimensionless variables are \(\mathrm{q}=[\mathrm{m}(\mathrm{r}) / \mathrm{M}]\), and \(\mathrm{x=r / R}\). The remaining terms in 4.1.11 can be normalized in a similar way by making use of the homologous dependence of P. That is

\[\mathrm{P}=\eta \mathrm{Gm}^2(\mathrm{r}) / \mathrm{r}^4,\label{4.1.17}\]

where \(\eta\) is a dimensionless scale factor. Therefore, we can let

\[P_1=\frac{1}{\mathrm{c}^2} \int_{\mathrm{v}} \frac{\mathrm{Gm}(\mathrm{r}) \mathrm{PdV}}{\mathrm{r}^2}=\frac{1}{\mathrm{c}^2} \int_0^{\mathrm{R}} \frac{4 \pi \eta \mathrm{G}^2 \mathrm{m}^3(\mathrm{r}) \mathrm{r}^2 \mathrm{dr}}{\mathrm{r}^5}=\mathrm{Mc}^2\left[\frac{2 \mathrm{GM}}{\mathrm{Rc}^2}\right]^2 \int_0^1 \pi \eta\left(\frac{\mathrm{q}}{\mathrm{x}}\right)^3 \mathrm{dx}.\label{4.1.18}\]

As in Equation \ref{4.1.16}, the integral in Equation \ref{4.1.18} is dimensionless and determined by the equilibrium model. Thus the remaining Equation is

\[\delta P_1=-2 \xi P_1.\label{4.1.19}\]

Replacing P with \(\mathrm{u}\left(\Gamma_1-1\right)\) as with the first term and letting

\[\mathcal{U}_1=\frac{1}{\mathrm{c}^2} \int_{\mathrm{v}} \frac{\mathrm{uGm}(\mathrm{r})}{\mathrm{r}} \mathrm{dV}=\frac{1}{\left\langle\Gamma_1-1\right\rangle} P_1.\label{4.1.20}\]

Equation \ref{4.1.11} then becomes:

\[\frac{1}{2} \frac{\mathrm{d}^2}{\mathrm{dt}^2}\left(\delta \mathrm{I}_{\mathrm{r}}\right)=3<\Gamma_1-1>\delta \mathcal{U}-\delta \Omega+<\Gamma_1-1>\delta \mathcal{U}_1-2 \delta \Omega_1.\label{4.1.21}\]

Now, since the internal energy \(\mathcal{U}\) is coupled with all other terms including the relativity terms we shall eliminate it in a somewhat different fashion than in Chapter III. Since the total energy must be constant, its variation is zero. Thus

\[\mathrm{\delta E=0=\delta \mathcal{U}-\delta \Omega+\delta \mathcal{U}_1-\delta \Omega_1},\label{4.1.22}\]

and Equation \ref{4.1.21} becomes

\[\frac{1}{2} \frac{\mathrm{d}^2}{\mathrm{dt}^2}\left(\delta \mathrm{I}_{\mathrm{r}}\right)=<3 \Gamma_1-4>\delta \Omega-2<\Gamma_1-1>\delta \mathcal{U}_1-<3 \Gamma_1-5>\delta \Omega_1.\label{4.1.23}\]

Substituting in the variations from Equations Equation \ref{4.1.13}, (4.1.15, and Equation \ref{4.1.19} into Equation \ref{4.1.23} and noting that two time differentiations of the perturbation will give a σ² in the first term, Equation \ref{4.1.23} becomes

\[\mathrm{\sigma^2 \mathrm{I}_{\mathrm{r}} \bar{y}=<3 \Gamma_1-4>\Omega_0-4<\Gamma_1-1>\mathcal{U}_1+2<3 \Gamma_1-5>\Omega_1}.\label{4.1.24}\]

Making one last normalization of Ω0 which for polytropes is

\[\Omega_0=\left(\frac{3}{5-\mathrm{n}}\right) \frac{\mathrm{GM}}{\mathrm{R}}=\left(\frac{1}{5-\mathrm{n}}\right) \frac{3}{2}\left(\frac{2 \mathrm{GM}}{\mathrm{Rc}^2}\right) \mathrm{Mc}^2,\label{4.1.25}\]

and calling the dimensionless integrals in Equation \ref{4.1.16} and Equation \ref{4.1.18}, \(\zeta_1\) and \(\zeta_2\) respectively, Equation \ref{4.1.24} becomes

\[\mathrm{\sigma^2 \mathrm{I}_{\mathrm{r}} \bar{y}=\mathrm{Mc}^2\left(\frac{2 \mathrm{GM}}{\mathrm{Rc}^2}\right)\left[\frac{3}{2(5-\mathrm{n})}<3 \Gamma_1-4>-\left(\frac{2 \mathrm{GM}}{\mathrm{Rc}^2}\right)\left(4 \varsigma_1<\Gamma_1-1>+2 \varsigma_2<5-3 \Gamma_1>\right)\right]}.\label{4.1.26}\]

Since the average relativity factor \(\mathrm{\bar{y}}\) is always positive, this expression can be used as a stability criterion as in Chapter III. That is

\[\mathrm{\left(\frac{R_s}{R}\right)\left(4 \varsigma_1<\Gamma_1-1>+2 \varsigma_2<5-3 \Gamma_1>\right)<\left(\frac{3<3 \Gamma_1-4>}{2(5-n)}\right)}.\label{4.1.27}\]

where we have used the fact that the Schwarzschild radius \(\left(\mathrm{R}_{\mathrm{S}}\right)\) is \(2 \mathrm{GM} / \mathrm{c}^2\). We can now use this to investigate the stability of white dwarfs as they approach the Chandrasekhar limiting mass. As this happens, the Equation of state approaches that of relativistically degenerate electron gas and the internal structure, that of a poly trope of index, n = 3.

As \(\mathrm{n\rightarrow3}\), \(\Gamma_1 \rightarrow 4 / 3\), and the system becomes unstable. Thus let

\[\Gamma_1=\frac{4}{3}+\varepsilon,\label{4.1.28}\]

and Equation \ref{4.1.27} becomes (using Fowler’s⁴ values for \(\zeta_1\) and \(\zeta_2\)),

\[\left.\begin{array}{l}
\frac{9 \varepsilon}{4}-\frac{\mathrm{R}_{\mathrm{S}}}{\mathrm{R}_0}\left(\frac{4}{3} \varsigma_1+2 \varsigma_2\right)=\frac{9}{4} \varepsilon-2.5 \frac{\mathrm{R}_{\mathrm{S}}}{\mathrm{R}_0} \geq 0 \\[4pt]
\text { or } \frac{\mathrm{R}_{\mathrm{S}}}{\mathrm{R}_0}>\frac{1.13}{\varepsilon}
\end{array}\right\}.\label{4.1.29}\]

Thus, if you imagine a sequence of white dwarfs of increasing mass, the value of \((\mathrm{R_o/R_s})\) will monotonically decrease as a result of the mass radius relation for white dwarfs and \((1 / \varepsilon)\) will monotonically increase as the configuration approaches complete relativistic degeneracy. Clearly, the point must come where the system becomes unstable and collapses. However, in order to find that point, we need an estimate of how ε changes with increasing mass. For that we turn to an interesting paper by Faulkner and Gribben⁵ who show^4.2

\[\mathrm{\varepsilon \cong \frac{2 x^{-2}}{3}},\label{4.1.30}\]

where x is the Chandrasekhar degeneracy parameter.

So, our instability condition can be written as

\[\frac{\mathrm{R}_{\mathrm{S}}}{\mathrm{R}_0}>1.7 \overline{\mathrm{x}^2}.\label{4.1.31}\]

All that remains is to estimate an average value of the degeneracy parameter x which we can expect to be much larger than 1. From Chandrasekhar⁶

\[\rho_{\mathrm{e}}=\mathrm{Bx}^3=\left(8 \pi \mathrm{m}_{\mathrm{e}}^4 \mathrm{c}^3 / 3 \mathrm{h}^3\right) \mathrm{x}^3.\label{4.1.32}\]

Now neglecting inverse \(\beta\) decay the local density will be roughly given by \(\mathrm{\rho=m_p \rho_{\mathrm{e}} / m_{\mathrm{e}}}\) and

\[\mathrm{x^3=\left[\frac{3 h^3}{\left(8 \pi m_e^3 c^3 m_p\right)}\right] \rho}.\label{4.1.33}\]

Let \(\rho\) be given by its average value so that

\[\mathrm{\overline{x^2}=\left[\frac{9 h^3}{\left(32 \pi^2 m_e^3 c^3 m_p\right)}\right]^{\frac{2}{3}} \frac{M^{2 / 3}}{R^2}}.\label{4.1.34}\]

Normalizing R by Schwarzschild radius we get

\[\overline{\mathrm{x}^2}=7 \times 10^6\left(\mathrm{M}_{\odot} / \mathrm{M}\right)^{4 / 3}\left(\mathrm{R}_{\mathrm{S}} / \mathrm{R}_0\right)^2.\label{4.1.35}\]

This can be rigorously combined Equation \ref{4.1.31} to provide value for (Ro/Rs). However, since this entire argument is illustrative we also assume that the mass is roughly the limiting mass for white dwarfs. Then Equation \ref{4.1.31} becomes:

\[\left(\mathrm{R}_0 / \mathrm{R}_{\mathrm{S}}\right)>228\left(\mathrm{M}_{\odot} / \mathrm{M}\right)^{4 / 9} \cong 200,\label{4.1.36}\]

which is in remarkable agreement with the more precise figure of Chandrasekhar and Tooper² of 246. It is most likely that the discrepancy arises from the rather casual way of estimating \(\overline{\mathrm{x}^2}\) since it will be affected by both the type of volume averaging to determine \(\rho\) and the details of the Equation of state used in relating \(\rho\) to \(\rho_\mathrm{e}\). However, it should be remembered that the result is also only correct in the post Newtonian approximation and is an inequality setting a lower limit on instability. The interesting result is that General Relativity becomes important, indicating that instability sets in, at many times the Schwarzschild radius. This is the same result that Fowler found for supermassive stars supported by radiation pressure and serves as some justification for using the post-Newtonian approximation. It should be noted that substituting \(\mathrm{R_0}\) into the mass-radius relation for white dwarfs suggests that the critical mass should only be reduced by 1.5 percent. Hence the Chandrasekhar limiting mass for white dwarfs, while being somewhat too large, is still an excellent approximation.