4.8: Construction of a Model Stellar Interior

The construction of stellar models in steady state is essentially a numerical procedure which has been the subject of study of a large number of astrophysicists since the early 1950s and the pioneering work of Harm and Schwarzschild². Basically two methods have been employed to solve the equations. The early work utilized a scheme described by Schwarzschild which amounts to a straightforward numerical integration of the differential equations of stellar structure. In the early 1960s, this procedure was superceded by a method due to Henyey which replaces the differential equations with a set of finite difference equations whose solution is carried out globally and enables one to include time-dependent phenomena in a natural way. However, since this method requires an initial solution which is usually obtained by the Schwarzschild procedure, we describe both methods.

a. Boundary Conditions

Using the functional relations given by equations \ref{4.6.2}, we may reduce the problem of solving the structure equations to one of finding solutions for the four differential equations given in equations \ref{4.6.1}. These constitute a set of four nonlinear first-order differential equations in four unknowns. In general, such a system will have four constants of integration which must be specified to guarantee a solution. In principle, two of these constants are specified by requiring that the model be physically reasonable. These are \[\begin{aligned}
M(r) \rightarrow 0 & \text { as } r \rightarrow 0 \\
L(r) \rightarrow 0 & \text { as } r \rightarrow 0
\end{aligned}\label{4.7.1}\]

At the other end of the range of the independent variable, \[\begin{aligned}
M(r) & \rightarrow M_* & & \text { as } r \rightarrow R_* \\
L(r) & \rightarrow L_* & & \text { as } r \rightarrow R_*
\end{aligned}\label{4.7.2}\]

However, five constants are specified by equations \ref{4.7.1}, and \ref{4.7.2}, if R* is included as a parameter. Only four of these can be linearly independent. Thus, if one specifies M* and R*, the solution will specify L*. Another aspect of the problem is that the constants are not all specified at the same boundary, and so it is not possible to treat the problem as an initial-value problem and to solve by straightforward numerical integration. Such problems are known as two-point boundary-value problems, and one must essentially guess the missing integration constants at one boundary, obtain the numerical solution complete to the other boundary, and adjust the guesses until the specified integration constants at the far boundary are obtained. A further problem arises from the fact that the equations of hydrostatic equilibrium and energy transport are numerically unstable as r → 0 because the derivatives require the calculation of "0/0" at the origin. However, the problem can be recast as a double-eigenvalue problem with the fitting (solution adjustment) taking place in the interior but away from the boundary. This is essentially the Schwarzschild approach.

b. Schwarzschild Variables and Method

When one is searching for the numerical solution to a physical problem, it is convenient to re-express the problem in terms of a set of dimensionless variables whose range is known and conveniently limited. This is exactly what the Schwarzschild variables accomplish. Define the following set of dimensionless variables: \[\begin{array}{ll}
x=\frac{r}{R_*} & \\
q=\frac{M(r)}{M_*} & p=P\left(\frac{G M_*^2}{4 \pi R_*^4}\right)^{-1} \\
f=\frac{L(r)}{L_*} & t=T\left(\frac{m_h \mu G M_*}{k R_*}\right)^{-1}
\end{array}\label{4.7.3}\]

Note that the first three variables are the fractional radius, mass, and luminosity, respectively, while the two at the right represent the pressure and temperature normalized by a constant which describes the way they vary homologously. In addition, let us assume that the opacity and energy generation rate can be approximated by \[\bar{\kappa}=\kappa_0 \rho^n T^{-s} \quad \epsilon=\epsilon_0 \rho^\lambda T^v\label{4.7.4}\]

The differential equations of stellar structure then become \[\begin{array}{l}
\text { (a) }& \frac{d p}{d x}&=-\frac{q p}{t x^2} & \text {hydrostatic equilibrium }\\
\text { (b) }& \frac{d q}{d x}&=\frac{p x^2}{t}& \text {mass conservation }\\
\text { (c) }& \frac{d f}{d x}&=D(\lambda, v) x^2 p^{\lambda+1} t^{v-\lambda-1}&\text {energy conservation }\\
\text { (d) }&\frac{d t}{d x}&=\frac{-C(n, s) f p^{n+1}}{x^2 t^{n+s+4}}&\text {radiative equilibrium }\\
&\frac{d \ln p}{d \ln t}&=2.5&\text{convective equilibrium}
\end{array}\label{4.7.5}\]

which are subject to the boundary conditions \[q(0)=f(0)=0 \quad q(1)=f(1)=1\label{4.7.6}\]

