Genetic-Algorithm-based Optimization of White Dwarf
Pulsation Models using an Intel/Linux Metacomputer

White dwarf asteroseismology offers the opportunity to probe the structure and composition of stellar objects governed by relatively simple principles. The observational requirements of asteroseismology have been addressed by the development of the Whole Earth Telescope (WET), but the analytical procedures need to be refined before this technique can yield the complete physical insight that the data can provide. I propose to apply an optimization method utilizing a genetic algorithm (GA) for fitting white dwarf pulsation models to WET data. I will use this global approach to optimization in two distinct ways: (1) I will investigate the uniqueness and objectivity of our solutions by combining existing models with a GA-based fitting routine, and (2) I will investigate the completeness and adequacy of our understanding of the principles governing white dwarf interiors by parameterizing the constitutive physics and using a GA to find the family of solutions that produce observationally indistinguishable behavior. I will configure a specialized computational instrument to perform the calculations -- a metacomputer (consisting of a network of minimal PCs running Linux) capable of ~4 GigaFLOPS (billion floating-point operations per second). The metacomputer is already operational, needing only the addition of identical nodes to those present to achieve the needed performance, and the funding for this addition is already available.

We study white dwarf stars because they are the end-points of stellar evolution for the majority of all stars, and their composition and structure can tell us about their prior history. We can determine the internal structure of pulsating white dwarfs by observing their variations in brightness as a function of time, using the techniques of high speed photometry to define their light curves, and then matching these observations with a computer model which behaves the same way. The parameters of the models are chosen to correspond one-to-one with the physical processes that give rise to the variations, so a good fit to the data leads us to believe that our model reflects the actual physics of the stars themselves.

Although this procedure is simple in outline, its realization in practice requires specialized instrumentation to overcome the practical difficulties we encounter in the process. The Whole Earth Telescope (WET) observing network was developed to provide the 200 or more hours of essentially gap-free data we require for the analysis. This instrument is now mature, and has provided a wealth of seismological data on the different varieties of pulsating white dwarf stars, so the observational part of the procedure is well in hand. Now we need to improve our analytical procedures to take full advantage of the possibilities afforded by asteroseismology.

In order to understand the oscillations that we actually observe, we need to consider non-radial spheroidal oscillations in general. We arrive at the form of the equations that we actually solve by perturbing the fluid equations, and keeping only the terms of lowest order. We assume a static, spherical equilibrium structure which is given by a theoretical evolutionary model. Because the surface gravity is high (log g ~ 8 in cgs units) and rotation rates are typically of order days, the spherical approximation is very good, and we can expand our solutions in terms of spherical harmonics, (Y_lm). The most relevant parameters which describe the equilibrium model are the Brunt-Vaisala frequency (which is just the difference between the actual and the adiabatic density gradients), and the Lamb or acoustic frequency (which is related to the local sound speed). We can use plots of the run of these two parameters through the star to determine where modes of given frequencies will propagate. Such diagnostic plots are known as propagation diagrams, and contain essentially all of the information we can obtain through asteroseismological analysis of the star. In a sense, the ultimate goal is to use the distribution of observed frequencies to empirically determine the propagation diagram for the star. Once this is known, all of the constitutive physics can be deconvolved -- at least in principle. This kind of analysis is called seismological inversion, and has been used successfully so far on only one star, the Sun.

We can gain a great deal of physical insight into the oscillations, and how they sample the star through a sort of local analysis (cf. Unno et al. 1989). If we assume a radial dependence of the form e^ik_rr, and wavelengths short compared to the relevant scale-heights for the physical quantities, then we arrive at a local dispersion relation (LDR). In order for a given mode to be locally propagating k_r² must be positive, and the form of the LDR reveals that this can only occur when the oscillation frequency is either greater than both the acoustic and Brunt-Vaisala frequencies, or is less than both. In the limits of large and small frequencies, the LDR yields two physically distinct kinds of solutions that represent the two principal classes of non-radial spheroidal modes. The first limit yields the class of solutions representing the p-modes (so called because pressure is the principal restoring force). Radial displacements are dominant, and for white dwarf stars these have timescales of seconds -- too short for the periods of the observed oscillations. In any case, one would not expect to observe large radial displacements on such a high-gravity object. The second class of solutions represent the g-modes, where gravity is the dominant restoring force. These have timescales of hundreds of seconds and longer -- just like the observed oscillations in the white dwarf stars.

