1. Summary of Proposed Research, Including Objectives and Goals

Over half of all stars in the sky are multiple star systems and most of these are binary stars. We recently have completed a review of the mechanisms that have been proposed to explain the origin of binary star systems (Tohline 2002). Approximately one half of the binary stars are close enough to one another for mass to be exchanged between the components at some point in their evolution (Trimble 1983). The components of these interacting binaries influence one another in a variety of ways including tidal forces, irradiation, and in the most extreme manifestation of the interaction, through mass transfer and mass exchange. As a result, the evolution and observational properties of the binary components and the binary itself are dramatically different from the behavior expected in the absence of mass transfer, even when this transfer proceeds at a leisurely pace.

Interacting binaries of one type or another, especially those in which the accretor is a compact object (Frank, King & Raine 2002), play central roles in many of the astrophysical sources of high energy photons, and are involved in the generation of gravitational and high-energy particle radiation. Most of the interacting binary types known observationally are presently undergoing long-term, stable mass-transfer on timescales dictated by secular angular momentum losses (e.g. gravitational radiation and magnetic braking) and/or nuclear and/or thermal evolution. Well-known examples of the above are cataclysmic variables (CVs) [Warner 1995, Beuermann 2002], low-mass X-ray binaries (LMXBs) [White et al. 1995; van Paradijs & McClintock 1995; Verbunt & van den Heuvel 1995], soft X-ray transients (SXTs) or X-ray novae [Chen, Shrader & Livio 1997], Algols (Batten 1989, Vesper Honeycutt & Hunt 2001; Zavala et al. 2002) and W Serpentis stars (Wilson 1989), contact binaries (W Ursae Majoris stars) [Rucinski 1989; Qian 2001], super-soft X-ray sources [Greiner 2000, Rappaport, Di Stefano & Smith 1994], and the ultra-luminous X-ray sources in nearby galaxies [King et al. 2001]. In these systems, the fraction of the donor star's mass that is transferred during one orbit is typically » 10^-12 to 10^-9, many orders of magnitude less than what current numerical 3-D hydrocodes can resolve. All of the above systems, however, must have descended from binaries in which the accretor of today was initially the more massive component who evolved first off the main-sequence. The mass transfer in these progenitor systems was in many instances dynamically unstable, yielding mass-transfer rates many orders of magnitude above the currently observed rates and leading in some cases to a common envelope phase (see e.g. Warner 1995; Verbunt & van den Heuvel 1995; Nelson & Eggleton 2001) or a thermal mass-transfer phase (King & Schenker 2002).

In addition, there is a wide class of binary star systems not presently undergoing mass-transfer for which the astrophysical scenarios that have been proposed to explain their origin or ultimate fate involve dynamical or thermal mass-transfer in a close binary. Examples of such systems are millisecond pulsars (Bhattacharya 1995), some central stars of planetary nebulae (Iben & Livio 1993), double degenerate white dwarf binaries, perhaps leading to supernovae of type Ia through a merger (Iben & Tutukov 1984), or subdwarf sd0, sdB stars (Iben 1990), and double neutron star binaries, perhaps yielding g-ray bursts in a fireball when the neutron stars coalesce (Paczynski 1986, Mészáros & Rees 1992, Ruffert et al 1997, Mészáros 2001). The evolutionary scenarios that are drawn upon to explain the existence of some of these systems call for events in which 10% or more of the donor's mass can be transferred during a single orbit.

We shall focus on these relatively rapid phases of mass-transfer, investigate the conditions required for stable/unstable mass-transfer, and follow numerically their evolution in unstable cases. The dynamical mass-transfer that ensues in unstable binaries rapidly reaches nonlinear amplitudes and may lead to the formation of a massive disk-like object, maybe a common envelope phase (Iben & Livio 1993), some mass loss from the system, and in some cases, end in a merger. These phenomena are truly three-dimensional in character and require large-scale numerical simulations to further our theoretical understanding of them. While the direct observation of these phases – with the possible exception of gravitational wave emission from mergers – remains unlikely in the near future, these investigations will help us understand the origin of these binaries and their mutual relationships. We plan to study both the tidal and the mass-transfer stability of close binaries in which the component stars obey a polytropic equation of state (EOS) but, in general, possess nonuniform entropy. As we argue below, these simple structural approximations enable us to treat, with some limitations, a very diverse set of objects including main sequence (MS) stars, WDs and NSs.

We distinguish two types of dynamical instability that can occur in close binaries:   1) A dynamical instability which arises in close proximity of extended binary components with sufficiently stiff EOS due to the presence of tidal forces, even before contact and onset of mass-transfer; and   2) A dynamical instability which arises upon contact, even when total angular momentum and mass are conserved, as mass is transferred from one component to the other. We shall refer to these instabilities as the tidal and mass-transfer instabilities, respectively.

The classical tidal instability, first studied by Darwin (1887), deals with the secular evolution away from synchronous rotation due to dissipation. More recent treatments of this secular instability include: an analysis of tidal stability of coplanar binaries (spin axes perpendicular to orbital plane) by Counselman (1973), who shows that in unstable systems evolution away from synchronism is accompanied by a secular decrease or increase in the orbital separation; and a more general treatment by Hut (1980) allowing for arbitrary orientations of spin and orbital angular momenta.

