We study the dynamical evolution of globular clusters using our 2D Monte Carlo code with the inclusion of primordial binary interactions for equal-mass stars. We use approximate analytical cross sections for energy generation from binary-binary and binary-single interactions. After a brief period of slight contraction or expansion of the core over the first few relaxation times, all clusters enter a much longer phase of stable "binary burning" lasting many tens of relaxation times. The structural parameters of our models during this phase match well those of most observed globular clusters. At the end of this phase, clusters that have survived tidal disruption undergo deep core collapse, followed by gravothermal oscillations. Our results clearly show that the presence of even a small fraction of binaries in a cluster is sufficient to support the core against collapse significantly beyond the normal core collapse time predicted without the presence of binaries. For tidally truncated systems, collapse is easily delayed sufficiently that the cluster will undergo complete tidal disruption before core collapse. As a first step toward the eventual goal of computing all interactions exactly using dynamical three- and four-body integration, we have incorporated an exact treatment of binary-single interactions in our code. We show that results using analytical cross sections are in good agreement with those using exact three-body integration, even for small binary fractions where binary-single interactions are energetically most important.

We present the results of a new study of mass segregation in two-component star clusters, based on a large number of numerical N-body simulations using our recently developed dynamical Monte Carlo code. Specifically, we follow the dynamical evolution of clusters containing stars with individual masses m_1 as well as a tracer population of objects with individual masses m_2=\mu m_1, using N=10^5 total stars. For heavy tracers, which could represent stellar remnants such as neutron stars or black holes in a globular cluster, we characterize in a variety of ways the tendency for these objects to concentrate in or near the cluster core. In agreement with simple theoretical arguments, we find that the characteristic time for this mass segregation process varies as 1/\mu. For models with very light tracers (\mu <~ 10^-2), which could represent free-floating planets or brown dwarfs, we find the expected depletion of light objects in the cluster core, but also sometimes a significant enhancement in the halo. Using these results we estimate the optical depth to gravitational microlensing by planetary mass objects or brown dwarfs in typical globular clusters. For some initial conditions, the optical depth in the halo due to very low-mass objects could be much greater than that of luminous stars. If we apply our results to M22, using the recent null detection of Sahu, Anderson, & King (2001), we find an upper limit of ~25% at the 63% confidence level for the current mass fraction of M22 in the form of very low-mass objects.

We summarize the results of recent theoretical work on the dynamical evolution of globular clusters containing primordial binaries. Even a very small initial binary fraction (e.g., 10%) can play a key role in supporting a cluster against gravothermal collapse for many relaxation times. Inelastic encounters between binaries and single stars or other binaries provide a very significant energy source for the cluster. These dynamical interactions also lead to the production of large numbers of exotic systems such as ultracompact X-ray binaries, recycled radio pulsars, double degenerate systems, and blue stragglers. Our work is based on a new parallel supercomputer code implementing Henon's Monte Carlo method for simulating the dynamical evolution of dense stellar systems in the Fokker-Planck approximation. This new code allows us to calculate very accurately the evolution of a cluster containing a realistic number of stars (N ~ 10^5 - 10^6) in typically a few hours to a few days of computing time. The discrete, star-by-star representation of the cluster in the simulation makes it possible to treat naturally a number of important processes, including single and binary star evolution, all dynamical interactions of single stars and binaries, and tidal interactions with the Galaxy.

We study the dynamical evolution of globular clusters using our new 2-D Monte Carlo code, and we calculate the lifetimes of clusters in the Galactic environment. We include the effects of a mass spectrum, mass loss in the Galactic tidal field, and stellar evolution. We consider initial King models containing N = 10^5 - 3x10^5 stars, and follow the evolution up to core collapse, or disruption, whichever occurs first. We find that the lifetimes of our models are significantly longer than those obtained using 1-D Fokker-Planck (F-P) methods. We also find that our results are in very good agreement with recent 2-D F-P calculations, for a wide range of initial conditions. Our results show that the direct mass loss due to stellar evolution can significantly accelerate the mass loss through the tidal boundary, causing most clusters with a low initial central concentration (Wo <~ 3) to disrupt quickly in the Galactic tidal field. Only clusters born with high initial central concentrations (Wo >~ 7) or steep initial mass functions are likely to survive to the present and undergo core collapse. We also study the orbital characteristics of escaping stars, and find that the velocity distribution of escaping stars in collapsing clusters looks significantly different from the distribution in disrupting clusters. We calculate the lifetime of a cluster on an eccentric orbit in the Galaxy, such that it fills its Roche lobe only at perigalacticon. We find that such an orbit can extend the lifetime by at most a factor of a few compared to a circular orbit in which the cluster fills its Roche lobe at all times.