The adjustable parameters in our computer models of white dwarfs presently include: the total mass, the H and He layer masses, composition, and transition zone thicknesses. Finding a proper set of these to provide a close fit to the observed data is difficult. The current procedure is a cut-and-try process guided by intuition and experience, and is far more subjective than we would like. More objective procedures are essential if asteroseismology is to become a widely-accepted astronomical technique. We must be able to demonstrate that, within the range of different values the model parameters can assume, we have found the only solution, or the best one if more than one is possible. We plan to apply a search-and-fit technique employing genetic algorithms, which can explore all of the myriad parameter combinations possible and select for us the best one (or ones).

An optimization scheme based on a genetic algorithm (GA) can avoid the problems inherent in more traditional approaches. Restrictions on the range of the parameter space are imposed only by observations and by the physics of the model. Although the parameter space so defined is often quite large, the GA provides an efficient means of consistently finding the model which results in the absolute minimum variance when compared to the observational data. While it is difficult for GAs to find precise values for the set of `best fit' parameters, they are very good at finding the region of parameter space that contains the global minimum. In this sense, the GA is an objective means of finding a good first guess for a more traditional method which can narrow in on the precise values and uncertainties of the `best fit'.

The underlying ideas for genetic algorithms were inspired by Darwin's (1859) notion of biological evolution through natural selection. The first comprehensive description of how to incorporate these ideas in a computational setting was written by Goldberg (1989). In the first chapter of his book, Goldberg describes the implementation of a simple GA. Initially, trial solutions are generated at random throughout the parameter space. Each solution is assigned a `fitness' inversely proportional to the variance, and a new set of solutions is constructed from the old one, weighted by the fitness. In order to explore new regions of the parameter space, changes are introduced by `breeding' and `mutating' the solutions. The fitness is evaluated in the new population, and the process continues until the majority of solutions are in one region of parameter space. In practice, this technique converges in a fraction of the time that would be required by other comparably global methods, such as a simple grid-search of the parameter space.

Genetic algorithms have been used a great deal in other fields, but until recently they have not attracted much attention in astronomy. The application of a GA to problems of astronomical interest was first promoted by Charbonneau (1995), who demonstrated the technique by fitting: the rotation curves of galaxies, a multiply periodic signal, and a magnetohydrodynamic wind model. Since then, several other applications of GAs to astronomical problems have appeared in the literature. Tomczyk et al. (1995) used a GA with helioseismological models to constrain solar core rotation. Hakala (1995) optimized the accretion stream map of an eclipsing polar. Lang (1995) developed an optimum set of image selection criteria for use in the detection of high-energy gamma rays. Lazio (1997) searched pulsar timing signals for the signatures of planetary companions. Most recently, Wahde (1998) used a GA to determine the orbital parameters of interacting galaxies. The applicability of GAs to such a wide range of astronomical problems is a testament to their versatility.

Although genetic algorithms are more efficient than other comparably global techniques, they are still quite demanding computationally. To be practical, the code that I am proposing to develop will require a dedicated instrument to perform the calculations. Over the past six months, we have designed such an instrument -- a collection of about 40 minimal PCs connected by an isolated network. Since the structure of a GA is very conducive to parallelization, this metacomputer will allow us to run our code much faster than would otherwise be possible. I have configured a fully operational 12-node version of the metacomputer, and I now need only acquire and assemble the rest of the hardware with existing funds.