In a series of papers Lai, Rasio & Shapiro (1993a; 1993b; 1994a; 1994b) studied the equilibrium and stability of binaries with polytropic EOS, making use of analytic techniques and a variational energy method. They were able to extend their work to unequal masses and allow for internal motions. They showed that in addition to the classical secular instability, a new dynamical instability exists not far from the point of secular instability, which leads to rapid binary inspiral under certain conditions. They showed that the instability occurred before contact was reached for stiff EOS (see also New & Tohline 1997).

In binaries with unequal components, in general one of the stars will reach contact with its Roche lobe first (the donor) and mass-transfer to its companion star (the accretor) will commence either because the donor star is evolving or because the binary separation is reduced by systemic angular momentum losses. What happens next, i.e. the stability of binaries to mass-transfer, depends crucially on the reaction of the donor star to mass-loss, and on the changes in Roche geometry that ensue from the mass-transfer.   Important factors determining the evolution of the Roche geometry are: 1) whether all the mass transferred by the donor is accreted (conservative mass-transfer); and   2) if some of the mass transferred is lost from the system, what are the consequential angular momentum losses (Webbink 1985; King & Kolb 1995). If upon contact and onset of mass-transfer, the radius of the donor and the radius of the Roche lobe vary in such a way that contact becomes deeper, mass-transfer will proceed at an ever increasing rate and the situation is unstable. If, however, the contact becomes shallower, the donor will detach and mass-transfer will cease. In the latter case, if some physical effect other than mass-transfer insists on driving the system back into contact, stable mass-transfer may be possible for as long as the driving is present. As mentioned above, in CVs and LMXBs, systemic angular momentum losses drive stable secular mass-transfer from a less massive donor onto a compact object.

It has become customary to characterize the behavior of the donor and the Roche lobe upon mass-transfer by introducing the following parameters (Webbink 1985, Ritter 1996): z_R, z_e, and z_s, which are respectively defined as the logarithmic derivatives of the donor's Roche radius, its thermal equilibrium radius and its adiabatic radius, with respect to the mass of the donor. So, for example,

It can be shown (Webbink 1985, Ritter 1996) that mass-transfer is dynamically (or adiabatically) stable if z_s > z_R, whereas thermal stability requires z_e > z_R. Stars with deep convective envelopes, such as lower MS stars and giants have z_s < z_e, indicating that the dynamical instability sets in first, whereas stars with radiative envelopes, such as upper main sequence stars, have z_e < z_s. In the latter case, if z_R < z_e < z_s, it is in principle possible for the star to shrink rapidly first, avoiding a dynamical instability but then expand slowly by thermal relaxation yielding a stable mass-transfer at a rate controlled by thermal relaxation.

MOVIE-1 Stable, Equal-Mass Binary

Quicktime (1,380K)

In stable mass-transfer situations, the depth of contact is relatively shallow, on the order of a few pressure scale heights, and does not change much with time. In those situations it is correct to adopt an isothermal atmosphere approximation in order to model the mass-transfer and the structure of the stream. Recent simulations by Blondin (1998) confirm that such a model can reproduce the behavior of the stream predicted by Lubow & Shu (1975; 1976). In our case, however, we are interested in dynamically unstable situations which lead to deeper degrees of contact and therefore a polytropic approximation for the envelope is better. In the absence of nonlinear effects, which we expect to limit the growth (see below), a polytropic atmosphere of uniform entropy will yield a divergent mass-transfer rate in a finite time (e.g. Webbink 1985).

Before outlining what evolutionary scenarios are likely to develop as we follow the nonlinear evolution of systems that are found to be unstable toward either a tidal or mass-transfer instability, it is important to point out that ...

In an effort to illustrate our early successes with this aspect of our investigation, we offer two movies which may be accessed via the web (if you're reading this proposal on-line) by clicking on the "Movie-1" icon shown here [URL: www.phys.lsu.edu/ astro/ cazes/ movies/ stable.db.mov], or the "Movie-4" icon as discussed in §3.b. Movie-1 [taken from New & Tohline 1997] shows that, for a sufficiently soft equation of state, an equal-mass (q º M_donor/M_accretor = 1) binary system can be stable against the tidal instability even when the two stars are so close to one another that they are sharing a common envelope. Movie-1 follows the evolution of one such system through 3.5 orbits. Although the system is being viewed with its orbital plane face-on, orbital motion is not apparent because the evolution is being shown in a rotating reference frame. Movie-4 demonstrates how well we are able to follow through more than five orbits the evolution of an unequal-mass (q = 0.84) system in which the less massive component (on the right in the Movie-4 icon) very nearly fills its Roche lobe. The movie displays equatorial-plane mass densities on a logarithmic color scale, as viewed from a frame of reference rotating with the binary system's initial orbital frequency. Additional animation sequences from this evolution can be accessed through URL: casa.colorado.edu/ ~motl/ research/ binary/ stable_binary_visc.html.

We expect that, depending on the assumed conditions, some binaries may become stable after a certain mass excess is transferred, some may exchange mass, and still others may evolve through a common envelope phase. This schematic division is adopted for simplicity since there may well exist cases where a combination of the above processes operate sequentially.

We anticipate that if the mass ratio is not extreme and the polytropic indexes are similar, it should be possible for the binary to transfer some mass and regain stability. In systems that have a mass ratio q < 1, for example, the binary separation should increase with time as material flows from the donor to the accretor. But, based on the 1-dimensional linear stability theory, if q is not too far from unity the Roche lobe should nevertheless shrink around the donor and thereby promote a continuous mild transfer event. Once q drops to a sufficiently small value, however, shrinkage of the Roche lobe will be discontinued and the mass transfer event terminated. If the donor has a polytropic index n = 3/2, linear theory predicts that mass transfer should stop when q £ 2/3.

When the donor is more massive than the accretor (q > 1) and has a deep convective envelope, the mass-transfer almost certainly will be dynamically unstable. This is because, as mass is transferred from the more massive star to the less massive one, the orbital separation should actually decrease with time while the radius of the donor is increasing. It is generally thought that once the rate of transfer is high enough, the accretor is swamped and driven out of thermal equilibrium, and the matter transferred forms a common envelope. Most calculations of the subsequent evolution have dealt with a relatively compact star (the accretor) orbiting or plunging into the convective (polytropic) atmosphere of a giant star (Rasio & Livio 1996; Terman, Taam & Hernquist 1994; 1995 and references therein). Some calculations specifically focus on the ejection of the envelope (Taam et al. 1997; Sandquist et al. 1998).

The reasons why it is important to investigate the onset and early stages of a common envelope evolution (CEE) are discussed by Rasio & Livio (1996): 1) the dynamical phase is probably the only phase that can be calculated accurately within the capabilities of current computers; 2) the results of this early phase can be used as an initial state for follow-up calculations using ad hoc approximations; and 3) in some cases the fate of the system is resolved during the dynamical phase.