We present the results of Monte Carlo simulations for the dynamical evolution of star clusters containing two stellar populations with individual masses m1 and m2 > m1, and total masses M1 and M2 < M1. We use both King and Plummer model initial conditions and we perform simulations for a wide range of individual and total mass ratios, m2/m1 and M2/M1. We ignore the effects of binaries, stellar evolution, and the galactic tidal field. The simulations use N = 10^5 stars and follow the evolution of the clusters until core collapse. We find that the departure from energy equipartition in the core follows approximately the theoretical predictions of Spitzer (1969) and Lightman & Fall (1978), and we suggest a more exact condition that is based on our results. We find good agreement with previous results obtained by other methods regarding several important features of the evolution, including the pre-collapse distribution of heavier stars, the time scale on which equipartition is approached, and the extent to which core collapse is accelerated by a small subpopulation of heavier stars. We briefly discuss the possible implications of our results for the dynamical evolution of primordial black holes and neutron stars in globular clusters.

We present a new parallel supercomputer implementation of the Monte-Carlo method for simulating the dynamical evolution of globular star clusters. Our method is based on a modified version of Henon's Monte-Carlo algorithm for solving the Fokker-Planck equation. Our code allows us to follow the evolution of a cluster containing up to 5x10^5 stars to core collapse in < 40 hours of computing time. In this paper we present the results of test calculations for clusters with equal-mass stars, starting from both Plummer and King model initial conditions. We consider isolated as well as tidally truncated clusters. Our results are compared to those obtained from approximate, self-similar analytic solutions, from direct numerical integrations of the Fokker-Planck equation, and from direct N-body integrations performed on a GRAPE-4 special-purpose computer with N=16384. In all cases we find excellent agreement with other methods, establishing our new code as a robust tool for the numerical study of globular cluster dynamics using a realistic number of stars.

We present a new theoretical analysis of the PSR B1620-26 triple system in the globular cluster M4, based on the latest radio pulsar timing data, which now include measurements of five time derivatives of the pulse frequency. These data allow us to determine the mass and orbital parameters of the second companion completely (up to the usual unknown orbital inclination angle $i_2$). The current best-fit parameters correspond to a second companion of planetary mass, $m_2 \sin i_2 \simeq 7\times10^{-3} M_\odot$, in an orbit of eccentricity $e_2\simeq 0.45$ and semimajor axis $a_2\simeq 60 $AU. Using numerical scattering experiments, we study a possible formation scenario for the triple system, which involves a dynamical exchange interaction between the binary pulsar and a primordial star-planet system. The current orbital parameters of the triple are consistent with such a dynamical origin, and suggest that the separation of the parent star-planet system was very large, $\go 50 $AU. We also examine the possible origin of the anomalously high eccentricity of the inner binary pulsar. While this eccentricity could have been induced during the same dynamical interaction that created the triple, we find that it could equally well arise from long-term secular perturbation effects in the triple, combining the general relativistic precession of the inner orbit with the Newtonian gravitational perturbation of the planet. The detection of a planet in this system may be taken as evidence that large numbers of extrasolar planetary systems, not unlike those discovered recently in the solar neighborhood, also exist in old star clusters.

We present a general method for determining the masses and orbital parameters of binary millisecond pulsars with long orbital periods (P_orb >> 1 yr), using timing data in the form of pulse frequency derivatives. We apply our method to analyze the properties of the second companion in the PSR B1620-26 triple system. We use the latest timing data for this system to constrain the mass and orbital parameters of the second companion. We find that all possible solutions have a mass m_2 in the range 2.4 10^-4 M_sun <= m_2 sin i_2 <= 1.2 10^-2 M_sun, i.e., almost certainly excluding a second companion of stellar mass and suggesting instead that the system contains a planet or a brown dwarf. Using Monte-Carlo realizations of the triple configuration in three dimensions we find the most probable value of m_2 to be 0.010(5) M_sun, corresponding to a distance of 38(6) AU from the center of mass of the inner binary (the errors indicate 80% confidence intervals). We also apply our method to analyze the planetary system around PSR B1257+12, where a distant, giant planet may be present in addition to the three well-established Earth-mass planets. We find that the simplest interpretation of the frequency derivatives implies the presence of a fourth planet with a mass of ~100 M_earth in a circular orbit of radius ~40 AU.

Previous timing data for PSR~B1620$-$26 were consistent with a second companion mass $m_2$ anywhere in the range $\sim10^{-3}-1\,M_\odot$, i.e., from a Jupiter-type planet to a star. We present the latest timing parameters for the system, including a significant change in the projected semi-major axis of the inner binary, a marginal detection of the fourth time derivative of the pulse frequency, and the pulsar proper motion (which is in agreement with published values for the proper motion of M4), and use them to further constrain the mass $m_2$ and the orbital parameters. Using the observed value of $\stackrel{\ldots.}{f}$, we obtain a one-parameter family of solutions, all with $m_2 \lo 10^{-2}\,M_\odot$, i.e., excluding stellar masses. Varying $\stackrel{\ldots.}{f}$ within its formal $1\sigma$ error bar does not affect the mass range significantly. However, if we vary $\stackrel{\ldots.}{f}$ within a $4\sigma$ error bar, we find that stellar-mass solutions are still possible. We also calculate the predicted rate of change of the projected semi-major axis of the inner binary and show that it agrees with the measured value.