Once the hardware is complete, I will need to develop the software for a generalized parallel genetic algorithm. To do so, I will simply integrate two previously developed programs which are both in the public domain. (1) PIKAIA is a general purpose optimization subroutine based on a genetic algorithm. Given an arbitrarily complicated function, this routine will seek the global maximum within the defined bounds of the parameter space. (2) PVM allows a collection of networked computers to cooperate on a problem as if they were a single multi-processor parallel machine. The software consists of a daemon which runs on each host in the virtual machine, and a library of routines that need to be incorporated into the code so that it can utilize the available resources. When the elements of these two programs are properly merged, the result will be a generalized parallel genetic algorithm which I will use for both of my approaches to optimization of white dwarf pulsation models.

My first application of the parallel genetic algorithm will be to establish an objective fitting procedure for existing white dwarf pulsation models (the forward problem). Since the GA begins with a uniform random sampling of parameter space, it will have no predisposition towards any particular combination of parameters except for those which it objectively finds to be favorable. As the evolution of the final solution progresses, the GA will reveal the structure of the parameter space -- local minima, degeneracies, and correlations. In the end, the GA will provide an informed first guess for the final non-linear least squares fit. This procedure will either ensure that our fit is unique, or it will reveal a family of solutions which are all sufficient to describe the observations.

If our representation of the physics of white dwarf interiors is incomplete or inadequate, and especially if there is new physics to be found in these extreme conditions, we will not find it with the forward approach. To address these important possibilities directly, I will parameterize the physics into two simple functions representing the acoustic and the Brunt-Vaisala frequencies at various depths (the inverse problem). With the appropriate approximations, these two functions are sufficient to determine the set of allowed pulsational modes, of which some subset would be observed. I will use the GA to evolve the family of functions which produce identical sets of observed modes. This will probe the completeness and adequacy of our understanding by revealing which specific sections of our description of the interior (e.g. the functions describing the acoustic and Brunt-Vaisala frequencies) are constrained by observations, and which are not. By using pulsating white dwarf stars in this way -- as a laboratory for studying matter under conditions not accessible in any other way -- I will either reveal some new physics, or I will demonstrate for the first time the remarkable accuracy of our present understanding.

The new approach to the forward problem will immediately allow us to make objective mode identifications for the archive of asteroseismological data on nearly two dozen oscillating white dwarf stars. This will provide an independent check of the identifications that were made previously by the traditional, more subjective procedure. Because each mode offers a different constraint on the interior, results from the new approach to the inverse problem will come most readily from those white dwarfs with the greatest number of observed modes. Although PG 1159 presently claims this distinction, it is classified as a DOV, so the internal physics are complicated by its relatively rapid evolution (Winget et al. 1991). As a consequence, the best candidate is the DBV white dwarf GD 358 which has more than ten modes of differing radial overtone number, as well as a plethora of modes that have been identified as combination frequencies. After finding constraints from the best, we will simply work our way down the list of white dwarfs with many observed modes (e.g. PG 2131, PG 0122) and those with modes probing regions of the interior which remain unconstrained by earlier analyses.

The primary concern for this project is computing time. We have designed the metacomputer with the specific demands of this project in mind. It should provide the resources and flexibility that we need for developing and testing our code. If the final problem requires more computing power than we have allowed for, we can either scale up the metacomputer (with additional funding), or run the code on local supercomputing facilities with a few minor modifications.

Our investigation of the inverse problem promises to specify those regions of white dwarf interiors which are unconstrained by observations. We may find that observations of certain pulsational modes would provide the necessary constraints, and this could motivate WET observations of a specific target. As we begin to learn more about how to map pulsational modes directly to the internal physics of a specific region of the interior, our choice of targets for future WET runs can be guided appropriately. My participation in these campaigns will require continued allocation of telescope time to the WET, and possibly funding for travel to participating observatories.