Movie-2 shows results from our first partial run of an n = 3/2 polytropic system in which the more massive star is the donor. Specifically, the illustrated system has an initial mass ratio q = 1.32 and, at the beginning of the simulation, the donor (shown on the left in the Movie-2 icon) is just marginally filling its Roche lobe, but eventually it develops a well-defined mass-transfer stream. It is interesting to note that, in this case, the center of mass of the system (and, hence, the orbital axis for the system) is positioned to the left of and L1 Lagrange point and is actually inside the donor envelope. Clicking on the Movie-2 icon [URL: casa.colorado.edu/ ~motl/ reseach/ binary/ movies/ driven/ streamq.mov] will bring up a Quicktime movie that follows the evolution of this system through approximately three orbits. The evolution is shown in a rotating frame of reference, however, so the primary orbital motion is not apparent. Additional movie sequences and accompanying diagnostic plots from this simulation can be obtained from the URL: casa.colorado.edu/ ~motl/ research/ binary/ driven.html.

MOVIE-2 Mass Transferring Binary

Quicktime (4,805K)

While it would be interesting to compare the results obtained for common envelope evolutions using Smooth Particle Hydrodynamics (SPH; the technique used by Rasio & Livio) and Finite-Difference Hydrodynamics (FDH; our selected technique) our intent is to explore some basic questions in a parameter regime quite different from what most CEE calculations assume. Therefore we are - at least initially - not aiming to reproduce the calculations of Rasio and Livio using a FDH code. In future developments (beyond the immediate scope of this proposed project) we would expect to incorporate more physics in such calculations (e.g., turbulent viscosity, radiation transport and convection). While this is relatively straightforward in FDH, it is not so in SPH.

Ultimately, the incorporation of additional physics – in particular, radiative transport or an ad hoc method of cooling the envelope material – is expected to be crucial to accurately model the CEE. Nevertheless, our proposed adiabatic simulation of the CEE represents, to our knowledge, the first attempt to simulate the dynamical evoluton of two well-resolved and self-consistent binary components into a common envelope phase. As such, it will serve as an important basis for future research in this area.

The coalescence of double white dwarf (WD) and double neutron star (NS) binaries is important to examine since this process may produce a number of astrophysically interesting objects and events. Double white dwarf binary mergers have been suggested as precursors to some Type Ia Supernovae (Iben & Tutukov 1984; Iben 1988, 1991; Iben & Webbink 1989; Branch et al.  1995) and to long g-ray bursts (Katz & Canel 1996). WD-WD mergers may also lead to the formation of massive single white dwarfs or neutron stars (Colgate & Petschek 1982; Saio & Nomoto 1985; Iben & Tutukov 1986; Kawai, Saio, & Nomoto 1987; Chen & Leonard 1993); to the formation of subdwarf stars; or to the formation of hydrogen deficient, highly luminous stars (Iben 1990 and references therein; Webbink 1984). NS-NS mergers may result in bursts of g-rays and neutrinos, in the production of r-process elements, and in the formation of black holes (BH) (Eichler et al. 1989; Meyer 1989; Narayan, Paczynski, & Piran 1992; Rasio & Shapiro 1992; Davies, Benz, Piran, & Thielemann 1994; Katz & Canel 1996; Lipunov et al. 1995; Ruffert, Janka, & Schäfer 1996; New & Tohline 1997; Swesty, Wang & Calder 2000; Rosswog, Davies, Thielemann & Piran 2000). Merging compact binaries are also expected to be relatively strong sources of gravitational radiation (Abramovici et al. 1992; Cutler et al. 1993; Thorne 1995).

MOVIE-3 Binary Merger

mpeg (701K)

