The dynamical evolution of dense star clusters is a problem of fundamental importance in astrophysics, but many aspects of the problem have remained unresolved in spite of several decades of theoretical work and many recent improvements in the quality of observational data (especially with the Hubble Space Telescope and the Chandra X-ray Observatory). Starting around 1998, my students and I embarked on a new theoretical study of dense star cluster dynamics using a parallel supercomputer code based on Monte Carlo techniques for solving the Fokker-Planck equation in spherical symmetry. Although the initial focus of our work was on globular clusters, our code is very general and has been used more recently for the study of other systems, such as young, dense star clusters in galactic nuclei and starbursts. Our code was developed and tested thoroughly during the first few years of this project [61] [74], and was used in several initial studies of systems containing single stars only. In particular, we studied cluster lifetimes in the Galactic tidal field and their dependence on the stellar initial mass function [74], and we investigated the development of the “Spitzer mass-segregation instability” in two-component systems [60].

A key feature of our new study, and the main present goal of this work, is the inclusion of primordial binaries in the dynamical calculations. The realization in the early 1990's that primordial binaries are present in globular clusters in dynamically significant numbers has completely revolutionized our theoretical perspective on these systems. Most importantly, dynamical interactions between hard primordial binaries and other single stars or binaries are now thought to be the primary mechanism for supporting a globular cluster against core collapse [88]. In addition, exchange interactions between primordial binaries and compact objects can explain very naturally the formation of large numbers of X-ray binaries and recycled pulsars in globular cluster cores [57]. Resonant interactions involving primordial binaries can also result in dramatically increased stellar collision rates [94]. Direct observational evidence for stellar collisions and mergers in globular clusters comes from the detection (mainly with HST) of large numbers of blue stragglers concentrated in the dense cluster cores (see below).

The main goal of our current work is to design Monte Carlo simulations of cluster dynamics treating both the cluster itself and all relevant interactions self-consistently, including direct physical collisions and dynamical interactions involving primordial binaries. This idea is particularly timely because the latest generation of parallel supercomputers now makes it possible to perform such simulations for a number of objects equal to the actual number of stars in a large cluster (N ~ 10⁵ – 10⁶ for a globular cluster; up to ~ 10⁸ for galactic nuclei). Most previous studies of globular cluster dynamics incorporating primordial binaries have used direct N–body simulations. These simulations have been limited to very small N ~ 10³ – 10⁴. Special-purpose computing hardware, such as the Japanese GRAPE, has recently made possible direct N–body simulations with up to N ~ 10⁵ stars, although the computing time for the largest simulations with binaries can be many months. In contrast, our parallel Monte Carlo code allows us to perform self-consistent dynamical simulations of clusters containing N ~ 10⁵ – 10⁶ stars with ~ 10% binaries in typically less than a day of computing time. Even higher binary fractions, which may be more realistic for globular clusters [98] can now be handled very efficiently through direct integrations of all 3–body and 4–body encounters [94].

A closely related new project, done in collaboration with Dr. Marc Freitag (ARI, Heidelberg, Germany) concerns the formation of massive black holes in dense clusters through runaway collisions and mergers of young, massive main-sequence stars during the first ~ 10⁷ yr of dynamical evolution [90], [92]. Runaway collisions may provide the most natural path for the formation of “intermediate-mass” black holes in globular clusters and super star clusters. Another related new project, led by Northwestern postdoc Natasha Ivanova, studies the effects of dynamical interactions in dense clusters on the stellar evolution of binary star populations [98].

