If some words seem slightly out-of-context, please keep in mind that this was a proposal for computer time, not just a summary of research. Also keep in mind that this text was 'converted' from a document that was originally in TeX format... so there are undoubtedly some apparent typos.
One of the major aspects of this proposal is the computation of the spectrum predicted to emerge from the complex composition, density, and temperature structure of an exploding model star. In this case, the extra dimension is frequency, rather than a spatial dimension. The result is a theoretical spectrum to compare to observed spectra of supernovae, and in so doing to map out the nature of the ejected material.
Other problems involve the computation of the explosion process in at least two, and preferably three dimensions. One class of supernovae involves the thermonuclear burning of white dwarf stars. The mechanism of ignition, propagation, and quenching of burning fronts is intrinsically multi-dimensional and we have developed the codes and techniques to do this problem properly in two and three dimensions for the first time. The other basic physical mechanism involved in supernova explosions is gravitational collapse of an evolved stellar core to produce a neutron star with the liberation of a huge flux of neutrinos. The suspicion has grown in recent years that convective flows in the newly formed neutron star are critical to transport the neutrinos outward to power the explosion. We have developed the techniques to handle both the complex fluid flow and the energy-dependent flux of neutrinos in two and three dimensions.
We thus propose a comprehensive, coordinated series of large scale computations to forge new ground in our understanding of the nature of supernovae. We will compute state-of-the-art multi-dimensional models of both thermonuclear explosions and core collapse and compute the light curves and spectra of the thermonuclear models to compare to observations. Some of these computations will strain the capacity of the C90, and we will actively explore the use of the T3D. The codes we have used are highly vectorized and in some cases ready to run in multi-processor or massively parallel mode.
There are still great uncertainties about the precise physical mechanism of the explosion, subsonic deflagration or supersonic detonation (see below). One dimensional detonation models tend to burn too much of the star to iron-peak matter, in violation of the observations and model spectra that show that the outer parts of the star remain only partially burned. Ad hoc one dimensional deflagration models have proved fairly successful in accounting for the spectral evolution, but the physics of the combustion front is not well understood. Other models have proposed a hybrid with transitions between deflagrations and detonations (Khokhlov 1991). Previous work (Harkness 1991) has shown that many models that fit the light curve (the time profile of the broad band light output) fail to agree with the observed spectral evolution. The critical diagnostic is the ability to compute the time-dependent spectrum of a given model to compare with the observations. The solution of the radiative transfer equations in this context is a major challenge since the ejecta have a complex composition profile and are strongly dominated by electron scattering.
The code SuperNova Atmosphere Program (SNAP) designed to address this problem has been under development by R. P. Harkness for well over a decade. SNAP will be used to compute the spectral evolution of (angle averaged) dynamical models produced in this program as well as other relevant models computed by us and others. We expect the results to sharply constrain the range of acceptable models and hence to be a key guide to the physical mechanism of the explosion. An automatic by-product of these calculations will be the intrinsic luminosity of the model. This is a key quantity in the sharpening debate on the distance scale of the Universe.
The quality and quantity of observational data on Type Ia supernovae has improved greatly in recent years. Of particular interest is the Type Ia SN 1992A. This event was well observed from the ground (Suntzeff 1993), but also was the first supernova observed by the Hubble Space Telescope in the ultraviolet and optical (Kirshner et al.1993). The HST observations were obtained by a consortium of which the PI is a co-investigator. The data on SN 1992A reveal the ultraviolet spectrum with unprecedented accuracy. As first investigated by Harkness (Wheeler et al.1986; Harkness 1991), the UV spectra constrain the density and velocity profiles of the outermost matter that are most sensitive to the environment in which the supernova explodes and to line-blocking effects in the atmosphere.
SN 1992A also provided important new data on a long-observed, but unexplained feature of Type Ia supernovae. They display a secondary maximum in the infrared bands that is not seen in the blue bands and in other classes of supernovae. Careful observations of SN 1992A show that the secondary peak shows faintly even in the optical V band, so it is quite broadband (Suntzeff 1993). In addition, the I band behaves rather differently than those at shorter wavelength, with its principle maximum earlier rather than later which is the general trend with wavelength. There is a recent suggestion (Spyromilio, Pinto, and Eastman 1993) that some of this modulation is related to the paucity of infrared Fe II lines in the atmosphere, but it is not clear that this can account for the broad band nature of the effect nor the particular transient time dependence observed. Light curves in the BVR and I bands are again an automatic product of the computed spectral evolution. Study of this secondary IR peak and the special behavior of the I band will complement the spectroscopic work to give deeper understanding of the nature of these explosions.
Type Ia supernovae are currently at the heart of a long-standing debate about the size and age of the Universe (Fukugita, Hogan and Peebles, 1993). Because of their great intrinsic brightness and apparent uniformity, Type Ia supernovae have been regarded as "standard candles" and used to measure distances to host galaxies. Current thermonuclear models that fit the spectra demand a large luminosity, large distances, and a small value of the Hubble constant (Arnett, Branch and Wheeler 1985). Many independent and apparently self-consistent techniques are yielding smaller distances to galaxies and a larger Hubble constant (de Vaucouleurs 1993). This is a fundamental disagreement. Either the apparently successful supernova models are wrong, or the alternative observational techniques are wrong. One of the issues that needs to be carefully rechecked is the relation between the total power deposited by radioactive decay in the explosion and the power that is radiated in limited observational band passes. We will explore the range in total luminosity and the luminosity in various observational band passes that can be compatible with the observed spectral evolution and hence put the use of Type Ia supernovae as distance indicators on a new quantitative basis.
Although many Type Ia supernovae seem to follow nearly identical spectral behavior, there have been significant deviations discovered recently. Among these are SN 1991T (Filippenko et al.1992a) that showed an excess of Fe in its early spectra, and SN 1991bg (Filippenko et al.1992b; Leibundgut et al.1993) that was distinctly subluminous. In addition to the intrinsic interest in determining the physical cause of these variations, this issue is also of significance for the determination of the extra-galactic distance scale. It is critical to determine the nature of the deviations and the intrinsic luminosity of all classes of models. We will explore the variations in burning mechanisms that can account for events like SN 1991T and SN 1991bg.
The physics of the detonation instability is the focus of of the PhD research of graduate student J. R. Boisseau. Boisseau has constructed a multi-dimensional hydrodynamic code (Explosion Properties Investigation Code: EPIC) to compute the detailed nature of the detonation fronts in exploding white dwarf stars. To do so, he chose to use the Flux-Corrected Transport method and has worked closely with Dr. Elaine Oran of the Naval Research Laboratory, an expert on both FCT and the nature of detonations and flame propagation. He has produced a highly efficient, vectorized code designed to run with maximum effect on a Cray. This code is now nearing final development stages and production runs are required to finish the project.
The collapse problem has been under intensive study for over two decades, but the theoretical picture remains highly uncertain. Recently the understanding has grown that a very important element of this picture is the thermal convection that developes during collapse inside the forming the proto-neutron star, and the complicated interaction of this convection with the neutrino flux (Burrows and Fryxell 1993; Wilson and Mayle 1988, Colgate, Herant, and Benz 1993). Turbulent motions can mechanically carry neutrinos from the deep, hot layers, and enhance the energy deposition into the outer layers. A self-consistent picture of this process is currently lacking. The neutrino flux is critical to the success of the mechanism, but the two-dimensional studies to date have included only a very rudimentary, quasi-one-dimensional treatment of the neutrino transport.
The numerical technique to treat this problem requires a combination of hydrodynamical and neutrino transport algorithms. During the last year, a multi-dimensional, multi-group neutrino transport code has been developed jointly by I. Lichtenstadt, our collaborator from the Hebrew University of Jerusalem, and A. Khokhlov, which is an adaptation of the hydrodynamical code ALLA. The code, ALLA-C, has recently been completed and tested and is ready for production runs. We will apply it to investigate several key aspects of the core collapse mechanism of supernova explosions.
The ionization and excitation equilibrium is computed assuming LTE, but the radiation field can and does depart radically from LTE due to scattering. Except at the very earliest times, supernova opacity is dominated by continuous electron scattering and by the scattering effects of many millions of spectral lines. Continuous absorption is important at very short (bound-free) and very long (free-free) wavelengths. The ratio of scattering to absorptive opacities can vary by more than six orders of magnitude, with corresponding changes in the coupling of the radiation field to the local sources. In a spherical atmosphere the radius of the photosphere is therefore a very strong function of frequency.
The standard explosion model will be W7, a white dwarf deflagration computed by Nomoto et. al. (1984). The input model is expanded hydrodynamically using a 1 dimensional code called BOOM, together with frequency averaged opacities computed using elements of SNAP and the line list (42 million lines) generated by Kurucz (1992). At var- ious points in the evolution from less than 1 day to ~40 days the output of the hydrocode is used as input to SNAP in order to compute the time dependent emergent spectrum of the explosion model.
The initial W7 model has 172 zones. This is re-zoned to generate a more suitable grid for radiative transfer calculations, usually providing 128 to 256 radial grid points. A similar number of angle points are generated by the solution of the characteristic equations defining curved 'rays' along which the comoving frequency is constant. The frequency grid is defined by the number of strong resonance lines which are treated explicitly. Previous experience has shown that up to 5000 lines may be necessary, generating a frequency grid with from 20,000 to 100,000 points. Fortunately the method of solution is linear with respect to the number of frequency grid points.
The solution requires ' 1 second per frequency point on a single CPU of a Cray YMP and runs at over 200 Mflops per CPU on this machine. SNAP is also highly parallel. It has achieved over 1.8Gflops on a YMP/8 and over 9Gflops on a C90/16. SNAP has recently been run on a T3D.
ALLA incorporates two direct summation algorithms for computations of the self-gravity of a star. One is based on the expansion of the Poisson integral in spherical harmonics, and is suitable when stellar deformations are moderate. This algorithm works in all geometries. The other is based on a Fast Fourier Transform convolution. It is suitable for arbitrary mass distributions, but works only in cartesian geometry. The advantage of direct methods is the fact that the solution is always guaranteed (in contrast to relaxation methods, no iterations are necessary). This part of the code takes ~ 20-40% of the hydrodynamical step.
For modeling thermonuclear supernovae, a degenerate matter equation of state is incorporated. This EOS inncludes contributions of ideal Feri-Dirac gases of electrons and positrons, non-ideal ions, and radiation. For core collapse simulations, the baryonic EOS is a Compressible Liquid-drop model for nuclei, and includes the effects of interactions and degeneracy of the nucleons outside nuclei (Swesty & Lattimer). The EOS is pre-computed in a tabular form with the internal energy, mass density, and chemical composition as independent entries. The EOS is called three times during each hydrodynamical step to evaluate pressure, adiabatic index, and other thermodynamic quantities. Using the internal energy instead of temperature as an independent variable, we avoid iterations of the equation of state. This results in a significant reduction of computing time. We split the internal energy and pressure on the 'cold' (T = 0), and 'thermal' counterparts (thermal = internal - cold), and then perform interpolation of cold and thermal quantities separately. The EOS interpolation is fully vectorized and parallelized. Time requirements are ~ 40% of the time to perform the hydrodynamical step.
Explosive burning of degenerate stars involves many (~ 102 ) nuclei from protons and neutrons to isotopes of iron group elements and ' 103 nuclear reactions. Implementation of such big nuclear reaction networks in multi-dimensional hydrodynamical computations is impractical. The nuclear kinetics is coupled to the hydrodynamics via three important quantities - the electron concentration (which is changed by weak interactions), the ion concentration (which is changed by strong and electromagnetic interactions) and the binding energy of matter (which is changed by both). In ALLA, the nuclear kinetics are described by model kinetic equations for these three quantities plus one equation for a leading reaction (C+C for CO mixtures). These equations accurately approximate all energetically important stages of nuclear burning and provide accurate enough energy release and rate of change of both electron and ion mole fractions necessary for hydrodynamical computations (Khokhlov 1991). Detailed chemical composition will be obtained by post-processing results with a large nuclear reaction network code.
In numerical simulations of deflagrations, one must guarantee the grid-invariance of the flame front. The flame speed relative to a fuel must not depend on the orientation of the flame on the grid, the grid cell size, nor fluid motions. These requirements are of crucial importance, especially if a feedback exists between the energy generated by the flame and hydrodynamical motions which advect the flame front. ALLA incorporates a simple flame capturing technique, which satisfies these requirements (Khokhlov 1993a). The deflagration flame is modeled by a single reaction-diffusion equation with the diffusion coefficient and reaction rate chosen such that the flame has a specified velocity D. This velocity is pre-computed separately taking into account ' 100 relevant nuclear species, ' 1000 nuclear reactions, and the thermodynamic and transport properties of stellar matter (Timmes and Woosley 1992). It is stored in a tabular form. The method is very robust, and capable of treating flame fronts with complex, varying topology. Coupled with the hydrodynamical code, it provides a good flame invariance even for very subsonic burning, while the PPM-based hydrodynamical code allows modeling of supersonic combustion as well.
The nuclear rates routines are based on an alpha chain reaction network developed by A. Khokhlov (Khokhlov 1991). This network is a 13-species (He, C, O, Ne, Mg, Si, S, Ar, Ca, Ti, Cr, Fe, Ni), 18-reaction network. The rate equations are integrated using a routine for solving multiple sets of simultaneous ODEs that takes advantage of the gather-scatter hardware of a Cray (Brun and Patnaik). This routine sorts the rate equations into three groups, depending upon the size of the loss terms, and solves these groups simultaneously via classical, asymptotic, or equilibrium methods as determined by the relation of the loss term timescale to the timestep. In this way no more time than necessary is spent on the rate equations of a particular cell.
The code inputs are the initial conditions for temperature, density, composition, and velocity; pressure, internal energy, total energy, and momentum are derived from these. A tabulated equation of state for degenerate electrons and the reaction rate coefficients are also read in. We choose background conditions (density, temperature, etc.) appropriate to a carbon-oxygen white dwarf, and start a detonation at one point. A perturbation is placed in the path of the detonation in order to initiate the creation of transverse waves. It is the interaction of the transverse waves that leads to the cellular structure we intend to observe.
A 1D version of EPIC is fully operational. The 2D version of EPIC is nearing completion; the hydrodynamics, equation of state, and I/O routines are finished, and the nuclear rates routines are being coded now. The immediate goal is to finish coding the nuclear rates routines in 2D, and to test the full 2D version of EPIC on some standard problems. We will then investigate the multi-dimensional structure of detonations in degenerate carbon-oxygen material.
Time requirements for the current version of the code are problem dependent. Usually, ~ 30 sub-cyclings are required for the lowest energy group, but this number rapidly decreases for high energy groups. The number of the neutrino-matter interaction subcyclings is usually kept to one or two per hydrodynamic time step. These numbers, however, may vary. >From 4 to 10 sec of C90 CPU time may be required for one full time step on the 128 x 64 grid.