Most of the numerical simulations of mergers that have been carried out in recent years have been focused on the NS-NS or NS-BH problem because of their direct relevance to LIGO (and its sister instruments, like TAMA and VIRGO); and most, although not all (cf. Rasio & Shapiro 1994; Zhuge et al. 1996; Rosswog et al. 2000), have focused on equal-mass (q = 1) systems. Even when relativistic effects are ignored in such simulations, the stars usually merge on a dynamical timescale because one star (or both simultaneously, in the equal-mass case) encounters a tidal instability before either star fills its Roche lobe (Rasio & Shapiro 1994; New & Tohline 1997). This occurs, in turn, because the EOS of a NS is relatively stiff. As Rasio & Shapiro (1994) have summarized, the existence of a Roche limit and any associated mass-transfer instability has little importance for most NS star binaries because this "dynamical instability always arises before the Roche limit." When mass transfer begins, the system is already in "a state of dynamical coalescence," so the mass-transfer rates can be much larger than 10% per orbit. One such evolution, taken from New & Tohline (1997), is shown here as "Movie-3" [URL: www.phys.lsu.edu/ astro/ merger.mpg]. Although we are interested in continuing to pursue this type of merger problem, the focus of our proposed research will be on systems in which the secondary star has a relatively soft EOS and the star fills its Roche lobe before the system encounters a tidal instability.

2. Computational Methodology/Algorithms

In order to accomplish the stated objectives of this project, we plan, first, to use a 3D self-consistent field (SCF) technique to generate accurate equilibrium structures of synchronously rotating binaries in circular orbits. Each rotationally flattened and tidally distorted binary model will then serve as input into (initial conditions for) a dynamical simulation. Each dynamical simulation will be performed using a 3D, finite-difference hydrodynamic (FDH) technique that follows in a self-consistent fashion the motion of every fluid element in both stars. Our present SCF and FDH codes, as well as forerunners of these codes, have been effectively used to study gravitationally driven instabilities in protostellar gas clouds (Durisen et al. 1986; Woodward, Tohline & Hachisu 1994; Cazes & Tohline 2000); rapidly rotating, compact stellar objects (New, Centrella & Tohline 2000; Lindblom, Tohline & Vallisneri 2001, 2002); and, as mentioned above, close binaries with equal mass (New & Tohline 1997). Over the past few years, we have made key modifications in both of these algorithms in order to ensure that they perform with a sufficiently high level of numerical accuracy to permit us to confidently model mass-transfer instabilities in unequal-mass binaries. A detailed description of both codes and their capabilities is given in Motl, Tohline, & Frank (2002). Here we briefly summarize the principal elements of both.

It is important that we have the capability to construct accurate equilibrium structures as initial input into each dynamical simulation so that the simulations are not dominated by dynamical transients that depend heavily on initial conditions. We use the iterative SCF technique described by Hachisu (1986a,b) in its most general form in order to construct unequal-mass polytropic binaries with the internal structure of both stellar components fully resolved. Models are constructed on a uniformly zoned, cylindrical mesh that conforms exactly to the mesh that is used in the FDH code (see below) so that no interpolation is needed when the converged initial state is introduced into the dynamical code. To date, we have used this SCF technique to generate a wide variety of synchronous, polytropic binaries with mass ratios 0.1 < q < 1.5 in which both stars have uniform (but different) entropy distributions (homentropic). We can readily build detached, semi-detached or, in the limit of identical components, contact systems. These models have been constructed on computational grids ranging in size form 128³ to 256³ as necessary to conform to the chosen resolution of the FDH code. As measured by the global virial error, these models typically converge to better than one part in 10⁵. As a further test, we have compared parameters from our equilibrium binaries with those calculated from the point mass Roche model by Mochnacki (1984). As expected, our results differ from Mochnacki's because we use a self-consistent gravitational potential as opposed to the point mass potential as is appropriate in the Roche model. Our results converge to Mochnacki's in cases where the mass ratio becomes small and/or the separation increases.

For unequal-mass binaries, we have discovered that it is not possible to construct an initial model in which the center-of-mass of the system lies precisely on the coordinate axis of the cylindrical grid. Since the FDH simulations are performed in a frame of reference that rotates with the orbital frequency of the binary and the coordinate axis of the grid serves as the axis of rotation of the grid, it is then necessary to impart a very small velocity "kick" to every fluid element in the initial model in order to compensate for this discrepancy and minimize motion of the center-of-mass through the grid (for details, see Motl et al. 2002).

In addition to polytropic binaries, we are able to construct self-consistent, Newtonian models of compact binaries by choosing an appropriate equation of state such as the zero temperature WD equation of state (Hachisu 1986) or any of several realistic NS equations of state as described by New & Tohline (1997). Binary models using homentropic equations of state are, however, inadequate for examining the stability of mass transfer in main sequence stars without deep convective envelopes (stars with masses greater than but on the order of one solar mass). For an homentropic polytrope of index n, z_s = (1 - n)/(3 - n) so that z_s = -1/3 when n = 3/2 for example. On the other hand, the proper main-sequence mass-radius relationship demands z_s » 1. By allowing the specific entropy to be a function of the density, we are able to construct equilibrium configurations that are both convectively stable and that reproduce the desired mass-radius relationship.