Postdoctoral fellows who worked with me on this include Simon Portegies Zwart (while a Hubble Fellow at MIT), and Northwestern postdocs Natasha Ivanova and John Fregeau. Three graduate students have been associated with this work: Kriten Joshi (MIT PhD '00), John Fregeau (MIT PhD '04), and Atakan Gurkan (Northwestern PhD '05). Several undergraduate students (MIT UROP, NSF REU, and senior thesis students) have also been involved, and in many cases co-authored some of our papers: Eric Ford [56], Boris Zbarsky [56], Wes Watters [60], Cody Nave [74], Pokman Cheung [94]. This work has been supported by several grants from the NASA Astrophysics Theory Program and the NSF Stellar Astronomy and Astrophysics program.

Coalescing compact binaries with neutron star or black hole components (NS-NS or NS-BH) are the most promising sources of gravitational waves for detection by current laser interferometers, such as LIGO and VIRGO. Most calculations of gravitational wave emission from coalescing binaries have focused on the waveforms emitted during the last few thousand orbits, as the frequency sweeps upwards from about ~10Hz to ~1000Hz. The waveforms in this regime can be calculated fairly accurately by performing high-order post-Newtonian (PN) expansions of the equations of motion for two point masses. However, at the end of the inspiral, when the binary separation becomes comparable to the stellar radii (and, ultimately, when the two stars merge), hydrodynamics becomes important and the character of the waveforms must change. In advanced LIGO, special purpose narrow-band detectors that can sweep up frequency in real time may be used to try to catch the corresponding final burst of gravitational waves. In this terminal phase of the coalescence, the waveforms contain information not just about the effects of relativity, but also about the interior structure of a NS and the nuclear equation of state (EOS) at high density. Extracting this information from observed waveforms, however, requires detailed theoretical knowledge about all relevant hydrodynamic processes. Since 1992 I have been developing both numerical and analytical tools that have allowed us to develop this knowledge, in preparation for the arrival of the first gravitational-wave data. Additional motivation for this work is provided by the many theoretical models of gamma-ray bursts (GRBs) that rely on coalescing compact binaries to provide the ~ 10 ⁵³ erg of energy required to power bursts observed over cosmological distances.

The final merger of the two stars is driven by a combination of relativistic effects and global hydrodynamic instabilities, which can drive the binary system to rapid coalescence once the tidal interaction between the two components becomes sufficiently strong. Using numerical hydrodynamic calculations in Newtonian gravity, Rasio & Shapiro [10, 21, 27] * demonstrated for the first time the existence of these global instabilities for double NS systems. Our results were later confirmed by many other groups, sometimes using very different numerical approaches. A detailed analytic study of dynamical and secular instabilities in compact binaries was also undertaken (extending the classic work of Chandrasekhar for an incompressible fluid), using a new energy variational method [11, 14, 17]. This analytic work confirmed the importance of hydrodynamic instabilities for the terminal evolution of coalescing compact binaries.

Until recently, most hydrodynamic calculations of compact binary coalescence had been performed, for simplicity, in Newtonian gravity. However, it is clear that the coalescence of NS-NS and, even more NS-BH binaries, will be strongly influenced by PN corrections to gravity. In particular, we have shown analytically that PN effects tend to make compact binaries even more unstable than their Newtonian counterparts [46]. Northwestern postdoc Josh Faber (formerly my MIT graduate student) and I have developed a new version of our 3-D hydrodynamic SPH (Smoothed Particle Hydrodynamics) code StarCrash incorporating a full treatment of PN gravity. We have used this code for a comprehensive new study of NS-NS binary coalescence in relativistic gravity [62] [73] [79]. We performed detailed comparisons with our previous Newtonian results [21, 27, 54] and we studied systematically the dependence of the gravitational-wave signals on the NS EOS, NS spins, and on the binary mass ratio.

Traditional PN expansions capture qualitatively many new relativistic effects (such as compactness effects and gravitational radiation reaction), but, quantitatively, they provide a rather poor approximation to the structure of a NS-NS binary.
Over the last few years, with Northwestern postdocs Josh Faber and Philippe Grandclément, we have developed the first fully relativistic version of our SPH code. Gravity is now treated using the equations of general relativity (rather than a PN expansion) in the so-called conformal flatness (CF) approximation, which suppresses some the dynamic degrees of freedom of the field. Our new CFSPH code [93] provides, for the first time, a completely consistent treatment of the initial conditions for close NS-NS binaries on quasi-circular orbits (always constructed in the CF approximation of GR) and the hydrodynamic merger calculation. A key result of our work is that NS EOS signatures may be present in the gravitational-wave signals at much lower frequencies than previously believed [81] [93].

Our new CFSPH code will also allow us for the first time to study realistically the hydrodynamic coalescence of NS-BH systems. These systems may be even more important than NS-NS binaries as sources of gravitational waves. Stable mass transfer from a NS to a BH, and the tidal disruption of a NS by a BH in a compact binary, have also been discussed recently by many authors in the context of GRBs. However, very little is known theoretically about these systems, especially when compared to the large amount of work done on NS-NS binaries. In the next few years of this project we plan to perform a systematic survey of the parameter space for NS-BH binaries, and, in particular, we will try to identify the boundary between stable mass transfer and rapid unstable coalescence. These hydrodynamic studies will also complement the many “Grand Challenge” efforts currently underway to try to simulate numerically the collision of two BH or the mergers of BH-BH binaries (which are purely gravitational problems in numerical relativity, without any fluid sources).

Northwestern postdocs Josh Faber and Philippe Grandclément, as well as several undergraduate (NSF REU) students have been working with me on this project, which has been funded by several grants from the NSF Gravitational Physics Program. REU students who have co-authored some of our papers include Justin Manor [73]. Movies of many of our 3D hydrodynamic calculations of binary mergers have been developed by REU student Rocky Jones [3d Astronomical Visualization].



* Numbers in square brackets refer to papers in the list of publications.

Blue stragglers are objects that appear as main-sequence (MS) stars above the turnoff point in the color-magnitude diagram of a globular cluster, i.e., they appear too massive for their environment. They are generally thought to be formed through collisions and mergers of lower-mass MS stars in the dense cluster core. Although this merger hypothesis has been with us for more than two decades, and is by now well accepted, detailed calculations of the merger process were lacking until a few years ago.

In 1995 I initiated a large collaborative effort to study theoretically the formation and evolution of blue stragglers resulting from stellar collisions, using a combination of 3D numerical hydrodynamics calculations with SPH [31, 34] and stellar evolution calculations. For the first time in a study of stellar collisions, the results of hydrodynamic calculations were combined with a state-of-the-art stellar evolution code (the YREC code developed at Yale by Demarque and collaborators). This allowed us, in particular, to study chemical mixing and stellar evolution of collision products in a completely self-consistent way [47]. This is crucial for obtaining accurate theoretical models of blue stragglers and quantitative predictions for their observable properties. The YREC code is used to construct realistic initial conditions (spherical stellar models) for our SPH collision calculations. In turn, the SPH code produces final merged configurations to be used as initial conditions for stellar evolution calculations of blue stragglers using YREC.

The initial phase of our work has focused on hydrodynamic mixing during the merger process as well as the subsequent evolution of merger remnants, first on a thermal timescale as they contract to the main sequence, and then on a stellar evolution timescale. We include a treatment of rotation and rotationally-induced mixing in the stellar evolution calculations (as provided by YREC). Our goal is to calculate accurately the entire evolution of blue stragglers from the moment of their formation in a dynamical merger, through their pre-main-sequence and main-sequence phases. This must be done for many different combinations of stellar masses and orbital parameters for the collisions. This initial phase of the project is now nearing completion [67, 68], and our results lead to detailed theoretical predictions for the observable characteristics of blue stragglers (colors, luminosities, numbers). Comparisons with observations (in particular, many large samples of blue stragglers detected with HST in several different clusters) can now be made to test dynamical models of star cluster evolution including stellar collisions, and to test the merger hypothesis for the origin of blue stragglers.

This ongoing collaboration includes J. Lombardi (Vassar) and A. Sills (McMaster), as well as Northwestern postdocs J. Fregeau and N. Ivanova.

A new era in astronomy began a few years ago with the first clear detections of several Jupiter-type planets around nearby solar-like stars. These discoveries of extrasolar planets will no doubt lead to significant improvements in our understanding of many processes related to planet formation, structure and evolution, as well as deeper questions such as the existence of extraterrestrial life in the Universe. NASA has identified the search for and characterization of planetary systems as one of its highest future priorities.

Beginning in 1996 I have been interested in exploring some of the theoretical implications of these recent discoveries. With Eric Ford (at the time a sophomore in Physics at MIT!) I wrote the first paper [44] proposing that the development of dynamical instabilities during the early stages of planet formation could provide a natural explanation for the very surprising orbital characteristics of the observed systems (which include giant planets in highly eccentric orbits, or in very tight circular orbits with periods as short as a few days). This idea is by now well accepted (and has received further observational support with the recent detection of three giant planets of comparable masses in a marginally stable configuration in the Upsilon Andromedae system). The long-term stability of the Solar System, in spite of its chaotic nature, may have been necessary for the development of intelligent life. However, it may also be very atypical, and may in fact require very special conditions during the early stages of planet formation. One of these special conditions may be the formation of a single, dominant giant planet (Jupiter in the Solar System). Our work, based on numerical integrations using simplectic methods, examines the dynamical consequences of the presence of several giant planets of comparable masses in the same system. As a first case, we have recently completed a systematic study of the dynamical evolution of unstable planetary systems containing two identical giant planets [69]. As a corollary, we have also studied the tidal interaction (including circularization of a highly eccentric orbit) between a giant planet and a solar-like star [40, 52].

Extrasolar planets have also been detected around radio pulsars. To this day, the only example we know of a system of earth-mass planets orbiting any star other than our own Sun is the system of planets around the millisecond radio pulsar PSR 1257+12, discovered by Wolszczan in 1992. Clear confirmation for at least two planets in this system came in 1994 when my theoretical prediction [8] for a near-resonant gravitational perturbation effect in the system was verified. More recent theoretical modeling of the pulsar timing data has revealed the presence of a more distant, giant planet in the system [45]. A giant planet has also been detected in orbit around PSR 1620-26, a binary millisecond pulsar located near the center of the globular cluster M4 [19, 26]. Among the many theoretical implications of this extraordinary system [56] is that planets must also be common in older stellar populations such as those found in globular clusters and the Galactic halo.