The parameters \(\mathrm{C}(n,\mathrm{s})\) and \(\mathrm{D}(\lambda,v)\) are the eigenvalues of the problem, and these values specify the type of star being considered. In physical variables they are \[\begin{aligned}
C(n, s) & =\frac{3}{4(4 \pi)^{n+2} a c}\left(\frac{k}{m_h G}\right)^{s+4} \frac{\kappa_0}{\mu^{s+4}} \frac{L_* R_*^{s-3 n}}{M_*^{s-n+3}} \\
D(\lambda, v) & =\frac{1}{(4 \pi)^\lambda}\left(\frac{m_h G}{k}\right)^v\left(\epsilon_0 \mu^v\right) \frac{M_*^{v+\lambda+1}}{L_* R_*^{v+3 \lambda}}
\end{aligned}\label{4.7.7}\]

Note that the ideal-gas law has been used to eliminate the density from the problem, and this may cause some problems with the solution at the surface where the pressure and temperature essentially go to zero, in addition to the numerical problems at the center when x → 0. However, Schwarzschild shows that near the surface one may approximate the dimensionless pressure \(p\) and dimensionless temperature \(t\) by \[\begin{aligned}
& t \approx \frac{n+1}{n+7.5} \frac{1-x}{x} \\
& p \approx(\text { const }) t^{(n+7.5) /(n+1)}\quad\text { for radiative envelopes }
\end{aligned}\label{4.7.8}\]

If the star has a convective core, then all the energy is produced in a region where the structure is essentially specified by the adiabatic gradient and so the energy conservation equation [equation \ref{4.7.5c}] is redundant. This means that the eigenvalue \({\mathrm{D}(\lambda,v)}\) is unspecified and the problem will be solved by determining \({\mathrm{C}(n,\mathrm{s})}\) alone. Such a model is known as a Cowling model. The additional constraints on the solution are specified by the mass and size of the convective core (q_c and x_c). These are determined by the value of x for which \(\mathrm{d}(\ln p) / \mathrm{d}(\ln t)<2.5\), and the star becomes subject to the radiative temperature gradient. The stellar luminosity is then and for the envelope \(\mathrm{f=1}\). While such a scheme works well for models with convective cores, numerical problems will generally occur at the center should it be in radiative equilibrium and the solution obtained numerically. However, a slightly different set of dimensionless variables can be defined where the pressure and temperature are scaled by their values at the center of the star. The differential equations of stellar structure become stable at \(\mathrm{r=0}\) since the dimensionless pressure and temperature are both unity at the center by definition. One then, integrates outward from the center with P_c and T_c as eigenvalues. The stellar mass, luminosity, and radius can be related to these new eigenvalues. That there are two distinct eigenvalues is demonstrated by the surface boundary condition that both the surface pressure and surface temperature must vanish at the same value of x.

Unfortunately the equations of stellar structure become numerically unstable near the surface for the same reasons that required the approximation of the solutions of equations \ref{4.7.5} by equation \ref{4.7.8}. Although the errors in the model can be made small with the aid of modern computers, it is bad practice to numerically solve equations which are inherently unstable. For that reason, the usual procedure is to integrate from both the outside and the inside and to make the fit at the boundary between the core and envelope. The approximations near the surface are still present, but their effect on the solution is minimized. In actual practice, the fitting can be accomplished in the U-V plane where the solutions are homologously invariant. The fitting procedure is similar to that described in Chapter 2.

Since Schwarzschild introduced this method of solution of the equations of stellar structure in the 1950s, many variants have been used by numerous investigators. In one form or another, all variants suffer from problems similar to those that plague the Schwarzschild procedure. In general, this approach to the numerical solution of two-point nonlinear boundary-value problems always suffers from the propagation of errors from one boundary to the other. The most serious of these errors are usually the truncation errors associated with the numerical integration scheme which tend to be systematic. However, this approach enabled the generation of stellar models which represented the steady-state aspect of stars for the first time. Although qualitative information about stellar evolution can be gained from polytropes (and we do so in Chapter 5), specific and detailed descriptions of stellar evolution require the generation of steady-state models. However, some aspects of stellar evolution happen on time scales which are very short compared to the thermal time scale, and in some instances short compared to the dynamical time scale. Often, substantive changes occur to the internal structure which produces only small changes at the surface. Thus, minor changes in the surface boundary conditions can reflect monumental changes in the internal structure of the star. In addition, we must include the time-dependent terms in the equations describing the conservation of momentum and energy. Specifically, if some of the generated energy does work on the star, causing it to expand as energy is liberated by contraction, then this energy must be included in the energy conservation equation relating the stellar luminosity to the sources of energy. This is usually accomplished by keeping track of the time rate of change of the entropy. The direct integration scheme does not readily lend itself to the inclusion of such terms. Such models are no longer merely steady-state models, and we will require more sophisticated tools to deal with them.