Our primary dynamical simulation tool is an mpi (message passing interface) based algorithm designed to perform multidimensional, gravitational CFD (computational fluid dynamic) simulations efficiently on parallel computing platforms. As has been described in detail by Motl et al. (2002), the evolution of a system is governed through a finite-difference representation of the multidimensional equations that describe the dynamics of inviscid, self-gravitating, compressible gas flows. More specifically, the employed finite-difference hydrodynamics (FDH) algorithm is based on the second-order accurate, van Leer monotonic interpolation scheme described by Stone and Norman (1992), but extended to three dimensions on a uniform, cylindrical coordinate mesh. Most of our simulations, to date, have been performed on a computational lattice containing 130´98´256 grid zones in the R, Z, and q directions, respectively, with no assumed geometric symmetries. But over the past year, utilizing the capabilities and capacity of Blue Horizon (the SP3) at the SDSC, we have begun to utilize a computational lattice with twice this resolution in each dimension: 258´194´512. The fluid equations are advanced in time via an explicit integration scheme so an upper limit to the size of each integration time step is set by the familiar Courant condition. At each time step, the Poisson equation is solved implicitly using a combined Fourier transformation and iterative ADI scheme (Cohl, Sun & Tohline 1997) in order to determine, in a self-consistent fashion, how fluid elements are to be accelerated in response to the fluid's own self-gravity. As is described more fully, below, boundary values for the Poisson equation are determined via a compact cylindrical Green's function technique (Cohl & Tohline 1999). By performing our simulations in a frame of reference that is rotating with an angular velocity W₀ equal to the initial angular velocity of the binary orbit, we are able to minimize the effects of numerical viscosity (see New & Tohline 1997; Swesty, Wang and Calder 2000), which invariably plays a role in any FDH algorithm.

We have constructed a heterogeneous computing environment through which we are able to routinely create volume-rendered images and 3D animation sequences of our evolving fluid configurations. Relatively small amounts of data are regularly passed across the network from the parallel machine that is performing the FDH simulation to a machine that is more suited to visualization tasks; images that will become part of an extended animation sequence are created from this data; then the 3D data sets are immediately destroyed in order to minimize the impact on disk storage. We are utilizing the scriptable attributes of Maya (formerly Alias-Wavefront) software on SGI systems in order to perform these critical visualization tasks. (The 3D "Movie" sequences referred to throughout this proposal have been produced in this fashion. A variety of other animation sequences that have been produced with this tool can be accessed on-line via the URL: www.phys.lsu.edu/ faculty/ tohline/ image.ensemble/ movies.html.) This visualization tool and the heterogeneous computing environment in which it functions will provide a critical component of the simulations being proposed here.

This is the computationally intensive segment of a multi-year research project that is designed to address many of the scientific objectives that have been outlined here in §1, and thereby contribute to a more complete understanding of the tidal and mass-transfer instabilities in close binary systems. In this context, we propose to carry out the following well-defined simulation tasks.

• Task I: Determine whether homentropic, polytropic binaries are dynamically unstable to a mass-transfer instability, as predicted by linear theory, when the less massive component is the donor star and the mass ratio q is close to unity. If any of the modeled systems is found to be unstable, follow its nonlinear evolution to determine whether the mass-transfer event is conservative.

• Task II: If the nonlinear phase of mass transfer is conservative, follow its evolution to completion to ascertain whether mass-transfer shuts off and the system restabilizes at a new (smaller) value of q, or if merger occurs.

We have laid the groundwork for these investigations by (a) constructing numerous equilibrium binary sequences with q close to but not equal to unity; and (b) following through more than five orbits the dynamical evolution of two detached systems and showing that these systems are perfectly stable against mass transfer (see Motl et al. 2002 and the discussion here associated with Movie-4).

• Task III: From equilibrium model sequences having q > 1 (i.e., systems in which the more massive star is the donor), select several models that are unstable toward mass transfer and follow their evolutions far enough in time to determine whether or not (a) mass-transfer rates are generally larger than in the q < 1 systems and (b) the mass-transfer event is violent enough to lead to a common-envelope evolution.

We have laid the groundwork for these investigations by following through several orbits the evolution of one n = 3/2, q = 1.32 model that we expected to be rather violently unstable toward mass transfer (see the discussion associated with Movie-2).

• Task IV: Construct and follow the evolution of at least one system whose stars initially have a nonuniform entropy distribution, in an attempt to more closely represent the structural properties of solar-type MS stars.

3. Preliminary Results and Progress to Date

Over the past five and a half years, we have been awarded a total of approximately 390,000 SUs at the SDSC and UT-Austin to study dynamically evolving astrophysical systems in connection with our NSF- and NASA-sponsored research into the birth and destruction of binary star systems. Our present allocation of 130,000 SUs on SDSC's SP3 (Blue Horizon) runs through October 31, 2002. Our current request for PACI resources will build directly upon this successful, ongoing NPACI project.

i) Parallel CFD Algorithm

Our initial simulations on the T3E at the SDSC were carried out with a version of our code that had been developed on other parallel platforms and written in high-performance-fortran (HPF). Using the available Cray HPF compiler (PGHPF, developed by the Portland Group), we found that relatively few unexpected code modifications were required in order to successfully compile and execute our code on the SDSC system. We discovered, however, that the PGHPF compiler was producing executable code that performed computational instructions on each T3E processor with extraordinarily low efficiency. In 1998, we successfully rewrote our entire FDH algorithm, incorporating explicit message passing instructions via mpi. We found that this new, mpi-based code executed four to five times faster than the old HPF version of the code! According to numbers drawn from "pat," a performance monitoring tool that ran on the T3E, we determined that the code's single-processor execution rating was 36 MFlops. We also found that the code scaled very well as the number of processors was increased from 4 to 128. (Details of the code's performance on the Cray T3E can be found in several of our previous NRAC proposals, which are available on-line: 1999, 2000, 2001.)

