We attempt to develop a code to describe the structure and evolution of a star.
The most basic models of stars are defined by the Russell-Vogt theorem, which states that a star's mass and chemical composition are sufficient to determine its structure and evolution.
The determinations of structure and evolution can be done using a set of tools, many of which are derivable from basic physical assumptions:
1. Stellar structure differential equations
2. Description of stellar material under conditions of particular composition, density, and temperature
3. Boundary conditions for structure equations, at stellar center and stellar surface
We will describe each of these in turn in order to produce a general picture of the physical problem. For our purposes we will assume a spherically-symmetric main-sequence star in complete equilibrium (i.e., no structural hydrostatic or thermal time dependence). We will also work within the limitations of a time-independent chemical composition. In essence, we will calculate a snapshot in time of a main-sequence star.
Gas constant
Rg = 8.3144 ⋅ 107 erg g-1 K-1
Radiation constant
a = 7.5657 ⋅ 10-15 erg cm-3 K-4
Gravitational constant
G = 6.672 ⋅ 10-8 cm3 s-2 g-1
Speed of light
c = 2.99792458 ⋅ 1010 cm s-1
Stefan-Boltzmann constant of radiation
The four time-independent Lagrangian stellar structure equations for a spherically-symmetric star describe changes in radius, luminosity, pressure, and temperature:
where
ρ is density
εn is nuclear energy generation rate
εν is energy loss to neutrino production
m is m(r), mass internal to radius r
∇ is temperature gradient, 
dr/dm is derived from the fact that 
.
For dl/dm, l(m), or luminosity at a mass point m within the star, is given by the total energy generation occurring internal to mass point m. Main-sequence stars produce energy via fusion of light elements. εn gives the nuclear energy generation rate per mass per time.
dP/dm is derived from the requirement of hydrostatic equilibrium. Stars support themselves against gravitational collapse by a pressure gradient; that is, for a mass element dm,
This translates to the requirement
,
.
Energy transport in stars occurs mainly by radiation and convection. dT/dm takes on the form
=
=
,
that is, ∇ describes variation in temperature with depth, where pressure is a proxy descriptor for "depth," since pressure grows with decreasing stellar radius (and internal mass). Although ∇ should be calculated directly by solving the equations of mixing length theory, such a technique is beyond the scope of our stellar code, and we will use only two forms for ∇, one adiabatic and the other radiative. Where the efficiency of convection is very high, adiabatic ∇ is a good approximation for actual ∇. Where convection is ineffective, energy transport occurs largely by radiation and we use radiative ∇. The test for dynamical stability is given by the Schwarzschild criterion:
∇ad < ∇rad implies stability such that convection will not occur
Given the proper boundary conditions at the center and the surface and initial guesses for those parameters which we cannot easily bound (T and P at the center, r and l at the surface), we should be able to integrate these ordinary differential equations both inward and outward from the surface and center of the star respectively, to check the end values of the integrated functions at some appropriate mass point between the center and the surface, and to adjust our initial guesses and reintegrate over as many iterations as necessary until the results of our inward and our outward integrations converge.
Our stellar structure equations contain several parameters whose values vary with composition, density, and temperature. We are assuming static composition, so will concern ourselves with variations only in density and temperature.
First, let us define our composition. Stellar composition is given in terms of mass fraction Xi:
, for element i
Because most elements are rare in stars, most Xi need not be specified. In a basic stellar model it is enough to give:
Specifying at least XC, XN, and XO would further refine the model because of their relevance to energy generation, but we will not concern ourselves with those here. Z is typically low, usually ranging from 0.04 to 0.001. X generally dominates Y by roughly a factor of 2.
The mean molecular weight μ of a free particle in the star is given in units of the mass of the hydrogen atom and in a simple star, wherein all atoms are fully ionized due to high temperatures,
Next we define a pressure equation of state. In our basic model, we will assume that our star is composed of an ideal gas, fully ionized, with possible pressure effects from radiation:
All of our calculations of density, ρ, will be directly derived from our pressure equation of state:
The opacity, κ, of stellar material depends on composition, density, and temperature. Calculating opacity is computationally intensive, usually done numerically for set values of ρ, T, X, Y, and Z. We will interpolate in ρ and T in a table of κ values done for set X, Y, and Z in the OPAL Opacity Code project (http://www-pat.llnl.gov/Research/OPAL/). Specifically, we use tables calculated by Grevesse and Noels (1993) for solar composition.
We need descriptions of ∇ for the adiabatic and radiative cases. First, we calculate the fraction of the total pressure that is given by ideal gas pressure (as opposed to radiation pressure):
Adiabatic ∇ can be calculated through manipulation of β (we will not derive this formula, as it requires extensive discussion):
Radiative ∇ is derived from radiative transfer equations which depend on a diffusive interpretation of photon travel (again, we will not discuss the derivation here):
Now we should investigate the nuclear energy generation rate, εn. The lion's share of energy generation comes from hydrogen burning, which produces helium nuclei. Although the contributions of helium burning and the burning of heavier elements (carbon, oxygen, etc.) are not insignificant, they require temperatures higher than those found in the cores of most main-sequence stars, so we will neglect them.
For hydrogen burning we must calculate the contributions of the proton-proton chain and the CNO cycle. The proton-proton chain begins by fusing two hydrogen nuclei and eventually results in the production of helium nuclei. The CNO cycle uses certain isotopes of carbon, nitrogen, and oxygen as catalysts in the production of helium nuclei from hydrogen nuclei. The contributions of each are as follows:
Each is generally determined by the rate of its slowest constituent reaction. The gjk are correction factors to the Gaussian of the Gamow peak, which describes the warring effects of Maxwellian distribution and quantum-tunneling probability in determining which particles have sufficiently high energies to overcome the Coulomb barrier and fuse. fjk is a screening factor that increases the eligible "reaction volume" of a particle. ψ is a correction factor to the reaction rate, dimensionless, which changes with T and differs for each branch of the proton-proton chain. It usually falls between 1 and 2 and we will treat it as a constant value in our model. T6 is merely T/106.
Therefore
in our model.
The simplest boundary conditions for our problem are
center: m = 0 ⇒ r = 0, l = 0
surface: T = 0, P = 0
We will refine these slightly.
At the center, using r = 0 and l = 0 in our stellar structure equations causes them to blow up to infinity. Therefore, we perform a series expansion of the structure equations about m = 0 in order to learn the behavior of our equations near the center:
dT4 is used in the radiative case, d ln T in the convective case.
The subscript c indicates the value at stellar center.
Because stars do not have precisely delineated "surfaces," taking T = 0 and P = 0 is not strictly true until a distance from the center of the star at which the stellar structure equations can no longer be applied with any accuracy in the vicinity of that radius. Therefore, we use the photosphere as the star's "surface." To find P, we set the optical depth equal to 2/3 (the definition of the photosphere) and integrate dP/dr from the stellar total radius R to infinity, with constant surface g and substituting the integral over ρ with that from the definition of optical depth:
To find T, we use the definition of the effective temperature at the photosphere:
Physical constants are provided at top.
Given the proper boundary conditions at the center and the surface and initial guesses for those parameters which we cannot easily bound (T and P at the center, r and l at the surface), we should be able to:
The mass fitting point should ideally be chosen to be at the interface between, for example, a convective core and a radiative envelope (K&W). However, for purposes of simplicity we have chosen a constant fitting point mF.
The adjustment of the initial guesses is done via simultaneous solution of four equations characterizing the differences between the original results of the inward and outward integrations, akin to a multidimensional Newton-Raphson method.
Three sets of trial integrations are done to determine the changes in the inward and the outward results with changes in the initial guesses; one integration for a preliminary result, one integration for change in inward integration with initial ltotal = L and change in outward integration with initial Tc, and one integration for change in inward with initial rtotal = R and change in outward with initial Pc.
The appropriate adjustments to the initial guesses are given by the following matrix solution of the four equations for the four unknown adjustments:
The adjustments are added numerically to the initial guesses. However, the correction assumes linear behavior of the changes in final values with changes in initial guesses; because the actual behavior may not be linear, we take only a fraction of the calculated correction in our actual adjustment.
Following the correction of the initial guesses, the integration proceeds again.
The stellar code is written in Python.
The function odeint from the module scipy (ScientificPython; http://www.scipy.org) is used to integrate the stellar structure equations. SciPy's odeint is a wrapper around the adaptive step-size integrator LSODA from the Fortran library ODEPACK, written at Lawrence Livermore National Laboratories (http://www.llnl.gov/CASC/odepack/).
The n-dimensional array object and some basic linear algebra functions (e.g., numpy.linalg.solve) from NumericalPython (NumPy; http://numpy.scipy.org) are used to solve the matrix formulation of the simultaneous solution problem for the initial guess corrections. NumPy requires full and preferably optimized versions of libraries for LAPACK, the linear algebra package written in Fortran77, and BLAS, the basic linear algebra subroutines package (see http://www.netlib.org/lapack and http://www.netlib.org/blas). SciPy requires NumPy.
The model used in the stellar code represents a simple solar-type star with M, R, L, and composition similar to those of the Sun. Because the integration has not yet converged, we will not yet attempt any qualitative or quantitative comparison of results with solar models found in the literature.