c. Henyey Relaxation Method for Construction of Stellar Models

To overcome some of the numerical instability problems described in the previous section, Louis Henyey et al³. developed a superior numerical scheme in the early 1960s. This method is the foundation for all modern stellar-model calculations. His approach was to transform the problem to a set of variables in which the nonlinearity of the differential equations was minimized. The differential equations of stellar structure were then replaced with a set of finite difference equations whose solution could be carried out simultaneously over the entire model. This tended to reduce the effect of truncation error by spreading it more or less evenly across the model. Furthermore, the addition of time-dependent terms proved to be relatively easy to incorporate in the structure equations. We do not describe all the details that make this method so powerful, but only sketch the principles involved.

We begin by replacing the independent variable r with M(r). Henyey noticed that the behavior of the equations was far more linear when the mass interior to r was used as the independent variable. The radial coordinate then becomes a dependent variable whose value must be found for any particular M(r). If we make this transformation, the four differential equations of stellar structure become \[\begin{array}{l}
\frac{\partial P}{\partial M(r)}=\frac{G M(r)}{4 \pi r^4} & \text { hydrostatic equilibrium } \\
\frac{\partial r}{\partial M(r)}=\left(4 \pi r^2 \rho\right)^{-1} & \text { conservation of mass } \\
\frac{\partial L}{\partial M(r)}=\epsilon(P, T, \mu)-T \frac{\partial S}{\partial t} & \text { conservation of energy } \\
\frac{\partial \ln T}{\partial \ln P}=f(P, T, \mu) & \text { energy transport }
\end{array}\label{4.7.9}\]

Here we have explicitly included the time dependent entropy term in the energy equation for purposes of example. In addition, we have written the energy transport term in a general manner which can accommodate either radiation or convection. Now we divide the star into N - 1 zones, starting with the center as the first point and ending at the surface or some outer point where the boundary conditions are known. By approximating the derivatives of equations \ref{4.7.9} by the difference of the parameters at adjacent points, we get the following finite difference equations: \[\begin{aligned}
& \mathrm{\frac{P_{i+1}-P_i}{M_{i+1}-M_i}=\frac{G M_{i+1 / 2}}{4 \pi r_{i+1 / 2}^4} }\\
& \mathrm{\frac{r_{i+1}-r_i}{M_{i+1}-M_i}=\left(4 \pi r_{i+1 / 2}^2\right)^{-1} P_{i+1 / 2}^{-1}} \\
& \mathrm{\frac{L_{i+1}-L_i}{M_{i+1}-M_i}=\varepsilon_{i+1 / 2}-\left.T_{i+1 / 2} \frac{\partial S}{\partial t}\right|_{i+1 / 2}} \\
& \mathrm{\frac{T_{i+1}-T_i}{P_{i+1}-P_i}=\frac{T_{i+1 / 2}}{P_{i+1 / 2}} f\left(P_{i+1 / 2}, T_{i+1 / 2}, M_{i+1 / 2}\right)}
\end{aligned}\label{4.7.10}\]

The subscript i + ½ is used exclusively on the right-hand side of equations \ref{4.7.10} to indicate that the value to be used is intermediate between the values at i and i + 1. It will turn out that we must have an initial guess of the model's structure in order to solve the finite difference equations. It is this guess which may supply the initial information for evaluating the parameters at the points i + ½. Since the mass points M_i represent the independent variable of this problem, the four equations given in equations \ref{4.7.10} contain eight unknowns. However, we have N - 1 systems of such equations with considerable overlap in unknowns among them. The situation at the outer zone will be handled somewhat differently since there is no N + 1 point. Thus if we count the total number of equations we have 4N - 4. But at each point there are only four unknowns, making the total number of unknowns of 4N. The remaining four constraints are essentially the boundary conditions of the problem. By analogy to the Schwarzschild problem, let us take the central boundary conditions to be r₁ = L₁ = 0, which removes two of the additional unknowns. Now if we choose two of the remaining unknowns at the surface, such as rN and LN, the problem is completely specified. Indeed, if we choose the surface pressure to be zero, then choosing a star of a particular mass and radius (and distribution of chemical composition) will specify the stellar configuration. One of the motivating notions that led Henyey to this type of technique was the ability to match a stellar interior to a model of the stellar atmosphere. This technique is ideally suited to do this. One simply takes as the outer zone that point where the physical parameters are known as the result of a separate study. In the second part of this book, we present a theory of stellar atmospheres which provides far more accurate surface boundary conditions than those of early investigators. In addition to improving the manner in which the surface boundary conditions are handled, it may be advisable to ignore the point at the center. A Taylor series expansion can be used to express the values of P₂, T₂, L₂, and r₂ in terms of the central temperature T_c and P_c. Because the system of equations is strongly diagonal, the solution is easier to come by if the central boundary conditions are expressed in this manner.