Over the past three years, we have turned our attention to problems associated with mass-transfer in binary star systems during the late stages of stellar evolution, as described in this proposal. Over the past year and a half in particular, we have invested considerable effort fine-tuning the capabilities of our SCF and FDH algorithms so that they will be able to perform the required modeling tasks with a high degree of quantitative accuracy. Motl et al. (2002) describe in detail how well our code presently performs on a variety of benchmark simulations. In one simulation, for example, a detached binary system with mass ratio q = 0.84 was evolved through just over five orbital periods on a grid having 130 × 98 × 256 zones in R, Z, q respectively. This evolution is illustrated by Movie-4 [URL: casa.colorado.edu/~motl/research/binary/movies/stable_binary_visc/eq.mov ], which shows how the equatorial-plane density distribution of the system changes with time. Other aspects of this system's evolution as modelled by the current version of our FDH code are illustrated in the middle two frames of Fig. 1 (see proposal p. 14). The solid curves in the left-hand frame of Fig. 1 show as a function of time the volumes occupied by material inside three different isodensity layers in the secondary star (the star closest to filling its Roche lobe and shown on the right in the Movie-4 icon), with the top-most solid curve tracing the stellar surface; the dashed curve shows as a function of time the volume enclosed by the secondary's self-consistently determined Roche lobe. All of the curves are essentially straight, horizontal lines, which means that the star (and its companion, which influences the Roche lobe volume) remains in hydrostatic equilibrium to a high degree of accuracy.

MOVIE-4 Stable, Unequal-Mass Binary

Quicktime (8,451K)

From this same benchmark evolution, the center-of-mass trajectories of both stars and of the system as a whole are shown in the middle, right-hand frame of Fig. 1. In this figure frame, which has been drawn at a scale (-1 < x < +1; -1 < y < +1) to include the entire mass of the system, the three separate trajectories appear to be small dots with little or no discernible structure. (When plotted in the inertial reference frame on this scale the trajectories of the two stars are indistinguishable from circles.) This illustrates that, even after five orbits, the centers-of-mass of the two stars and of the system as a whole essentially have not moved from their initial positions. This provides strong confirmation that our SCF code produces excellent initial states and that the hydrodynamical equations are being integrated forward in time in a physically realistic manner. In the bottom three plots of Fig. 1, we have magnified a small region around each of the center of mass trajectories – expanding the linear scale of the previously described plot by approximately a factor of 50. These magnified views reveal that, although it is very small, there is measurable motion of the centers of mass. In the plot on the left that shows the magnified trajectory of the primary, we also have shown the size of our radial grid spacing, DR = 7.87 × 10^-3. This indicates the characteristic size of our discretization and emphasizes how very small the motion of each center of mass is; the motion of all three centers has been confined within a single grid cell through five full orbits. In the magnified plots, the trajectory of the center of mass of each individual star shows a small oscillatory motion in the x-direction. This is very low-amplitude, epicyclic motion which tells us that the orbits of the stars are circular to the level of two parts in 10⁴.

The results shown in the top two frames of Fig. 1 display information that is directly analogous to the middle two frames, but the displayed results come from an earlier version of our SCF and FDH codes; a version only slightly different from the one used by New & Tohline (1997) to examine the stability of equal-mass binaries against tidal instabilities. As the two figures from this evolution illustrate, with this earlier version of our simulation tools, we observed a slow expansion of the secondary star; the orbit itself developed a significant epicyclic amplitude; and after about three orbits, the Roche lobe was encroaching on the surface of the secondary. A comparison of the top two frames of Fig. 1 with the middle two indicates that significant improvements have been made over the past few years in our ability to accurately model the evolution of binary stars. These include: Taking advantage of the message-passing interface (MPI) libraries on parallel computers so that the number of azimuthal zones could be doubled – from 128 to 256 zones over the full 2p radians; SCF iterations are performed in double precision to obtain better initial states; as mentioned earlier, a small velocity is imparted to the stars as they are introduced into the FDH code in order to compensate for the slight off-axis initial location of the system center of mass; the FDH algorithm was modified to make the integration scheme more properly time-centered; a more accurate scheme has been developed to determine the gravitational potential on the boundary of the computational grid (see § 3.a.ii); and artificial viscosity was included to mediate a variety of weak shocks that appear at the surface of the stars.

These improvements have put us in a position to study mass-transfer instabilities in unequal-mass binaries with a degree of quantitative accuracy that heretofore has not been possible. We not only have improved on the simulation capabilities of New & Tohline (1997), but as Motl et al. (2002) detail, the epicyclic amplitude of our binary orbits are an order of magnitude smaller and angular momentum conservation is an order of magnitude better than in the recently published equal-mass binary simulations of Swesty, Wang, and Calder (2000). As Motl et al. (2002) have concluded, with our current simulation tools we should be able to follow mass-transfer events in which the fraction of the donor's mass that is transferred is greater than a few × 10^-5 per orbit for approximately 20 orbits. In particular, we are in a position to accurately determine mass-transfer instability boundaries of the type defined elsewhere in this project proposal.

