1.   Introduction

As guest editor for the 1999 March-April issue of the new interdisciplinary (AIP & IEEE) magazine Computing in Science and Engineering (CiSE), the PI recently was afforded the opportunity to highlight the astrophysics community's accomplishments in the area of "Computational Cosmology" by assembling a set of expert theme articles on that topic. In that issue of CiSE, Bond et al. (1999) discussed progress that is being made with regard to measuring fluctuations in the cosmic microwave background (CMB); White and Springel (1999) reviewed how large N-body simulations have been used to effectively simulate the natural development of CMB-related fluctuations into nonlinear-amplitude, dark matter structures (galaxies); Bryan (1999) described efforts by the GC³ collaboration to properly simulate the dynamics of gaseous (baryonic) material as it flows into the dark matter structures; and, on behalf of the Sloan Digital Sky Survey collaboration, Szalay (1999) reviewed efforts that are underway to document in digital form a "complete" and detailed description of the universe's present nonlinear structure at optical wavelengths. After reading such a collection of articles on modern cosmological investigations, one gathers the distinct impression that a fairly accurate and unified picture is emerging regarding the process by which galaxies form.

In contrast to this, the astrophysics community is still struggling to find a broad, general explanation of the process by which stars form from the gas that resides in galaxies. (This, of course, also presents a monumental barrier to cosmologists because we won't really understand how galaxies form until the gas dynamical simulations of, for example, the GC³ collaboration, actually produce stars.) It is generally accepted, of course, that gravity is responsible for gathering low-density clouds of gas in the interstellar medium into clumps that are dense enough to be called protostars, and it is generally believed that once gravity begins to control the process, it happens fairly rapidly -- on the order of a dynamical time t_d ~ [ Gr₀ ]^-1/2 as measured by the mean density r₀ of the initial cloud state. (For interstellar molecular clouds with n ~ 10⁴ - 10⁶ cm^-3, for example, t_d ~ 10⁶ - 10⁵ years.) But there is still a debate as to how the process of dynamical collapse usually gets started, a debate over what sets the typical mass scale for stars -- a scale which must ultimately make sense in terms of observed properties of the stellar initial mass function (IMF) -- and related debates over what processes are principally responsible for the breakup of single clouds into binary or small multiple star systems, or into larger groupings of stars that we call clusters.

In the context of the formation of individual, low-mass stars and their associated circumstellar disks Shu, Adams & Lizano (1987) offer a compelling model in which collapse proceeds from initially centrally condensed, magnetically supported molecular cloud cores. (For related models and ideas along these lines, see Mestel & Spitzer 1956; Nakano 1979; Lizano & Shu 1989; Tomisaka et al. 1989; Basu & Mouschovias 1994; Li & Shu 1997; and Galli et al. 1999. But for a recent cautionary discussion, see Nakano 1998.) In the model as outlined by Shu, Adams, & Lizano (1987), a star's mass is not primarily determined by the physical properties of a molecular cloud core but, rather, by processes associated with the onset of nuclear fusion, envelope convection, and winds. But based on the clustering and binary statistics of stars in the nearby Taurus-Auriga region, Larson (1995) has found renewed support for the idea that it is the thermal Jeans mass in molecular cloud cores that sets the mass scale for young, low-mass stars. In this context it is interesting to note that, although it is widely acknowledged that magnetic fields play an important role in defining the observed properties of molecular cloud complexes, most research groups who have attempted over the past decade to simulate the process by which binary and small multiple star systems form have used conditions set by the thermal Jeans mass to define the onset of collapse (cf., Boss 1993, 1996, 1998; Boss & Myhill 1995; Burkert & Bodenheimer 1993, 1995; Matsumoto & Hanawa 1999; Monaghan 1994; Myhill & Kaula 1992; Nelson & Papaloizou 1993; Sigalotti & Klapp 1994, 1996, 1997; Truelove et al. 1997, 1998). Recently, however, Boss (1997, 1999) has tried to bring the two ideas together.

Given that star formation is ongoing in the solar neighborhood -- that is, relatively speaking, it is a process that can be studied in our own backyard -- it might seem strange to suggest that models of galaxy formation are more mature than are models of star formation. But the following few points illustrate why the star formation problem is in many respects the more challenging one:

• A model that satisfactorily explains how normal stars (r₀ ~ 100 g cm^-3) form from interstellar molecular gas clouds (r₀ ~ 10^-16 - 10^-18 g cm^-3) must cover a dynamic range of 20 (!) or more orders of magnitude in density whereas, even at the present epoch, the mean density of a typical galaxy is only 10 orders of magnitude larger than the mean density of the universe (and it was significantly less than 10 orders of magnitude during the epoch of galaxy formation).

• Whereas galaxies originate from linear-amplitude fluctuations in an otherwise homogeneous and isotropic medium, the molecular clouds from which stars form in the solar neighborhood are extremely clumpy -- that is, the clouds exhibit markedly nonlinear-amplitude inhomogeneities. Hence, it is difficult to base models of star formation in one-dimensional, linear perturbation theory. (As already mentioned, the clumpy structure of molecular clouds is almost certainly sustained by magnetic fields, which further complicates models of star formation.)

• Historically it has been very difficult to study protostellar configurations observationally except at the two density extremes, i.e., in the molecular cloud state or after the star has essentially formed. This is in part because star formation is a relatively rapid process so it is difficult to "catch" the system at intermediate states, and in part because the cloud cores where all the action is taking place are surrounded by dense layers of gas and dust that are optically thick to all but the longest wavelengths of electromagnetic radiation. And even within the solar neighborhood, one finds that either the relevant photon flux is very low at these longer wavelengths, or it has been difficult to resolve the relevant spatial scales observationally.

With the advent of large digital infrared array detectors and the development of larger and more sophisticated millimeter-wave radio telescope arrays, observational issues associated with this third point are being partially alleviated. In particular, millimeter-wave observations are providing sufficient spatial resolution and signal-to-noise to permit mapping of the structural and dynamical properties of star-forming gas clouds at relatively high volume densities and with linear scales approaching the size of our own solar system (cf., Sargent & Welch 1993; André, Ward-Thompson, & Motte 1996; Ohashi et al. 1997).

Faced with the tremendous dynamic range that separates the mean densities of molecular clouds from the mean densities of stars, and evidence that nearby star forming regions show structure at virtually all scales (Larson 1995), it would be foolish to believe that any single model that is based on a relatively simple geometry explains in full how stars form directly from conditions that are specified in molecular clouds. That's not to say that we haven't gained considerable insight into the process by which stars form by comparing observations to early spherically symmetric (e.g., Larson 1969) or more recent axisymmetric (e.g., Shu, Adams, & Lizano 1987) models of protostellar cloud evolution, but rather to emphasize that such models cannot possibly be painting for us the complete picture.

Taking advantage of the rapidly expanding computational resources that have emerged this decade, a significant number of research groups have attempted to model the star formation process in full three-dimensional generality, as mentioned above. But the dynamic range of such models is still necessarily limited by their finite numerical resolution so, at best, these models are also only able to paint for us a portion of the picture. Generally speaking (see Matsumoto and Hanawa 1999 for a recent overview), published models of cloud fragmentation have focused on the early stages of collapse when it is generally understood that protostellar clouds cool very effectively and, hence, "direct fragmentation" based on the thermal Jeans criterion is not difficult to achieve. However, as we have been reminded most recently by Boss (1999), at densities above 10^-13 - 10^-12 g cm^-3, protostellar clouds become heated by compression and this increase in temperature serves to discourage further fragmentation.

Over the past 15-20 years, the PI also has devoted a significant fraction of his research efforts toward modeling the dynamical evolution of protostellar clouds in full three-dimensional generality with an eye toward understanding why stars preferentially form in pairs. This has been motivated by the observational investigations of, for example, Abt & Levy (1976), Abt (1983), Duquennoy & Mayor (1991) and Mathieu (1994) which quite convincingly indicate that, as Mathieu states, ''binary formation is the primary branch of the star-formation process.'' The orbital periods of binary stars span an extremely wide range -- from fractions of days to thousands of years -- and they include systems with a wide variety of component mass ratios and a wide range of orbital eccentricities. Until we at least understand how binary systems with such a variety of physical properties are formed, we will be hard pressed to claim that we understand the star formation process.

Some of the PI's earliest simulations focused on the (lowest density) isothermal phase of collapse and demonstrated how direct fragmentation can produce binary stars with relatively long orbital periods (Tohline 1980; Bodenheimer, Tohline, & Black 1980). Subsequent work by numerous groups, as indicated above, have shown that it is possible to specify initial conditions in such a way that a binary system with practically any orbital period can be formed via direct fragmentation. But under realistic molecular cloud conditions, it seems most natural to expect that direct fragmentation will work effectively only at densities r₀ < 10^-13 - 10^-12 g cm^-3 -- i.e., below the point when heating due to adiabatic compression begins to occur -- and produce multiple systems whose orbital periods are of the same order as the dynamical time t_d that is associated with those densities -- i.e., greater than a few hundred years. (See the related discussion connected with Table 1 in § C.3.c of this proposal.)

With this in mind, in the late '80s and early '90s the PI's simulation efforts turned toward modeling the (higher density) adiabatic phases of protostellar cloud evolution in an effort to understand how binary stars with relatively short orbital periods form. One advantage of focusing on models with mean densities r₀ > 10^-12 g cm^-3 is that ionization fractions are expected to be extremely low, so magnetic fields are unlikely to significantly influence the cloud's evolution. Another is that the effective adiabatic exponent of the gas is large enough to permit the construction of configurations that are stable against further dynamical collapse. [In the protostellar models of Larson (1969) and Bodenheimer et al. (1990), for example, there is always a central core of material that is in hydrostatic balance. The high-density cloud core continues to contract fairly rapidly as it cools (and/or adds additional mass through accretion), but its subsequent contraction occurs quasi-adiabatically rather than dynamically.] If cloud cores of this type are to break into multiple pieces, however, it will have to occur as a result of an instability that arises spontaneously from an equilibrium state rather than via a Jeans-type instability that drives the "direct" fragmentation process.

In this context, progress has been made in following the nonlinear development of bar-like (often spiral-shaped), nonaxisymmetric instabilities in rapidly rotating, equilibrium gas clouds (Tohline, Durisen & McCollough 1985; Durisen et al. 1986; Williams & Tohline 1987, 1988). Also, building on ideas developed by Papaloizou & Pringle (1984), Goodman and Naryan (1988), and Adams, Ruden & Shu (1989), and tools developed by Hachisu (1986a,b), we have helped establish a better understanding of the behavior of gravitationally driven nonaxisymmetric instabilities that arise in relatively massive protostellar disks (Tohline & Hachisu 1990; Woodward, Tohline, & Hachisu 1994; Andalib, Tohline, & Christdoulou 1997). However, until very recently neither of these broad investigations into dynamical instabilities that arise in quasi-adiabatically evolving systems has developed into a theory of binary star formation that is competitive with models of direct fragmentation. As a result, recent reviews have concluded that direct fragmentation during the collapse of molecular cloud cores provides the best explanation for the formation of most binary stars (Boss 1993a; Bodenheimer, Ruzmaikina, & Mathieu 1993; Boss 1999).

However, as we explain in §C.2 of this proposal, improvements in modeling and simulation tools over the past few years have made it possible for the PI's group to examine the fission hypothesis of binary star formation in a more realistic manner than has heretofore been possible. We have done so in the context of the hypothesis as recently formulated by Lebovitz (1987) for quasi-statically contracting, inviscid protostellar gas clouds. Our results offer the most convincing evidence, to date, that binary stars can form during the quasi-adiabatic phase of a protostellar cloud's evolution by a process of fission, rather than via a Jeans-type fragmentation process. Because fission should work best at relatively high gas densities when protostellar gas clouds do not cool effectively under compression, it offers a mechanism for forming binary (or small multiple) star systems that complements the mechanism of direct fragmentation. Over the next few years we propose to extend our study of the structural and stability properties of quasi-adiabatically evolving protostellar gas clouds with an emphasis on the fission problem.

Many classes of interesting astronomical phenomena are known for which the proposed astrophysical scenarios involve dynamical or thermal mass-transfer in a close binary star system, sometimes leading to a common envelope phase and/or the complete destruction of the binary system through merger of its two components. Well-known examples of such systems are cataclysmic variables (CVs) [Warner 1995], low-mass X-ray binaries (LMXBs) [White et al. 1995; van Paradijs & McClintock 1995; Verbunt & van den Heuvel 1995], millisecond pulsars (Bhattacharya 1995), soft X-ray transients (SXTs) or X-ray novae [Chen, Shrader & Livio 1997], Algols (Batten 1989) and W Serpentis stars (Wilson 1989), contact binaries (W Ursae Majoris stars) [Rucinski 1985], some central stars of planetary nebulae (Iben & Livio 1993), double degenerate white dwarf (WD) binaries, perhaps leading to supernovae (SNe) of type Ia through a merger (Iben & Tutukov 1984), or subdwarf sdO, sdB stars (Iben 1990), and double neutron star (NS) or NS-black hole binaries, perhaps yielding gamma-ray bursts in a fireball when the system coalesces (Mészáros & Rees 1992; Ruffert et al. 1997; Kluzniak & Lee 1998).

In many cases, the present-day systems are known to be undergoing long-term, stable mass-transfer driven by angular momentum losses or nuclear evolution. Although such long-term, secular phases of evolution may be easier to study observationally, they are too slow to be amenable to present numerical simulations. But when a system becomes dynamically unstable toward mass-transfer, the transfer rate can rapidly reach 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 large-scale numerical simulations are required to further our theoretical understanding of them. Realizing that the numerical tools that we have constructed to examine fission and gravitational fragmentation processes in protostellar clouds and protostellar disks may in some cases be used effectively to simulate such phenomena, we have begun to examine tidal and mass-transfer instabilities in short period, compact binary star systems. As is documented in §G (Current and Pending Support) of this proposal, we recently have received partial support through NASA's Astrophysics Theory Program to pursue this investigation; the research is being conducted in collaboration with J. Frank [LSU] who is lead PI on the NASA project.

Although, in the context of this NSF proposal, the study of mass-transfer in evolved binary systems will not represent the principal focus of our research efforts over the next few years, we have decided to include a brief discussion of the subject here for several reasons. First, as is explained more fully in § C.2.e, below, our investigation into the relative stability against merger of close binary star systems began during the period of time that was covered by our most recent NSF award. Second, by demonstrating that our simulation tools can be used very effectively to follow through many orbits the motion of dynamically stable, close binary star systems, we have gained considerable additional confidence in the physical correctness of our simulations that address binary star formation processes. Finally, as we have grown to appreciate, dynamical processes that prove to be important in the context of evolved, close binary star systems often prove to be important in the binary formation problem as well. Hence, our plans to expand our studies of tidal and mass-transferring instabilities in evolved stars over the next few years will both overlap and complement the research activities being proposed herein.

2.   Results from Prior NSF Support

More than three years have passed since the starting date of award AST-9528424 from the MPS division of the NSF, and four years have passed since the PI last submitted a full proposal requesting renewed NSF support. Over this period of time, support from the NSF has been acknowledged in eighteen separate publications (7 refereed journal articles, 6 articles in conference proceedings, and 5 doctoral dissertations, as itemized by the "Jnn" reference notation in §D of this proposal). As itemized by the "Mnn" reference notation in §D, eighteen Quicktime movies also have been produced in an effort to illustrate more completely our various numerical simulation results. Our research activities have focused on the following five separate, but overlaping areas:

1. The construction of rapidly rotating, [Newtonian] self-gravitating, gas dynamical configurations that are in detailed force balance and therefore provide accurate, reliable, and realistic models for dynamical stability studies. [J11, J13, J16, J18; M3, M4, M5, M8, M9, M10, M11, M12, M16]
2. The study of global gas dynamics in galaxies with an eye toward understanding what mediates star formation in normal, gas-rich systems. [J1, J4, J5, J8, J9, J12]
3. The study of nonaxisymmetic instabilities in equilibrium self-gravitating gas dynamical systems in an effort to better understand how small multiple (primarily binary) star systems form. [J3, J13, J16, J18; M6, M7, M13, M14, M15]
4. The study of tidal and mass-transfer instabilities in close binary star systems. [J2, J7; M1, M2, M17, M18]
5. The development of accurate and efficient numerical techniques and algorithms for the simulation and visualization of gas dynamical systems in astrophysical settings. [J6, J14, J15, J17]

Over the past four years, five graduate students have completed their doctoral dissertation research under the PI's direction [student names and dissertation titles are cited in §D as references J2, J11, J12, J15, J16]. The PI also has supervised one postdoctoral research associate (D.M. Christodoulou; a position funded primarily through state-allocated funds to the LSU Department of Physics and Astronomy) and the research projects of three undergraduate students, one of whom (D. Sherfesee) was awarded a 1998 NSF Graduate Research Fellowship and is presently enrolled in Berkeley's graduate astronomy program.

Over the past four years, the PI also has been heavily involved in a significant science outreach project in the Parish of East Baton Rouge, Louisiana. He was lead investigator on a project that secured $115,000 in funds from the Louisiana State Board of Regents and $26,000 from the LSU Foundation to purchase and install a modern, 20" optical telescope in a new public park observatory building that was constructed at a cost of approximately $300,000 by the Park and Recreation Commission of the Parish of East Baton Rouge (BREC). Since the opening of "Highland Road Park Observatory" in November, 1997, in close collaboration with the Baton Rouge Astronomical Society (BRAS), LSU and BREC have hosted an open house and public viewing session every Friday evening; have hosted numerous monthly K-12 field trips; have run summer astronomy camps and teacher training sessions; and have provided an exceptional optical instrument for use by local amateur astronomers. For example, BRAS members have discovered over 30 new asteroids using the telescope and its CCD camera. (See URL http:// www.phys.lsu.edu/ observatory and [J10] for more details.) This science outreach project was initiated and completed without direct funding support from federal agencies, but is mentioned here as an indication of the scope of the PI's activities outside of normal research and university-level instruction.

In connection with the research area itemized as i, above, two particularly significant results have emerged. First, Andalib [J11] has developed a modified self-consistent-field technique for constructing equilibrium models of rapidly rotating, self-gravitating, gas dynamical configurations with compressible equations of state, nontrivial internal motions, and a variety of different nonaxisymmetric structures. To date, the technique has been applied only toward the construction of two-dimensional systems (i.e., either infinitesimally thin disks or systems having infinite vertical extent) with nonaxisymmetric structures, but we understand in principle how to extend the technique to fully three-dimensional systems (see further discussion in §C.3.b, below). Andalib has used the technique, for example, to construct infinitesimally thin disks with steady-state elliptical, dumbbell, and "binary" shapes, as viewed from a frame that is rotating with a specific system pattern frequency [J11, J13]. Invariably, as viewed from that same rotating frame, the fluid exhibits nontrivial internal flows [M3, M4, M5] -- somewhat analogous to the classical incompressible Riemann S-type ellipsoids. The Quicktime movie [M5] associated with Fig. 1 illustrates the differential internal flow in three of Andalib's steady-state models of centrally condensed, elliptical disks (one with flow that is entirely retrograde with respect to the figure motion; one that is entirely prograde; and one with two off-axis vortices separating a region of prograde flow from one of regrograde flow), and illustrates the flow in one of his models of an equal-mass, common-envelope binary. Although Andalib's models are only two-dimensional, his work represents a significant breakthrough in the sense that, historically, it has not been possible to construct equilibrium models of rapidly rotating nonaxisymmetric systems with nontrivial internal flows and equations of state that are sufficiently compressible to be of astrophysical interest in the context of protostellar clouds, gas-rich galaxies, and most normal stars. When it is not possible to construct individual equilibrium models of such systems, of course, it is impossible to diagnose the relative stability of such sytems or to say anything quantitatively meaningful about related evolutionary sequences.

Figure 1 Four Andalib Models [M5]

Second, Cazes [J16, J18; see also J13] has successfully constructed two homentropic, compressible (g = 5/3), fully three-dimensional, self-gravitating gas configurations that, in virtually all respects, appear to be compressible analogs of Riemann ellipsoids (CAREs) with fairly extreme principle axis ratios (10:5:2). Although both of Cazes' models have been constructed with a nonlinear hydrodynamic (rather than, say, a self-consistent-field) technique from initially unstable axisymmetric equilibrium models [M6, M7], both clearly have settled into nonaxisymmetric configurations that are in steady-state, in the sense that their structures are unchanging on a dynamical time scale, and they are dynamically stable [M8, M11]. These bar-like configurations spin about their shortest axis with a well-defined pattern speed and, as viewed from a frame rotating with the overall pattern, both models exhibit strongly differential internal motions [M9, M10, M12]. These internal flows are gratifyingly similar to the flows discovered by Andalib in 2-D disk models constructed with an SCF technique (discussed above), but in both of Cazes' models, a "violin-shaped mach surface" and a pair of weak standing shock fronts appear to be integral components of the steady-state flow. These results are significant because, as explained above, previous attempts to construct CAREs have been discouragingly unrewarding, yet understanding the properties of such configurations is extremely important to studies of the structure and evolution of protostellar gas clouds, rapidly spinning compact stellar objects, and self-consistent models of barred spiral galaxies.

In connection with the research area itemized as v, above, we are continuously working (with the manpower at our disposal) to improve our numerical simulation tools in order to be able to more accurately and efficiently simulate realistic astrophysical fluid flows. To date, our efforts have been focused on the accurate, three-dimensional representation and analysis of rapidly rotating, self-gravitating dynamical flows, while incorporating relatively simple equations of state and virtually ignoring the effects of magnetic fields, radiation transport, and relativistic fields. This approach has been reasonable because it is the competing effects of rotation and Newtonian self-gravity that are expected to dominate the dynamics of the systems in which we have been most interested, and it is still computationally impractical to treat all of these physical processes simultaneously in a realistic manner. Athough there are certainly improvements that remain to be made in our techniques and tools -- some of which are discussed in §C.3, below -- indications are that our present tools and techniques are competitive. Through the peer review process that has been established to evaluate proposals for high-performance-computing time on hardware maintained by NSF's national computing alliances, over the past two years we have been allocated and have utilized over 100,000 service units on Cray T3E platforms. We recently have submitted a proposal to NRAC requesting an additional 65,000 SUs for the coming year in support of our ongoing research activities.

Our primary simulation tool is a three-dimensional, finite-difference hydrodynamics code patterned after ZEUS-2D (Stone & Norman 1992). Our specific treatment of advection and source terms in the principal dynamical equations has been detailed by New [J2] and Cazes [J16]; our technique for efficiently solving the global Poisson equation has been detailed by Cohl [J6, J15, J17]; and the heterogeneous computing environment that we have developed to permit the routine visualization of complex, time-dependent fluid flows (see, for example, the Quicktime movie sequences [M1 - M18]) has been described by Cazes et al. [J14]. The mpi (message passing interface) version of our code -- recently designed and implemented by LSU graduate student Patrick Motl -- can now routinely simulate self-gravitating flows of the type discussed throughout this proposal on Eulerian grids containing 16 MegaCells. (Given the continuing upward trend in hardware, networking, and parallel compiler technologies, we comfortably expect this capability to increase by a factor of 4 - 8 over the next few years.)

Of all the improvements that we have made in our simulation tools over the past few years, the one that is most likely to make a long-term impact on the astrophysical community as a whole is the one reported in Cohl & Tohline [J17; see also J15]. As the abstract of that paper summarizes, we have discovered that an exact expression for the Green's function in cylindrical coordinates (v, f, z) is,

Q_m-1/2 is the half-integer degree Legendre function of the second kind, and the summation over the index "m" is from minus infinity to infinity. This expression is significantly more compact and easier to evaluate numerically than the more familiar cylindrical Green's function expression which involves infinite integrals over products of Bessel functions and exponentials. It also contains far fewer terms in its series expansion -- and is therefore more amenable to accurate evaluation -- than does the familiar expression for |x - x'|^-1 that is given in terms of spherical harmonics. This compact cylindrical Green's function (CCGF) expression is well-suited for the solution of potential problems in a wide variety of astrophysical contexts because it adapts readily to extremely flattened (or extremely elongated), isolated mass distributions. We now use it exclusively in our algorithm that calculates the value of the gravitational potential on the boundary of our computational grid; these values are then used as boundary conditions for our solution of the global Poisson equation in all of our dynamical simulations.

The idea that the quasi-adiabatic contraction of a rotating protostellar gas cloud may lead in a very natural way to the formation of binary stars through a process of "fission" was suggested over 100 years ago. (See, for example, the reviews by Lyttleton 1953; Chandrasekhar 1969; Durisen & Tohline 1985; and Lebovitz 1987.)

Figure 2 Rotating Fluid Drop
	Results from a 1992 space shuttle (STS-50) Drop Dynamics Experiment.

In a qualitative sense, this fission mechanism is envisioned to be similar to the process by which a drop of fluid breaks into two roughly equal pieces if it is induced to spin rapidly enough. (See the images in Fig. 2 from a drop dynamics experiment performed in the microgravity laboratory [USML-1] of the space shuttle.) In an effort to test this important hypothesis, we have slowly "cooled" one of Cazes' dynamically stable, rotating bar-like CARE models (described above in § C.2.b) in a very simple way, continuously following the model's dynamical flow in a self-consistent fashion throughout the evolutionary cooling phase [J16]. Specifically, in our polytropic equation of state, we have forced the polytropic constant to decrease linearly with time at a rate that would have reduced a spherically symmetric cloud to half its initial radius in 17.5 dynamical times. In effect, this reduced the specific entropy of the gas at a steady rate uniformly throughout the model. Although this rate of cooling is faster than would be expected during the adiabatic contraction phase of a real protostellar gas cloud, it was slow enough to permit the contraction to occur in a quasi-static fashion (i.e., the model remained in excellent virial balance throughout the cooling evolution) yet fast enough to permit us to complete the cooling simulation without demanding an unreasonable amount of high-performance computing resources.

Initially as the model cooled, it became even more elongated and its overall spin pattern frequency increased -- in concert with a shrinking principal moment of inertia -- but, it remained centrally condensed. After the polytropic constant dropped to 55% of its initial value, however, the overall configuration began to oscillate between two well-defined states: a highly elongated, centrally condensed bar, and an equal-mass, "common-envelope" binary. This phase of the slow cooling evolution can be most fully appreciated by viewing several of our Quicktime movie sequences [M13, M14, M15], but the two states between which the oscillation occurred are also illustrated here in Fig. 4 (see the last page of §C of this proposal). In the right-hand column of Fig. 4, the binary state is most clearly delineated by the nested 3-D isodensity surfaces (middle frame) and by the surface plot (bottom frame) of the effective potential in the equatorial plane, as viewed from a frame rotating with the overall pattern frequency of the configuration -- i.e., the orbital frequency of the binary. The effective potential (density profile) in the equatorial plane of the system is also illustrated by the solid (dashed) nested contour lines in the top frame of Fig. 4. (The solid circle identifies corotation.)

Notice the similarities between the effective potential shown on the right-hand-side of Fig. 4 and the familiar effective potential one derives in the binary Roche problem: two off-axis potential minima; two maxima associated with the "L4 and L5" Lagrange points; two saddle points (at which envelope material is paritally spilling out of the system) associated with the "L2 and L3" Lagrange points; and, finally, an "L1" saddle point separating the two potential minima. As the velocity vectors in the top frame of Fig. 4 also illustrate, when the model is in its binary state, a portion of the flow continues to stream along the full length of the "bar-like" configuration, but another portion of the flow is isolated around the separate binary components, effectively giving the components a net spin. These three physical features -- off-axis density maxima, close correspondence with the Roche potential, and fluid flow that is isolated around each off-axis component -- give us considerable confidence that the cloud configuration shown on the right in Fig. 4 is indeed a fully self-consistent (common-envelope) binary state.

Demonstration that a binary configuration does become available to a quasi-adiabatically contracting protostellar gas cloud as envisioned by proponents of the fission hypothesis of binary star formation over the past 100 years (see especially the modified scenario outlined recently by Lebovitz 1987) is the most important scientific result to have emerged from the PI's NSF-sponsored research efforts over the past few years. This result forms a principal component of Cazes' Ph.D. dissertation [J16] (degree awarded August, 1999), and is presently being written up for publication.

In connection with the research area itemized as iv, above, we have constructed equilibrium sequences of synchronously rotating, equal-mass binaries in circular orbit with a single parameter -- the binary separation -- varying along each sequence. Sequences have been constructed with various polytropic as well as realistic white dwarf and neutron star equations of state. Then, using our Newtonian, gravitational hydrodynamics code, we have examined the dynamical stability of individual models along these equilibrium sequences [J2, J7].

Our simulations indicate that no

Figure 3 Binary Mass-Transfer [M17]

points of instability exist on the sequences we analyzed that had relatively soft equations of state (polytropic sequences with polytropic index n = 1.0 and 1.5 and two white dwarf sequences). As is well illustrated by one of our Quicktime movie sequences [M1], such systems remain stable even when they share a sizeable common envelope. This result gives us added confidence that the common-envelope binary depicted on the right-hand-side of Fig. 4 can be a long-lived, stable state. We did, however, identify dynamically unstable binary models on sequences with stiffer equations of state (n = 0.5 polytropic sequence and two neutron star sequences). The evolution of one such unstable system is illustrated in another of our Quicktime movie sequences [M2]. We have thus inferred that binary systems with soft equations of state are not driven to merger by a dynamical instability. For the n = 0.5 polytropic sequence, the separation at which a dynamical instability sets in appears to be associated with the minimum energy and angular momentum configuration along the sequence. Our simulations suggest that, in the absence of relativistic effects, this same association will also hold for binary neutron star systems.

Very recently, graduate student P. Motl has begun to conduct a similar investigation involving unequal-mass polytropic binary star systems. This investigation is both computationally more challenging and broader in physical complexity than the study of equal-mass systems. For example, for a given polytropic index and component separation, there is an infinite range of component mass ratios from which to choose; and even when the system is dynamically stable against a global tidal instability, it may be unstable toward an entirely different type of mass-transferring instability. One of our Quicktime movie sequences [M16] illustrates the degree to which we have been able to follow a detached system through more than four full orbits. Two other movie sequences [M17, M18] (one frame of which is shown here as Fig. 3) illustrate the early phase of a mass-transfer instability in a system with a mass ratio q = 0.88.

3.   Proposed Research

Cazes' [J16] demonstration (see the discussion in §C.2.d, above) that a quasistatically contracting gas cloud can spontaneously evolve to a binary state breaths new life into the fission hypothesis of binary star formation, as most recently outlined by Lebovitz (1987), and thereby represents a milestone in modern star formation research. As explained earlier, because it can operate in high density cloud cores that heat up under compression, fission offers an alternate but complementary mechanism for the formation of binary stars to the Jeans-type process of direct fragmentation which works most effectively in relatively low-density clouds. Since the binary state arises spontaneously from an otherwise dynamically stable, equilibrium configuration, this mechanism for forming binary stars brings with it the hope (expectation) that the mass-spectrum of binary systems with orbital periods less than a few hundred years can be explained by processes unrelated to the spectrum of initial fluctuations in interstellar clouds.

However, with the tools in hand and the time available, Cazes was only able to study one slow cooling evolution, and that particular evolution did not actually end with the cloud in a binary state. After watching the system oscillate back and forth several times between the two states depicted here in Fig. 4, Cazes stopped cooling the model in order to determine which configuration represented the absolute minimum energy state at that particular point in the evolution. As the Quicktime movie [M13] associated with Fig. 4 illustrates, in response to this query the system settled back into the bar-like configuration. But we are confident that, had the cooling been continued, and had Cazes been able to maintain adequate spatial resolution to accurately follow the cloud's continued radial contraction, the system would have eventually settled into the binary equilibrium state. Clearly this speculation needs to be verified if the fission hypothesis is to be fully resurrected. Other examples, and examples employing more realistic cooling processes, need to be given as well.

We should point out that the classical idea that fission can only produce binary systems with equal-mass components is almost certain to be dispelled by carrying to completion this type of simulation. As we appreciate from models of mass-transferring binary systems, when the components do not share a common envelope, the equal-mass configuration is not preferred. Hence, although the fission process may initially create two equal-mass off-axis condensations, we suspect that as those condensations individually contract (via cooling processes) and begin to separate from one another, an equal-mass configuration also is unlikely to be preferred.

This is one principal avenue of research that we propose to follow over the next few years. Employing the more efficient mpi-based hydrodynamics algorithm, mentioned earlier, that has recently been developed in the PI's group and an adaptive grid technique like the one proposed by Berger & Colella (1989) and implemented by Truelove et al. (1997,1998) and Bryan (1999), we propose to follow Cazes' cooling evolution with higher sustained grid resolution, further in time to document how and when the system settles permanently into the binary state. As a second example, we propose to slowly cool either the second of Cazes' CARE models (his Model B), or a new 3-D CARE constructed via a self-consistent-field technique (see § C.3.c, below). Third, in order to determine whether or not the process of fission is sensitive to the rate at which the cloud contracts, we propose to repeat at least one of these simulations but introduce a rate of cooling that is two or three times slower that the rate employed by Cazes. Finally, as discussed more fully below, we propose to incorporate a radiation transport algorithm into our primary hydrodynamic simulation code in an effort to model spatially differential cooling processes more realistically.

Because they are very broadly applicable to studies of the structure and evolution of protostellar clouds, compact stellar objects, and gas-rich galaxies, we also propose to pursue further the development of techniques for constructing compressible analogs of Riemann ellipsoids (CAREs) with a wide range of physical characteristics. We expect to build directly on the self-consistent-field technique that has been successfully applied by Andalib to 2-D systems, extending his technique to fully three-dimensional configurations. In order to illustrate how this might be done, we must first outline certain features of Andalib's 2-D technique.

In a frame of reference that is rotating with a constant angular velocity kW_f, Euler's equation governing the motion of a barotropic gas may be written in the form (cf., Lynden-Bell & Katz 1981; Papaloizou & Savonije 1991),

where H is the gas enthalpy, F is the gravitational potential, v is the gas velocity as measured in the rotating frame, z º Ñ ´ v is the fluid vorticity, and v is cylindrical radius. For a steady-state system in which the gas is confined to an infinitesimally thin disk, the left-hand-side of Euler's equation may be replaced by,

where, making use of the continuity equation and following the lead of Papaloizou & Savonije (1991), we have related the product of the surface mass density S and velocity to a pseudo-stream function Y through the expression,

Hence, in such a steady-state system it must be possible to express the vortensity [ ( z + 2 W_f ) / S ] (the inverse of what Lynden-Bell & Katz 1981 refer to as the "load") as a function only of Y (Papaloizou & Savonije 1991) and then move the term that is customarily on the left-hand-side of Euler's equation to the right-hand-side and group it with all the other terms under the gradient operator. As a reasonable first pass at the problem, Andalib [J11] assumed that the vortensity was either uniform in space or, at worst, only a linear function of Y, i.e.,

in which case Euler's equation can be transformed into the following relatively simple algebraic expression defining the spatial relationship among the various key physical variables:

Following the procedure laid out in traditional self-consistent-field techniques (Ostriker & Mark 1968; Hachisu 1986a,b), therefore, in Andalib's technique it is this algebraic expression that must be satisfied simultaneously with the selected barotropic equation of state and with the Poisson equation relating the gravitational potential to the surface density in order to construct a steady-state, equilibrium configuration.

But in the case of nonaxisymmetric flows, an additional mathematical expression must be used to constrain the functional form of Y in a manner that properly satisfies the continuity equation. Based on the definition of vorticity and the relationship given above between the momentum density Sv and Y, Andalib [J11] realized that it is the following self-adjoint PDE that governs the spatial behavior of the pseudo-stream function in the case where the vortensity is, at worst, a linear function of Y:

Andalib has devised an iterative algorithm that successfully solves this second-order PDE in concert with the other three traditional SCF equations for a variety of different nonaxisymmetric surface boundary conditions, a variety of different degrees of gas compressibility (polytropic indices 0.1 < n < 1.3 ), and a variety of different values of the constants C₀ and C₁. It is with this tool that he has been able to build the various 2-D, nonaxisymmetric equilibrium models described above.

In extending the technique to nonaxisymmetric structures with finite vertical extent, we expect to be able to break the problem down into a series of vertical slices perpendicular to the figure's overall spin axis, in a manner similar to that used by Eriguchi & Hachisu (1985). We should be able to solve for the structure of each 2-D slice using Andalib's technique and adopt slice-to-slice variations in the various "constants" of the flow such that the slices blend together vertically in a physically reasonable manner. In order to ascertain what is physically reasonable, we will take certain cues from the properties of the two CAREs that have been constructed by Cazes using dynamical simulation tools, as well as from the somewhat simpler "Dedekind-like" configurations that have been examined recently by Uryu and Eriguchi (1999). For example, both Cazes [J16, J18] and Uryu & Eriguchi (1999) have noted that, locally throughout the fluid, the vorticity vector should be slightly missaligned with respect to the spin axis of the figure. Also, as we move from 2-D to 3-D structures, it will be important to relate the pseudo-stream function to the load, as defined by Lynden-Bell & Katz (1981), rather than to the vortensity which has been defined earlier only in terms of a surface density.

Finally we note that, although Andalib successfully introduced gas compressibility into his equilibrium models using a polytropic prescription for the equation of state, he was unable to converge to models with polytropic indices n > 1.3. Based on Cazes' results, we suspect that it will be necessary to introduce a "violin-shaped mach surface" and standing shock fronts into model configurations that have relatively large degrees of compressibility. This will be a nontrivial, but we think not insurmountable, task. For example, in Chapter XII of Landau & Lifshitz's (1959) classical text on Fluid Mechanics, a discussion is presented on how to use the method of characteristics to determine the location and strength of weak shocks in curved flows that resemble in many respects the flows observed in Cazes' CAREs. As long as the shocks are weak, treating the flow as homentropic will introduce relatively little error in the determination of the primary flow variables because the entropy jump across a weak shock goes as the third power of the (small) pressure jump. This approximation should simplify the model construction task, then afterward an estimate of the timescale for secular evolution of the converged model could be acquired by accounting for the actual (small) amount by which entropy is created by each passage of the fluid through one of the standing shocks.

Observational data from arrays of millimeter-wave radio telescopes now provide sufficient spatial resolution and signal-to-noise to permit mapping of the structure and dynamical properties of star-forming gas clouds with linear scales approaching the size of our own solar system (cf., Sargent & Welch 1993; Ohashi et al. 1997). Hence, there may be an opportunity to directly compare the properties of our models of quasi-adiabatically contracting protostellar cloud cores -- both as CAREs and as common-envelope binaries -- with the observed properties of star-forming clouds having comparable scales. To illustrate this point, Table 1 (drawn from reference [J16]) scales one of Cazes' dimensionless 3-D, bar-like CAREs to a cloud of one solar mass having five different sizes ranging between 1 and 100 AU. These represent a possible protostellar cloud core configuration at different stages of its contraction (mean densities between 10^-13 and 10^-7 g cm^-3). Although Cazes' model employs a simple n = 3/2 polytropic equation of state, when scaled to the sizes given in Table 1 it predicts mean cloud temperatures and values of the total angular momentum that are quite realistic.

During the first year of this project, we propose to incorporate radiative transfer techniques and realistic equations of state into our modeling algorithm along the lines of those developed and implemented by Boss (1984), Bodenheimer et al. (1990), Boss & Myhill (1992), and Sigalotti (1998). In addition to helping us simulate the slow contraction of our gas clouds more realistically, as mentioned above, such techniques will permit us to produce spatially resolved surface brightness and velocity maps from our dynamical models for comparison with radio-frequency observations of star forming regions.

Although directed primarily toward a better understanding of how binary star systems are formed, as has been discussed in § C.1.b of this proposal most of the PI's recent modeling efforts have been of a sufficiently general nature that they also can provide insight into the structural and stability properties of other self-gravitating systems, such as galaxies or highly evolved stars with or without accompanying accretion disks. In particular, if our CARE models are scaled to the size of compact stars, as long-lived nonaxisymmetric configurations they could be a source of continuous-wave gravitational radiation. We also propose, therefore, to determine (in a post-Newtonian approximation) what the observational signature of such a continuous-wave source of gravitational radiation would be. Building on the techniques outlined by Finn & Evans (1990) and Rasio & Shapiro (1992), New [J2] has already incorporated into the PI's basic hydrodynamic algorithm the tools that are necessary to compute this signature. What remains to be determined is the most realistic size scale to which our CARE model should be scaled and whether it would be expedient to incorporate a more realistic equation of state into the model.

In connection with these plans to compare our steady-state and evolving cloud models to observations of molecular cloud cores and compact stellar objects, the PI has arranged to spend his upcoming sabbatical leave (Spring, 2000) at Caltech, working in collaboration with both Anneila Sargent and Kip Thorne. As Director of the Owens Valley Radio Observatory (OVRO), Sargent leads a strong radio astronomy group whose efforts are largely focused on studies of star forming regions in our Galaxy. By spending an extended period of time interacting with this radio astronomy group, the PI hopes to gain a much better appreciation of the variety and quality of data that is being collected in connection with ongoing star formation processes in the solar neighborhood, and can effectively begin to establish a relationship between his models and the observations. Under Kip Thorne's direction, Caltech is widely recognized as the home base for one of the world's leading theory groups whose efforts are directed largely toward modeling sources that should be detectable by the laser interferometer gravitational-wave observatory (LIGO). By participating in the activities of this group on a regular basis while on leave, the PI hopes to achieve the second observationally-related objective as outlined above. No NSF funds are being requested in support of the PI's sabbatical leave at Caltech during the Spring of 2000. It is mentioned here, however, because the PI's planned research activities during this leave period will significantly complement the research activities for which NSF support is being requested.