The Henyey approach shown in equations \ref{4.7.10} represents the derivative of the structure equations by first-order finite differences. Thus the errors of the approximation are second order in those derivatives. This necessitates the use of the large number of zones to accurately represent the model, and it is this large number of zones that represents the primary computational burden in the construction of the model. Although increasing the order of the finite difference equations would improve the stability, it would also increase the density (i.e., the number of terms) of the resulting linear algebraic equations, slowing their solution and decreasing their stability. Budge⁴ has shown that an improvement in the accuracy of the approximation can be achieved by using a Runge-Kutta fourth-order approximation for the derivatives without increasing the resulting linear equation density. Although there is some increase in the computational burden for obtaining the coefficients, this is more than offset by being able to reduce the number of zones in the model.

We still must solve these linear equations. It is not uncommon in the standard Henyey scheme to choose up to 500 zones in the star, which will yield some 2000 nonlinear equations in as many unknowns. Now it is clear why we need an initial solution. If we have a solution which is close to the correct one, we may express the correct solution in terms of the initial solution and a small linear correction to that solution. This will reduce the system of nonlinear equations to a linear system where the corrections are the unknowns. Such a scheme is known as a Newton-Raphson iteration scheme. Since the system is sparse (each equation contains only 8 of the 2000 unknowns) and the independent variable was chosen so as to make the equations somewhat linear, the iteration scheme is usually stable. However, the stability also depends on the quality of the initial solution. This is normally obtained by means of a Schwarzschild-type integration or a previously determined model.

It is clear that the Henyey method lends itself naturally to the problem of stellar evolution. In this case the initial model is a model calculated for an earlier time. Thus the procedure would be to start at some initial time with a Schwarzschild model, allow a small interval of time for the model to pass by, calculate the changes in the chemical composition resulting from nuclear processes, and modify the model accordingly. This serves as the initial first guess for the Henyey scheme, and a new model is produced. The effects of time are again allowed for, and the next Henyey model is constructed, etc. In this way an entire sequence of stellar models representing the life history of the star can be constructed. One generally starts the sequence when the star is well represented by a steady-state model, and the Schwarzschild solution gives an accurate description of the stellar structure. Such a time is the arrival on the main sequence and the accompanying onset of hydrogen burning. The resulting life history of the star is as good as the microphysics which goes into the solution and the quality of the computer and the associated numerical techniques used to obtain the solution.

At this point, we have covered the fundamentals required to construct a model of the interior of a star. However, we should not leave the impression that such a model would reflect the accuracy of contemporary stellar interior models. There are many complications and refinements which should be treated and included to produce a model with modern sophistication. We have said nothing about the small departures of the equation of state from the ideal-gas law which occur at quite modest densities due to electron screening. Nor have we dealt with many of the vagaries of the theory of convection, such as semi-convection, convective overshoot, or mixing-length determination. These result largely from the primitive nature of the existing theory of convection, and while they do pose significant problems at certain points in a star's evolution, they do not affect the conceptual picture of stellar structure. It seems almost criminal not to devote more attention to the efforts of those who have labored to provide improved opacities and nuclear energy generation rates. But again, while these improve the details of the models and enhance our confidence in the predictions based on them, they do not conceptually change the basic physics upon which the models rest. While we have outlined the numerical procedures necessary to actually solve the structure equations, there is much cleverness and imaginative numerical analysis required to translate what we outlined to a computer program which will execute to completion in an acceptable time. Do not forget that the early models of Schwarzschild and Harm were calculated basically by hand, aided only by a desk calculator whose capabilities are far exceeded by even the cheapest pocket calculator of the present. It is no accident that the rapid advance of our knowledge of stellar structure parallels the explosive advance in our ability to carry out numerical calculations.

An understanding of the refinements of contemporary models is essential for any who would choose to do meaningful research in stellar interiors. It is not essential for those who would understand the results and their physical motivation, and it is to those people that this book is addressed. With the knowledge of the physical processes that determine the structure of stars, let us now turn to the crowning achievement of the study of stellar structure - the theory of stellar evolution.