We also have followed the evolution of several semi-detached systems in order to begin our investigation of both stable and unstable mass transfer. An up-to-date listing of the simulations (with corresponding links to movies of each) can be found at URL casa.colorado.edu/~motl/research/binary/binary.html. For example, as described earlier in connection with Movie-2 (see §1.c.ii) we have followed the evolution of a semi-detached system in which the donor is the more massive component and therefore the mass transfer rate should rapidly grow with time. It is generally difficult to generate initial conditions for semi-detached binaries in which the contact is deep enough so that the resultant mass transfer can be calculated correctly. To generate such initial conditions we introduced driving by artificially removing orbital angular momentum at the level of about 2% per orbit. We expected that once the system is driven into sufficiently deep contact and the mass transfer is growing, we might then be able to remove this artificial driving and continue the evolution without it. These simulations did eventually lead to a merger of the components of the binary but also revealed a subtle problem with the conservation of angular momentum which is being investigated.

4. Justification of Resources Requested

i) Required Platforms

To complete the tasks outlined for this project, we are requesting an allocation of 480,000 SUs on the IBM SP3 at the SDSC (see Table 2 for details). From our past work at the SDSC and elsewhere we feel that this platform is quite suitable for our simulation work for the following reasons. First, the mpi implementation of our FDH code has been developed and thoroughly tested on the SP3. As illustrated by the numbers in Table 1, our FDH code scales fairly well on the SP3 architecture. Second, the SP3 has sufficient capacity to permit us to simulate systems with significantly higher spatial resolution than previously has been possible (or at least practical). Third, we have a small SP3 system at LSU, so numerous tests can be conducted on this local machine in direct concert with production runs at SDSC. Finally, we have developed a heterogeneous computing environment that smoothly integrates the capabilities of the SP3 with our available visualization tools.

ii) Code Requirements and Optimization Issues

Since we have discussed code optimization issues fairly thoroughly in earlier sections of this proposal, we will not comment further on them here. From the information provided above -- in §3.a.ii and Table 1, for example -- we can estimate how much computing time will be required to make significant progress on this project. Typically, each binary simulation must be followed through 5 or more orbits before the physical outcome of the evolution can be ascertained. Utilizing 128 processors on the SP3, and a grid resolution of 258²´512, each evolution will therefore demand 885 hours, that is, 1.25 months (!) of wall-clock time to complete; and this assumes that we have uninterrupted access to 128 processors. This is not very practical. Instead, we propose to utilize 256 (or 512) processors which, according to the numbers shown in Table 1, will demand only 95 wall-clock hours = 4 days (or 59 hours = 2.5 days) per orbit. Thus, each 5-orbit evolution can be completed in as few as 20 (or 12) days of dedicated execution. In a single year, we can realistically expect to complete only a few evolutions of this magnitude. Quite a few 'trial' evolutions can be completed on lower-resolution grids, however.

Our proposed simulations at a resolution of 258²´512 will require at least 32 GB of RAM, hence, at least 128 processors on the SP3. More specifically, according to the AIX "size" command, if we utilize 128, 256, or 512 processors, our simulations will require 0.25 GB/proc., 0.14 GB/proc., or 0.08 GB/proc., respectively. As just discussed, due to execution time constraints, we will most likely run these high-resolution simulations on 256 or 512 processors. Our estimated SU requirements (calculated below) assume we will use 512 processors.

iii) Proposed Production Runs

In §2d of this proposal we described the simulation tasks that must be finished in order to satisfactorily complete this multi-year project. We will postpone Task IV to future years and will focus this year on Tasks I, II and III. As Table 2 indicates, in connection with each task we plan to carry out several runs at low spatial resolution (requiring 64 or fewer compute nodes) and for only a few orbits each in order to determine which initial state is most worthy of an extended, high-resolution run. The number of service units required to complete each high-resolution run has been estimated assuming we utilize 512 compute nodes. (If we use fewer nodes, somewhat fewer SUs would be required but, as explained earlier, the wall-clock time required to complete each evolution becomes prohibitive.)

The Task numbers that have been identified in connection with each proposed simulation in Table 2 follow the outline used in our earlier discussion of the respective problems in §2d. The column labeled N_models specifies the number of different initial model configurations that will be examined and the column labeled N_orbits lists the number of orbits through which we estimate we will need to follow each binary system in order to complete the specified task. Along with each Task we have specified the (R´Z´q) grid resolution that will be used and the approximate number of SUs that will be required to simulate one complete orbit at that grid resolution (see §3.a.i for details). Finally, the number of SUs that we estimate will be required to complete each proposed simulation task (last column of Table 2) is obtained by multiplying the required number of SUs/orbit by N_orbits and by N_models. Our Table 2 total is approximately 480,000 SU. This is the amount of computing resources that we are requesting over the coming twelve-month allocation period.

v) Special Requirements

Our primary simulation tools (the SCF code and the FDH code) have been developed entirely by our group at LSU and require no special packages other than standard mathematics and FFT utilities on the SP3. In the past, we have relied heavily on visualization software running on the SDSC Vislab facilities, but this will no longer be a requirement because the necessary software (specifically, Maya) now resides on our local machines at LSU.

5. Local Computing Environment

Through LSU's Concurrent Computing Laboratory for Materials Simulation we have routine access to a dual-processor SGI Octane and an SGI Onyx that drives an ImmersaDesk system. On this SGI equipment we have installed Maya (TM), the Alias-Wavefront software that will be used to routinely render 3D images and generate movie sequences from our dynamical simulations, as explained in §2.c. LSU operates a 32-processor IBM SP3 that is dedicated to research computing. On this machine we are able to run dynamical evolutions of modest size and length. We expect to continue to utilize this machine to complete a variety of short test runs, with the idea that full production runs will be carried out on NPACI facilities.

This past year, the Louisiana legislature funded a new, statewide information technology initiative through which LSU received $7M in new, recurring funding. The PI (Tohline) has been serving as interim director of the Center for Applied Information Technology and Learning (LSU CAPITAL) that has been established at LSU through this initiative. A portion of this year's funding has been used to purchase a 1024-processor, Beowulf-class supercomputer: 1.8 GHz Intel P4 Xeon processors; myrinet network; 1 GByte RAM per processor (see www.lsu.edu/ lsucapital/ beowulf.html for details). The hardware for this system arrived in May, 2002 and, within the past week (1^st week of July), the system has been installed and integrated successfully. Early linpack benchmarks indicate that the system as a whole clocks at 2 TFlops. In preparation for the arrival of this machine, two students from our research group (Shangli Ou and Mario D'Souza) attended a workshop at NCSA in October, 2001 that was designed to introduce users to the Pentium III-based Beowulf cluster that had just been installed at NCSA. During this workshop our gravitational CFD code was successfully ported to the NCSA cluster. This has positioned us well to take advantage of the new cluster that is being installed at LSU. Early, short test runs indicate that, on a given number of physical processors (16, for example) our code executes at least as fast on LSU's P4 Xeon-based system as it does on the SP3 at SDSC. We expect, therefore, that the system at LSU will ultimately be a tremendous resource for our group. Since the hardware system has not yet been thoroughly tested, however, we are not sure how robust it will be and, in particular, we are not sure how well our code will scale to large numbers of nodes or will perform on it in a production mode. It is therefore important that we continue to secure access to the SP3 at the SDSC to support our research program.

6. Other Supercomputer Support

Over the past 2 years, the PI has collaborated with Thorne and Lindblom (Caltech) to model the nonlinear growth of an unstable "r-mode" in isolated, rotating neutron stars, and to simulate the tidal disruption of a neutron star that is in orbit about a more massive, Kerr black hole. Through this collaboration, the PI has gained access to the HP Exemplar at CACR. However, the intent of this allocation (to Thorne) has been to provide a mechanism by which students in Thorne's group can begin to develop the numerical tools necessary to study these general relativistic processes. The request being made here for NSF PACI resources is intended to support our own primary investigation into mass-transferring, nonrelativistic binary systems; none of the simulations being proposed herein will be performed on CACR machines. We should note that our initial application of the PI's simulation tools (as described here in §2) to the "r-mode" problem have proven to be very successful (Lindblom, Tohline, & Vallisneri 2001, 2002). The results of this work (and images produced with our visualization tools) were featured in the April-June 2001 issue of NPACI & SDSC's EnVision Magazine.

7. Project Team Qualifications

A major component of the PI's research work over the past twenty years has involved the large-scale numerical simulation of multidimensional gas flows in astrophysical systems utilizing a variety of different high- performance- computing platforms. In the early 1980s, the PI held a postdoctoral research position in group T-6 at Los Alamos National Laboratories at which time he broadened his exposure to state-of-the-art CFD techniques and learned to develop efficient algorithms for vector supercomputers. From the mid-1980s through the early 1990s, Tohline was PI on a large account at the Cornell Theory Center during which time he and his students experimented with the parallel implementation of CFD algorithms on the IBM 3090 and various accompanying FPS Array Processors. When an 8K-node MasPar MP-1 system arrived at LSU in 1991, the PI abandoned the CTC facilities because his group was having considerable success utilizing the MasPar system for its largest simulation efforts. As has been detailed elsewhere in this proposal, in recent years the PI and his students have gained experience developing gravitational CFD algorithms for execution on a variety of parallel platforms including MasPar's MP-2, the CM-5, the Cray T3E, and the IBM SP2 & SP3. The group also has developed and implemented a heterogeneous computing environment (at LSU and at the SDSC) that tightly couples its parallel simulation efforts to sophisticated visualization tools.

Co-PI J. Frank has broad experience in studies of the astrophysics of interacting binary systems, particularly when such interactions involve fluid flow and accretion events. He is, for example, lead author of a widely referenced technical CUP publication entitled, "Accretion Power in Astrophysics." His expertise within this project helps ensure that the results of the proposed complex multi-dimensional and nonlinear simulations are physically correct and quantitatively meaningful when placed in the context of the extensive literature on this subject. Frank also is lead investigator on the NASA-funded, Astrophysics Theory Program project in partial support of which this request for NPAC resources is being made. (Tohline is Co-PI on this NASA-funded project.)

Shangli Ou and Mario D'Souza are graduate students enrolled in the Ph.D. physics program at LSU. They both will be heavily involved in the simulation work outlined in this proposal. Both students are entering the fourth year of their graduate studies and will be deciding the focus of their dissertation research based in large part on our success with this year's simulations of mass-transferring binaries, as described in this proposal.

	A1:	Bibliography
	A2:	Publications Resulting from NPACI Research
	B1:	Curriculum Vitae for J. E. Tohline
	B2:	Curriculum Vitae for J. Frank