The Astrophysical Journal, 568:939-953, 2002 April 1
© 2002. The American Astronomical Society. All rights reserved. Printed in U.S.A.

Stellar Collisions and the Interior Structure of Blue Stragglers

James C. Lombardi, Jr., and Jessica S. Warren ¹
Department of Physics and Astronomy, Vassar College, 124 Raymond Avenue, Poughkeepsie, NY 12604-0745; lombardi@vassar.edu, jesawyer99@alum.vassar.edu
Frederic A. Rasio ²
Department of Physics, Massachusetts Institute of Technology 6-201, Cambridge, MA 02139; rasio@northwestern.edu
Alison Sills ³
Department of Physics and Astronomy, Leicester University, Leicester LE1 7RH, England, UK; asills@mcmail.cis.mcmaster.ca
and
Aaron R. Warren ¹
Department of Physics and Astronomy, Vassar College, 124 Raymond Avenue, Poughkeepsie, NY 12604-0745; aawarren00@alum.vassar.edu

Received 2001 June 29; accepted 2001 December 7

ABSTRACT

Collisions of main-sequence stars occur frequently in dense star clusters. In open and globular clusters, these collisions produce merger remnants that may be observed as blue stragglers. Detailed theoretical models of this process require lengthy hydrodynamic computations in three dimensions. However, a less computationally expensive approach, which we present here, is to approximate the merger process (including shock heating, hydrodynamic mixing, mass ejection, and angular momentum transfer) with simple algorithms based on conservation laws and a basic qualitative understanding of the hydrodynamics. These algorithms have been fine-tuned through comparisons with the results of our previous hydrodynamic simulations. We find that the thermodynamic and chemical composition profiles of our simple models agree very well with those from recent SPH (smoothed particle hydrodynamics) calculations of stellar collisions, and the subsequent stellar evolution of our simple models also matches closely that of the more accurate hydrodynamic models. Our algorithms have been implemented in an easy-to-use software package, which we are making publicly available.⁴ This software could be used in combination with realistic dynamical simulations of star clusters that must take into account stellar collisions.

     ¹ Present address: Department of Physics and Astronomy, Rutgers University, Piscataway, NJ 08854.
     ² Present address: Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208.
     ³ Present address: Department of Physics and Astronomy, McMaster University, Hamilton, ON L8S 4M1, Canada.
     ⁴ For the foreseeable future, this package can be downloaded from http://vassun.vassar.edu/~lombardi/mmas.

1. INTRODUCTION AND MOTIVATION

     Blue stragglers are stars that appear along an extension of the main sequence (MS), beyond the turnoff point in the color-magnitude diagram (CMD) of a cluster. All observations suggest that blue stragglers are indeed more massive than a turnoff star and are formed by the merger of two or more parent MS stars. In particular, Shara, Saffer, & Livio (1997) and Sepinsky et al. (2000) have directly measured the masses of several blue stragglers in the cores of 47 Tuc and NGC 6397 and confirmed that they are well above the MS turnoff mass, some even with masses apparently above twice the turnoff mass. Furthermore, Gilliland et al. (1998) have demonstrated that the masses estimated from the pulsation frequencies of four oscillating blue stragglers in 47 Tuc are consistent with their positions in the CMD.

     Stellar mergers can occur through either a direct collision or the coalescence of a binary system (Leonard 1989; Livio 1993; Stryker 1993; Bailyn 1995). Single-single star collisions occur with appreciable frequency only in the cores of the densest clusters (Hills & Day 1976), but in lower density clusters collisions can also happen indirectly, during resonant interactions involving binaries (Leonard 1989; Leonard & Fahlman 1991; Sigurdsson, Davies, & Bolte 1994; Sigurdsson & Phinney 1995; Davies & Benz 1995; Bacon, Sigurdsson, & Davies 1996). Merger rates depend directly on cluster properties, such as the local density, velocity dispersion, mass function, and binary fraction. When mergers do occur, all of these cluster properties are affected. The dynamics of a cluster, including mass segregation and the rate of core collapse, are consequently influenced, leading to an intricate relation between individual stellar interactions and cluster evolution (Hut et al. 1992; Rasio, Fregeau, & Joshi 2002). By studying stellar mergers, we are therefore able to probe the dynamics of globular clusters. Results from ongoing Hubble Space Telescope surveys of nearby globular clusters continue to expand the statistics of blue straggler populations, making it timely for a detailed comparison between observations and theory.

     The final structure and chemical composition profiles of merger remnants are of central importance, because they determine the subsequent observable properties and evolutionary tracks of merger products in a CMD (Sills & Bailyn 1999). Three-dimensional hydrodynamic simulation is one means by which we can focus on fluid mixing during stellar mergers and determine the structure of a remnant. Many such simulations of stellar mergers have been presented in the literature (Lombardi, Rasio, & Shapiro 1996; Sandquist, Bolte, & Hernquist 1997; Sills & Lombardi 1997; Sills et al. 2001). The problem with these simulations, if they were to be coupled with calculations of the cluster as a whole, is the prohibitive computing time: high-resolution hydrodynamic simulations can typically take hundreds or even thousands of hours to complete.

     In this paper, we develop a method for computing the structure and composition profiles of zero-age blue stragglers without running hydrodynamic simulations. Because our method takes considerably less than a minute to generate a model on a typical workstation, we are able to explore the results of collisions in a drastically shorter time. Our approach can be generalized to work for more than two parent stars, simply by colliding two stars first and then colliding the remnant with a third parent star. Most importantly, such algorithms will make it possible to incorporate the effects of collisions in simulations of globular clusters as a whole.

2. PROCEDURE

     We begin with two (nonrotating) parent-star models, specifying initial profiles for the stellar density , pressure P, and abundance of chemicals as a function of mass fraction. The profile for the entropic variable A P/, a quantity closely related (but not equal) to thermodynamic entropy, can also be calculated easily and is of central importance. Here is the adiabatic index of the gas. Because the quantity A depends directly upon the chemical composition and the entropy, it remains constant for each fluid element in the absence of shocks.

     Fluid elements with low values of A sink to the bottom of a gravitational potential well, and the A profile of a star in stable dynamical equilibrium increases radially outward. Indeed, it is straightforward to show that the condition dA/dr > 0 is equivalent to the usual Ledoux criterion for convective stability of a nonrotating star (Lombardi et al. 1996). The basic idea here can be seen by considering a small fluid element inside a nonrotating star in dynamical equilibrium. If this element is perturbed outward adiabatically (that is, with constant A), then it will sink back toward equilibrium only if its new density is higher than that of its new environment. If, instead, the element is less dense than its surroundings, it will continue to float away from the equilibrium, an unstable situation. Likewise, if an element is perturbed inward, its density must be less than the environment's in order to return toward equilibrium. Because pressure equilibrium between the element and its immediate environment is established nearly instantaneously, the ratio of densities satisfies _elem/_env = (A_elem/A_env)^-1/, by the definition of A. Therefore, if the perturbed element has a higher A than its new environment, it has a lower density, and buoyancy will push the element outward. Similarly, a fluid element with a lower A then its surroundings will sink. As a result, a stable stratification of fluid requires that the entropic variable A increase outward: dA/dr > 0. In such a star, a perturbed element will be subject to restoring forces that cause it to oscillate about its equilibrium position. For a detailed discussion of the stability conditions within rotating stars, see § 7.3 of Tassoul (1978) or Tassoul (2000). In practice, even in rapidly rotating stars, fluid distributes itself in such a way that the entropic variable A increases outward.

     During a collision, the entropic variable A of a fluid element can increase as a result of shock heating (see § 2.2). However, the relative impact speed of two MS stars in a globular cluster is comparable to the speed of sound in these parents: both of these speeds are of order (GM/R)^1/2, where G is Newton's gravitational constant and M and R are, respectively, the mass and length scales of a parent star. Consequently, the resulting shocks have Mach numbers of order unity and shock heating is relatively weak. Therefore, to a reasonable approximation, a fluid element maintains a constant A throughout a collision.

     The underlying principle behind our method exploits the two special properties of A that were just discussed. Namely, the entropic variable A will (1) increase outward in a stable star and (2) be nearly conserved during a collision. Therefore, to a good approximation, the distribution of fluid in a collisional remnant can be determined simply by sorting the fluid from both parent stars in order of increasing A: the lowest-A fluid from the parent stars is placed at the core of the remnant and is surrounded by shells with increasingly higher A. In this paper, we will further improve upon this approximation by also modeling mass loss, shock heating, fluid mixing, and the angular momentum distribution.

     Our algorithms are calibrated from the results of smoothed particle hydrodynamic (SPH) calculations presented in Lombardi et al. (1996) (for collisions of polytropic stars), as well as in Sills & Lombardi (1997) and Sills et al. (2001) (for collisions of more realistically modeled stars). For details and tests of our SPH code, see Lombardi et al. (1999). For reviews of SPH, see Monaghan (1992) or Rasio & Lombardi (1999). Characteristics of the various parent stars used in our calculations are summarized in Table 1, with thermodynamic profiles shown in Figure 1. The various collision scenarios we have considered are listed, along with mass-loss information, in Table 2.

Table 1   Parent-Star Characteristics

Fig. 1   Thermodynamic profiles of A, P/5/3, pressure P, and density as a function of enclosed mass m. Our three realistically modeled parent stars are displayed in the left panels, and the four polytropic models are displayed on the right. The dotted, short-dashed, long-dashed and solid curves refer to parent stars of 0.8, 0.6, 0.4, and 0.16 M, respectively. Logarithms are base 10 and units are cgs.

Table 2   Mass Loss

     The realistically modeled parent stars are based on calculations performed with the Yale Rotating Evolution Code (YREC), as discussed in Sills & Lombardi (1997). In particular, we evolved nonrotating MS stars with a primordial helium abundance Y = 0.25 and metallicity Z = 0.001 for 15 Gyr, the amount of time needed to exhaust the hydrogen in the center of a 0.8 M star. We note that P/^5/3 decreases slightly in the outermost layers of the 0.4 and 0.6 M stars modeled by YREC (see Fig. 1). The adiabatic index is actually less than the ideal gas value of 5/3 in these regions, because of the relative importance of ionization and radiation pressure. In this paper, however, we neglect these effects and instead simply force the A profile to increase by some negligibly small amount in these regions. Figure 2 displays chemical abundance profiles of these parent stars.

Fig. 2   Fractional chemical abundance (by mass) as a function of enclosed mass m for various chemical elements in our three realistically modeled parent stars. Line types are as in Fig. 1.

2.1. Mass Loss

     The velocity dispersion of globular cluster stars is typically only 10 km s^-1, which is much lower than the escape velocity from the surface of an MS star; for example, a star of mass M = 0.8 M and radius R = R has an escape velocity (2GM/R)^1/2 = 552 km s^-1. For this reason, collisional trajectories are well approximated as parabolic, and the mergers are relatively gentle: the mass lost is never more than about 8% of the total mass in the system (mass loss with hyperbolic trajectories is treated by Lai, Rasio, & Shapiro 1994). Furthermore, most MS stars in globular clusters are not rapidly rotating, and it is a good approximation to treat the initial parent stars as nonrotating.

     Given models for the parent stars (see Table 1), we first determine the mass lost during a collision. Inspection of hydrodynamic results for collisions between realistically modeled stars, as well as for collisions between polytropes, suggests that the fraction of mass ejected can be estimated approximately by



where c₁ and c₂ are dimensionless constants that we take to be c₁ = 0.157 and c₂ = 1.8, M₁M₂/(M₁ + M₂) is the reduced mass of the parent stars, R_{i, 0.5} and R_{i, 0.86} are the radii in parent star i enclosing a mass fraction m/M_i = 0.5 and 0.86, respectively, and r_p is the periastron separation for the initial parabolic orbit. While developing equation (1) we searched for a relation that accounted for the mass distribution (not just the total masses and radii) of the parent stars in some simple way. The more diffuse the outer layers of the parents, the longer the stellar cores are able to accelerate toward each other after the initial impact: the sum of half-mass radii, R_{1, 0.5} + R_{2, 0.5}, in the denominator of equation (1) accounts for this increased effective collision speed for parents whose mass distributions are more centrally concentrated. The dependence on in equation (1) arises from the expectation that the mass loss will be roughly proportional to the kinetic energy at impact and from the fact that a simple rescaling of the stellar masses (M_ikM_i) in a hydrodynamic simulation leaves f_L unchanged.

     The final (postshock) value of a fluid element's entropic variable will be higher than its initial value, as discussed in § 2.2. The mass loss must be distributed between the two parent stars in such a way that the outermost fluid layers retained from each parent have the same final entropic variable A, so that the layers can merge together into a stable equilibrium. We solve for this maximum value of A in the remnant by requiring that the mass of the fluid with higher final A be the desired ejecta mass. This constraint determines what fraction of the ejecta comes from each of the parent stars.

2.2. Shock Heating

     Shocks increase the value of a fluid element's entropic variable A = P/. The distribution and timing of shock heating during a collision involve numerous complicated processes: each impact generates a recoil shock at the interface between the stars, the oscillating merger remnant sends out waves of shock rings, and finally the outer layers of the remnant are shocked as gravitationally bound ejecta fall back to the remnant surface. For off-axis collisions this may be repeated several times. Our goal is not to derive approximations describing the shock heating during each of these stages, but rather empirically to determine physically reasonable relations that fit the available SPH data.

     Let A and A_init be, respectively, the final and initial values of the entropic variable for some particular fluid element. We used the results of hydrodynamic calculations to examine how the change A - A_init, as well as the ratio A/A_init, depended on a variety of functions of P_init (the initial pressure) and A_init. Our search for a simple means of modeling this dependence was guided by a handful of features evident from hydrodynamic simulations: (1) fluid deep within the parents is shielded from the brunt of the shocks; (2) in head-on collisions, fluid from the less massive parent undergoes less shock heating than fluid with the same initial pressure from the more massive parent; (3) in off-axis collisions with multiple periastron passages before merger, fluid from the less massive parent undergoes more shock heating than fluid with the same initial pressure from the more massive parent; and (4) the shock heating within each parent clearly must be the same if the two parent stars are identical. In all the hydrodynamic calculations considered, we model the fluid in our system using an adiabatic index = 5/3, corresponding to an ideal gas equation of state.

     We find that when log₁₀(A-A_init) is plotted versus log₁₀ P_init, the resulting curve for each parent star is fairly linear (see Fig. 3), with a slope of approximately c₃ = -1.1 throughout most of the remnant in the 25 simulations we examined:



Here the intercept b_i(r_p) is a function of the periastron separation r_p for the initial parabolic trajectory, as well as the masses M₁ and M₂ of the parent stars. Higher values of b_i correspond to larger amounts of shock heating in star i, where the index i = 1 for the more massive parent and 2 for the less massive parent (M₂ < M₁). For simplicity of notation, we have suppressed the index i on the A, A_init, and P_init in equation (2).

Fig. 3   Change in entropic variable A as a function of initial pressure Pinit on a log plot for the SPH remnant of cases e and g (see Table 2 for data describing these cases). Pentagons: Fluid from parent star 1 that has reached dynamical equilibrium by the end of the simulation and has been binned by enclosed-mass fraction; triangles: corresponding fluid from star 2. Logarithms are base 10 and units are cgs.

     The SPH data suggest that the intercepts b_i(r_p) can be fitted according to the relations







where c₄ = 0.5, c₅ = 5, c₆ = 2.5, and b₁(0) is the intercept for a head-on collision (r_p = 0) between the two parent stars under consideration.

     Although equations (2)(4) describe how to distribute the shock heating, the overall strength of the shock heating hinges on the value chosen for b₁(0). To determine b₁(0), we consider the head-on collision between the parent stars under consideration and exploit conservation of energy: more specifically, we choose the value of b₁(0) that ensures that the initial energy of the system equals the final energy during a head-on collision. Because we are considering parabolic collisions, the orbital energy is zero and the initial energy is simply E_tot = E₁ + E₂, the sum of the energies for each of the parent stars. The final energy of the system includes energy associated with the ejecta and the center-of-mass motion of the remnant, in addition to the energy E_r of the remnant in its own center-of-mass frame. In this paper we consider nonrotating parent stars, so the remnant of a head-on collision also is nonrotating, and its structure quickly approaches spherical symmetry. The values of E₁, E₂, and E_r are therefore simply the sum of the internal and self-gravitational energies calculated while integrating the equation of hydrostatic equilibrium. Because the energy E_r depends on the thermodynamic profiles of the remnant, it is therefore a function of a shock-heating parameter b₁(0) (see § 2.3 for the details of how the remnant's structure is determined).

     The energy E_r of the remnant is nearly equal to the initial energy E_tot of the system. However, the ejecta do carry away a portion of the total energy, suggesting that the energy conservation equation be written as



where the coefficient c₇ is of order unity and f_L is the fraction of mass lost during the collision (see § 2.1). We use a value of c₇ = 2.5, which is consistent with the available SPH data (see Table 3). In equation (5), the left-hand side is the initial energy of the system, and the right-hand side is its final energy. The second term on the right-hand side accounts for the energy associated with the ejecta and with any center-of-mass motion of the remnant (note that this term is positive because E_tot < 0). In practice, we iterate over b₁(0) until equation (5) is solved. It is necessary to solve equation (5) only once for each pair of parent-star masses M₁ and M₂: once b₁(0) is known, we can model shock heating in a collision with any periastron separation r_p by first calculating b₁(r_p) and b₂(r_p) from equations (3) and (4) and by then using these values in equation (2).

Table 3   Total Energy

2.3. Merging and Fluid Mixing

     As with any star in stable dynamical equilibrium, the remnant will have an A profile that increases outward. In our model, fluid elements with a particular postshock A-value in both parent stars will merge to become the fluid in the remnant with the same value of the entropic variable. Furthermore, if the fluid in the core of one parent star has a lower A-value than any of the fluid in the other parent star, the former's core must become the core of the remnant, because the latter cannot contribute at such low entropies. When merging the fluid in the two parent stars to form the remnant, we use the postshock entropic variable A, as determined from equation (2).

     Within the merger remnant, the mass m_r enclosed within a surface of constant A must equal the sum of the corresponding enclosed masses in the parents,



It immediately follows that the derivative of the mass in the remnant with respect to A equals the sum of the corresponding derivatives in the parents: dm_r/dA_r = dm₁/dA₁ + dm₂/dA₂, or dA_r/dm_r = [(dA₁/dm₁)^-1 + (dA₂/dm₂)^-1]^-1. In practice, we calculate these derivatives using simple finite differencing. If we partition the parent stars and merger remnant into mass shells, then two adjacent shells in the remnant that have enclosed masses that differ by m_r will have entropic variables that differ by



The value of A at a particular mass shell in the remnant is then determined by adding A_r to the value of A in the previous mass shell.

     In the case of the (nonrotating) remnants formed in head-on collisions, knowledge of the A profile is sufficient to determine uniquely the pressure P, density , and radius r profiles. While forcing the A profile to remain as determined from sorting the shocked fluid, we integrate numerically the equation of hydrostatic equilibrium with dm = 4r² dr to determine the and P profiles [which are related through = (A/P)^3/5]. This integration is an iterative process, because we must initially guess at the central pressure. Our boundary condition is that the pressure must be zero when the enclosed mass equals the desired remnant mass M_r = (1-f_L)(M₁ + M₂). During this numerical integration we also determine the remnant's total energy E_r and check that the virial theorem is satisfied to high accuracy. The total remnant energy E_r appears in equation (5), and if this equation is not satisfied to the desired level of accuracy, we adjust our value of b₁(0) accordingly and redo the shocking and merging process.

     As was done in Sills et al. (2001), the structure of a rotating remnant can be determined by integrating modified equations of equilibrium (see eq. [9] of Endal & Sofia 1976), once the entropic variable A and specific angular momentum j distribution are known (see § 2.4). To do so, one can implement an iterative procedure in which initial guesses at the central pressure and angular velocity are refined until a self-consistent model is converged upon. Even for the case of off-axis collisions and rotating remnants, the chemical composition profiles can still be determined, even without solving for the pressure and density profiles, as we will now discuss.

     Once the A profile of the remnant has been determined, we focus our attention on its chemical abundance profiles. Not all fluid with the same initial value of A_init is shock-heated by the same amount during a collision, because, for example, fluid on the leading edge of a parent star is typically heated more violently than fluid on the trailing edge. Consequently, fluid from a range of initial shells in the parents can contribute to a single shell in the remnant. To model this effect, we first mix each parent star by spreading out its chemicals over neighboring mass shells, using a Gaussian-like distribution that depends on the difference in enclosed mass between shells. Let X_i be the chemical mass fraction of some species X in a particular shell i, and let the superscripts "pre" and "post" indicate pre- and postmixing values, respectively. Then











where m_i is the mass of shell i, m_i is the mass enclosed by shell i, M is the total mass of the parent star, c₈ is a dimensionless coefficient that we choose to be c₈ = 5, and A_max and A_min are the maximum and minimum postshock entropic variables of fluid that will be gravitationally bound to the remnant. We have suppressed an additional index in equations (8)(10) that would label the parent star. The summand in equation (8) is the contribution from shell i to shell k. The second term in the distribution function, equation (9), is important only for mass shells near the center of the parent, while the third term becomes important only for mass shells near the surface; these two correction terms guarantee that an initially chemically homogeneous star remains chemically homogeneous during this mixing process (X = X = constant, for any shell k). The dependence of on A_max/A_min ensures that stars with steep entropy gradients are more difficult to mix (see Table 4 of Lombardi et al. 1996).

     Consider a fluid layer of mass dm_r in the merger remnant with a postshock entropic variable A that ranges from A_r to A_r + dA_r. The fraction of that fluid dm_i/dm_r that originated in parent star i can be calculated as (dA_r/dm_r)/(dA_i/dm_i). Therefore, the composition of this fluid element can be determined by the weighted average



where all derivatives are evaluated at A_r, the postshock value of the entropic variable under consideration. With the postshock A profiles given by equation (2) and the smoothed composition profiles given by equation (8), equation (11) allows us to merge the parent stars and determine the final composition profile of the remnant.

2.4. Angular Momentum Distribution

     To estimate the total angular momentum J_r of the remnant in its center-of-mass frame, we use angular momentum conservation in the same way that energy conservation was used in § 2.2. In particular, because the parent stars are nonrotating, the total angular momentum in the system is just the orbital angular momentum,



which we set equal to J_r plus a contribution due to mass loss (compare with eq. [5]):



The SPH simulations demonstrate that J_tot is always slightly larger than J_r, and the choice c₉ = 2 leads to good agreement with the SPH results. Equation (13) can be solved for J_r, and the results are compared with those of SPH simulations in Table 4. The agreement (between the numbers in the last two columns) is excellent for all cases with M₁/M₂ 2. For cases V and W, which have a relatively large mass ratio (M₁/M₂ = 5), the approximation begins to falter, although the predicted J_r-value still agrees with SPH results to within 10%.

Table 4   Total Angular Momentum

     The structure of the rotating remnants formed in off-axis collisions depends on the distribution of the specific angular momentum within the remnant. Even though the collisional remnants are axisymmetric around the rotation axis, with angular velocities that are constant on isodensity surfaces, the specific angular momentum distribution can nevertheless be quite complicated (see Fig. 12 of Lombardi et al. 1996 or Fig. 3 of Sills et al. 2001). The goal here is to simplify this complicated distribution into an average one-dimensional profile. The specific angular momentum j for the SPH remnants increases outward and is typically concave upward throughout most of the remnant when averaged over isodensity surfaces and plotted against enclosed mass.

     Once an approximate analytic form for the average j profile is specified, the profile can be constrained to satisfy



where m corresponds to the mass enclosed within a constant-density surface and J_r is determined from equation (13). We find that the relation



with c₁₀ = 0.6, is able to reproduce the important features of the specific angular momentum profile. Here, M_r is the remnant mass, c_s = (P/)^1/2 is the local sound speed, r is the local radius, and G is Newton's gravitational constant. Other forms for j(m) could also be used and normalized through equation (14). One advantage of equation (15) is that for m near 0 the specific angular momentum j(m) scales as m^2/3, in agreement with both simple analytic treatments and SPH results. Equation (15) is chosen because the rotation in the innermost regions (m < k₁M_r) of remnants is not strongly affected by r_p and because of the equation's loose resemblance to the j profile of a Keplerian disk for m close to M_r. The c_s and r profiles used in equation (15) are evaluated for a nonrotating equilibrium star with the same A profile as the star under consideration, a simplification that both eases and quickens the necessary computations. The coefficients k₁, k₂, and k₃ are not free parameters, but instead are determined by equation (14) and by the two additional constraints that j(m) and its derivative be continuous at m = k₁M_r. We always choose the smallest positive value of k₁ that meets these constraints. If no solution with k₁ > 0 exists, we set k₁ = k₃ = 0 and solve for k₂. For a nonrotating star, clearly k₁ = k₂ = k₃ = 0.

3. RESULTS

3.1. Comparison with Three-Dimensional Hydrodynamic Calculations

     To test further the accuracy of our simple models, we compare their structure and composition against models generated directly from the results of SPH calculations. These SPH calculations include both collisions between polytropic parent models (referred to with a capital letter as the case name) and collisions between realistically modeled parent stars (referred to with a lowercase letter).

3.1.1. Shock Heating

     Clearly, expressions for describing shock heating, such as equations (2)(4), are rather crude approximations that lump together complicated effects from the various stages of the fluid dynamics. However, these expressions do work quite well for parent stars of similar mass. To demonstrate this point, Figure 4 compares the shock heating of this method against the heating experienced by the individual SPH particles in six different collisions of identical parent stars, while Figure 5 presents similar data for four collisions between unequal-mass parents. The agreement between prediction and simulation is excellent for mass ratios M₁/M₂ from 1 to approximately 2, regardless of the periastron separation r_p. Even for a mass ratio as high as 5, the prescription continues to work well, at least for intermediate values of the periastron separation r_p (see case W in Fig. 5). For head-on collisions with high mass ratios, the predicted shock heating is an underestimate in the smaller star and in the center of the larger star (see case U in Fig. 5). It is worth noting that in collisions between two MS stars, a remnant must be relatively far past the turnoff on a CMD if it is to be identified as a blue straggler. It is therefore unlikely that collisions involving parent stars with a mass ratio more than 5 will produce a true blue straggler, simply because the remnant mass could not be considerably more than the turnoff mass.

Fig. 4   Change in the entropic variable A vs. the initial enclosed mass fraction m/M for collisions involving identical parent stars. Each point represents an individual SPH particle, and the solid curve shows the typical increase in A predicted by eq. (2). The dashed curve shows the log of the average change in A, where the averaging is done on the SPH data in 25 equally spaced bins along the m/M-axis. Logarithms are base 10 and units are cgs.

Fig. 5   Same as Fig. 4, but for collisions involving parent stars of unequal mass. For each case, the left and right panel displays data for the larger and the smaller star, respectively.

     The SPH data in Figures 4 and 5 clearly show that fluid from the same initial enclosed mass fraction m/M can be shock-heated by different amounts. This effect can be easily understood because, for example, fluid on the impact side of a star will be heated more than fluid with the same m/M on the back side. The treatment of mixing in § 2.3 does allow us to model this spread in shock heating, by redistributing fluid to positions with slightly higher or lower final A-values than what is given by equation (2). The extent of the redistribution is set by the smoothing parameter (see eq. [10]), taken to be constant over each parent star. Higher values of correspond to a smaller range of mass fractions over which the fluid is distributed. The parameter = 248, for example, for the parents in case a, while for case g, ₁ = 229 for the 0.8 M parent and ₂ = 143 for the 0.4 M parent. This approach does well at mimicking the overall effects of the spread in shock heating. However, mixing of fluid is somewhat overestimated in the core and underestimated in the outer layers, affecting predictions for the central concentration of helium (see § 3.1.3) and for the surface concentration of trace elements such as lithium (see § 3.1.2). Future treatments could perhaps improve results by implementing a position-dependent .

     Note that our shock-heating prescription does necessarily imply the desirable qualitative features that are discussed in § 2.2 and evident in the SPH data of Figures 4 and 5: (1) fluid with high initial pressure P_init (the fluid shielded by the outer layers of the star) is generally shock-heated less; (2) b₂(0) b₁(0), so that a less massive star is shock-heated less in head-on collisions; (3) b₁(r_p) decreases with r_p, while b₂(r_p) increases with r_p, so that for sufficiently large r_p we have b₂(r_p) > b₁(r_p), and the less massive star is shock-heated more; and (4) b₁(r_p) = b₂(r_p) whenever M₁ = M₂, so that identical parent stars always experience the same level of shock heating.

3.1.2. Mass Loss

     Although mass loss is very small in a parabolic collision, it is nevertheless important to understand well if one is interested in tracking trace elements that may only be found in outermost layers of the parent stars. The method of § 2.1 yields remnant masses that are seldom more than 0.01 M different from those yielded by a hydrodynamic simulation (see the last two columns of Table 2); this is clearly a significant improvement over neglecting mass loss completely, which sometimes can overestimate the remnant mass by more than 0.1 M. Table 2 also lists the fraction f₂ of the ejecta that originated in star 2 for each collision scenario, both as calculated from our simple method and as calculated from an SPH simulation. Our simple method allows the fraction of ejecta from each parent star to be determined simply by requiring that the outermost fluid retained from each parent has the same postshock entropy. Again, the agreement between prediction and simulation is excellent, especially for cases with mass ratios M₁/M₂ 2.

     The mass-loss and merging procedures implemented in this paper allow one to identify not only how much of the ejecta originated in each parent star, but also the original location of the fluid within each parent. The expression g_ikm_im_k/ g_jim_j gives the amount of mass in shell i that is transported to shell k (see eq. [8]). Therefore, the fraction of the mass from shell i that ultimately becomes ejecta is given by



where k_max corresponds to the shell with the highest-entropy fluid still gravitationally bound to the remnant. In equation (16), the sum over shells k includes only those shells with higher postshock entropies than this maximum, that is, only those shells associated with ejecta. Figure 6 displays the f_ejecta curves for the parent stars in a variety of collision scenarios, both as determined by this method and as determined by an SPH calculation.

Fig. 6   Local fraction of mass ejected by the collision, as a function of the initial enclosed mass m inside each parent star. The dotted curves represent the results of a three-dimensional SPH simulation, while the dashed curves represent the method of this paper. For cases a, k, and A, only a single curve is necessary for each line type, because the parent stars are identical and experience the same mass-loss distribution.

     Lithium is a particularly interesting element to follow during a collision. Lithium is burned during stellar evolution, except at low temperatures, and therefore can be used as an indicator of mixing. If a star has a deep enough surface convective layer, there will be essentially no lithium, because the convection mixes any lithium from the outer layers into the hot interior, where it is burned. A small amount of lithium does exist in the outer few percent of, for example, a 0.8 M turnoff star (see Fig. 2) and would consequently become part of the ejecta during a collision, resulting in a remnant with very little lithium.

3.1.3. Structure and Composition

     Thermodynamic (Fig. 7) and chemical composition (Figs. 8, 9, 10, and 11) profiles show that our remnant models are quite accurate. In case g, our remnant displays the kink in the A profile near m/M = 0.1 (see Fig. 7), inside of which the fluid originates solely from the 0.8 M star. Our models also reproduce the chemical profiles of the SPH remnant very well: the peak values in the chemical abundances are often accurate to within 20%, and the shapes of these profiles, although sometimes peculiar, are followed closely. Helium distribution is particularly important to model well because it will help determine the MS lifetime of the remnant. As mentioned in § 3.1.1, the core of the remnant is usually somewhat overmixed, flattening out the helium profile in that region. Nevertheless, the central value of the fractional helium abundance Y given by our models typically underestimates the SPH result by only about 5%.

Fig. 7   Thermodynamic profiles of A, pressure P, and density as a function of enclosed mass fraction m/Mr for the remnants of cases a and g, where Mr is the total bound mass of the remnant. The dotted line refers to the remnant resulting from a three-dimensional SPH simulation, and the dashed line refers to the remnant generated by the method of this paper. Logarithms are base 10 and units are cgs.

Fig. 8   Fractional chemical abundance (by mass) as a function of enclosed mass fraction m/Mr for the case a remnant. Line types are as in Fig. 7.

Fig. 9   Same as Fig. 8, but for the case g remnant

Fig. 10   Same as Fig. 8, but for the case e remnant

Fig. 11   Same as Fig. 8, but for the case k remnant

     Although near the remnant's surface our method sometimes yields a large fractional error in lithium abundance (see Figs. 8 and 9), this is simply because the overall abundance is so close to zero. For example, the predicted surface fractional ⁶Li abundance of 6.9 × 10^-10 for our case a remnant is an overestimate, but it is nearly 20 times smaller than the surface fractional abundance in the 0.8 M parent. Except in the extreme case of grazing collisions (when mass loss is exceedingly small), collisional blue stragglers should be severely depleted in lithium, a prediction that can be tested with appropriate observations (see Shetrone & Sandquist 2000 and Ryan et al. 2001).

3.1.4. Angular Momentum Distribution

     Our previous hydrodynamic simulations involving polytropic stars (Lombardi et al. 1996) implemented the classical form of the artificial viscosity (AV), which introduces a significant amount of spurious shear in our differentially rotating remnants. The effects of shear are discussed in § 4.2 of Lombardi et al. (1996) and studied in detail in Lombardi et al. (1999). Shear viscosity tends to weaken differential rotation, transporting angular momentum outward. This angular momentum transport acts on the viscous timescale, which is comparable to the total time of a typical simulation in our collisions between polytropic stars. We therefore avoid comparisons involving the angular momentum profiles of remnants from our polytropic SPH calculations.

     Our collisions between realistically modeled stars implemented the Balsara AV (Balsara 1995), with a viscous timescale that is significantly larger. In particular, Lombardi et al. (1999) show that the viscous timescale scales approximately with N for the classical AV (where N_N is the neighbor number) and N_N for Balsara AV. Because we used N_N = 64 in our polytropic simulations and N_N = 100 in our realistic simulations, the viscous timescale is larger in our collisions involving realistic parent stars by a factor of approximately 100/64^1/2 10. Consequently, in our SPH calculations using the Balsara AV, only a relatively small amount of specific angular momentum is spuriously transported out from the core.

     Specific angular momentum profiles, averaged over surfaces of constant density, are compared in Figure 12 for the three realistically modeled cases with rotating remnants: cases e, f, and k. We find excellent agreement between our simple models and their SPH counterparts. Our procedure for determining the angular momentum distribution, as described in § 2.4, yields values of k₁, k₂, and k₃ of, respectively, 0.161, 0.427, and 3.16 × 10¹⁶ cm² s^-1 for case e, 0.023, 0.570, and 5.93 × 10¹⁵ cm² s^-1 for case f, and 0.146, 0.415, and 4.56 × 10¹⁶ cm² s^-1 for case k. Note that the hump in the SPH j profile in the outer few percent of the remnant results from the need to terminate the simulation before all the gravitationally bound fluid has fallen back to the merger remnant: this artifact is gradually diminishing during the final stages of the SPH calculation.

Fig. 12   SPH specific angular momentum profiles (dotted curves), compared with the approximate profiles (dashed curves) generated from eq. (15). Logarithms are base 10 and units are cgs.

3.2. Stellar Evolution of Remnant Models

     A rigorous test of the validity of the simple models, which we have performed using YREC, is to compare their subsequent stellar evolution against that of SPHgenerated models. YREC evolves a star through a sequence of models of increasing age, solving the stellar evolution equations for interior profiles such as chemical composition, pressure, temperature, density, and luminosity. All relevant nuclear reactions (including p-p chains, the CNO cycle, triple- reactions, and light-element reactions) are treated. Recent opacity tables are used (ensuring that the remnant's position in a CMD can be accurately determined), and mixing mechanisms are incorporated. For blue stragglers, the various mixing processes can potentially carry fresh hydrogen fuel into the stellar core and thereby extend the MS lifetime of the remnant. Furthermore, any helium mixed into the outer layers affects the opacity and hence the remnant's position in a CMD. The free parameters in YREC (e.g., the mixing length) are set by calibrating a solar mass and solar metallicity model to the Sun.

     Using the method described in Sills et al. (1997), we used two of our simple models (cases a and g) as starting models in YREC and evolved the collision products from the end of the collision to the giant branch. Figure 13 shows the evolutionary tracks for these simple models (solid lines) and the SPHgenerated models (dotted lines). The tracks of the SPHgenerated models are discussed in Sills & Lombardi (1997).

Fig. 13   Evolutionary tracks for cases a and g. Solid lines: tracks of the simple models described in this paper; dotted lines: tracks for models generated directly from the output of SPH calculations. Symbols mark identifiable evolutionary phases in the evolution of the collision products: position of the star immediately following the collision (diamonds), zero age MS (circles), MS turnoff (squares), and the base of the giant branch (triangles).

     The agreement between the two sets of models is very good. Although there are some differences on the "preMS" portion of the tracks, this stage only lasts for approximately a Kelvin-Helmholtz timescale, a very small fraction of the total lifetime of the cluster. Indeed, for case a, the SPH and simple models reach the MS after 0.4 and 0.8 Myr, respectively, and the corresponding contraction times for the case g models are only 4 and 2 Myr. Because the stars contract to the MS so quickly, the exact preMS track may not be directly important for generating synthetic CMDs. Nevertheless, a reasonable model of a preMS star can be of interest: the radius of a newly born remnant, and hence its collisional cross section in a multistar interaction, is strongly dependent on how the fluid has been shocked, and, furthermore, the surface abundances are strongly dependent on how it has been mixed.

     Once the collision product reaches the MS, the two methods show very good agreement, and the subgiant and giant branch evolution of these stars are virtually identical. The MS lifetimes for the two different methods agree reasonably well. For case g, the SPH results give a lifetime of 850 Myr, while the simple models give 650 Myr. For case a, the SPH results give a lifetime of 80 Myr, compared to 180 Myr from the simple models. It should be noted that the case a remnant has central helium abundance near 100%, accounting for its short MS lifetime. For such remnants, even a slight inaccuracy in the core's helium profile has a large relative effect on how long the star remains on the MS. While the MS lifetime resulting from our method can be off by more than a factor of 2 for remnants with intact helium cores, the lifetime of such remnants is nevertheless a very small fraction of the lifetime of globular clusters, and therefore the simple models can still be useful for incorporating stellar collisions into dynamical models of globular cluster evolution.

4. CONCLUDING REMARKS

     An important question in the study of globular clusters is, what are the necessary features of stellar collisions that must be modeled in order to synthesize reliable theoretical CMDs? While detailed studies of the tracks and evolutionary timescales of rapidly rotating collision remnants will be necessary to answer this question fully, we do regard the features modeled in this paper (mass loss, shock heating, hydrodynamic mixing, and angular momentum distribution) as essential components to consider. Our treatment of ejecta allows for an accurate estimate of the total mass of the remnant, upon which the subsequent stellar evolution and MS lifetime depend sensitively. Furthermore, shock heating during a collision not only influences the structure of the remnant, but also helps determine if convective regions can develop during the future evolution: the entropy gradients implied by SPH results and by the shock-heating method of this paper tend to stabilize a contracting preMS star against convection. Previous stellar evolution studies of remnants have found that the tracks are sensitive to the assumptions made about how the fluid is mixed during a merger (Bailyn & Pinsonneault 1995; Sills & Bailyn 1999); our simple fluid-mixing algorithms give a compromise between previously used approximations that tend to bracket the actual amount of mixing during a collision. Finally, although the treatment of rotation in stellar evolution is a challenging problem, it is clear that rotational support and induced mixing must be considered, because they have profound consequences on stellar evolution (Sills, Pinsonneault, & Terndrup 2000; Sills et al. 2001); the form for the angular momentum distribution presented in this paper provides a simple means of generating very reasonable initial profiles for future studies of collisional remnants.

     The algorithms we have developed are implemented in a publicly available FORTRAN software package named "Make Me A Star."⁵ Researchers should be aware of the limitations of this method in terms of the structure, rotational properties, and evolutionary timescales of the models created, as outlined in this paper. This software does produce accurate models for a variety of collision scenarios, and we hope that it will be used in combination with realistic dynamical simulations of star clusters that must take into account stellar collisions.

     ⁵ For the foreseeable future, this package can be downloaded from http://vassun.vassar.edu/~lombardi/mmas.

     We thank J. Faber for his contributions to the SPH code, A. Thrall for testing our software package and for helpful comments, and the anonymous referee for valuable comments. J. C. L. acknowledges support from the Keck Northeast Astronomy Consortium, from a grant from the Research Corporation, and from NSF grant AST 00-71165. F. A. R. acknowledges support from NSF grants AST 96-18116 and PHY 00-70918, NASA ATP grant NAG 5-8460, and a Sloan Research Fellowship. This work was also partly supported by the National Computational Science Alliance under grant AST 98-0014N and utilized the NCSA SGI/Cray Origin2000 parallel supercomputer.

REFERENCES

• Bacon, D., Sigurdsson, S., & Davies, M. B. 1996, MNRAS, 281, 830 First citation in article | Crossref | ADS
• Balsara, D. S. 1995, J. Comput. Phys., 121, 357 First citation in article | Crossref | ADS
• Bailyn, C. D. 1995, ARA&A, 33, 133 First citation in article | Crossref | ADS
• Bailyn, C. D., & Pinsonneault, M. H. 1995, ApJ, 439, 705 First citation in article | Crossref | ADS
• Davies, M. B., & Benz, W. 1995, MNRAS, 276, 876 First citation in article | Crossref | ADS
• Endal, A. S., & Sofia, S. 1976, ApJ, 210, 184 First citation in article | Crossref | ADS
• Gilliland, R. L., Bono, G., Edmonds, P. D., Caputo, F., Cassisi, S., Petro, L. D., Saha, A., & Shara, M. M. 1998, ApJ, 507, 818 First citation in article | IOPscience | ADS
• Hills, J. G., & Day, C. A. 1976, Astrophys. Lett., 17, 87 First citation in article
• Hut, P., et al. 1992, PASP, 104, 981 First citation in article | IOPscience | ADS
• Lai, D., Rasio, F. A., & Shapiro, S. L. 1994, ApJ, 423, 344 First citation in article | Crossref | ADS
• Leonard, P. J. T. 1989, AJ, 98, 217 First citation in article | Crossref | ADS
• Leonard, P. J. T., & Fahlman, G. G. 1991, AJ, 102, 994 First citation in article | Crossref | ADS
• Livio, M. 1993, in ASP Conf. Ser. 53, Blue Stragglers, ed. R. A. Saffer (San Francisco: ASP), 3 First citation in article | ADS
• Lombardi, J. C., Jr., Rasio, F. A., & Shapiro, S. L. 1996, ApJ, 468, 797 First citation in article | Crossref | ADS
• Lombardi, J. C., Jr., Sills, A., Rasio, F. A., & Shapiro, S. L. 1999, J. Comput. Phys., 152, 687 First citation in article | Crossref | ADS
• Monaghan, J. J. 1992, ARA&A, 30, 543 First citation in article | Crossref | ADS
• Rasio, F. A., Fregeau, J. M., & Joshi, K. J. 2002, in The Influence of Binaries on Stellar Population Studies, ed. D. Vanbeveren (Dordrecht: Kluwer), in press First citation in article
• Rasio, F. A., & Lombardi, J. C., Jr. 1999, J. Comput. Appl. Math., 109, 213 First citation in article | Crossref
• Ryan, S. G., Beers, T. C., Kajino, T., & Rosolankova, K. 2001, ApJ, 547, 231 First citation in article | IOPscience | ADS
• Sandquist, E. L., Bolte, M., & Hernquist, L. 1997, ApJ, 477, 335 First citation in article | IOPscience | ADS
• Sepinsky, J. F., Saffer, R. A., Pilman, C. S., DeMarchi, G., & Paresce, F. 2000, AAS, 196, 41.06 First citation in article
• Shara, M. M., Saffer, R. A., & Livio, M. 1997, ApJ, 489, L59 First citation in article | IOPscience | ADS
• Shetrone, M. D., & Sandquist, E. L. 2000, AJ, 120, 1913 First citation in article | IOPscience | ADS
• Sigurdsson, S., Davies, M. B., & Bolte, M. 1994, ApJ, 431, L115 First citation in article | Crossref | ADS
• Sigurdsson, S., & Phinney, E. S. 1995, ApJS, 99, 609 First citation in article | Crossref | ADS
• Sills, A., & Bailyn, C. D. 1999, ApJ, 513, 428 First citation in article | IOPscience | ADS
• Sills, A., Faber, J. A., Lombardi, J. C., Jr., Rasio, F. A., & Warren, A. R. 2001, ApJ, 548, 323 First citation in article | IOPscience | ADS
• Sills, A., & Lombardi, J. C., Jr. 1997, ApJ, 484, L51 First citation in article | IOPscience | ADS
• Sills, A., Lombardi, J. C., Jr., Bailyn, C. D., Demarque, P. D., Rasio, F. A., & Shapiro, S. L. 1997, ApJ, 487, 290 First citation in article | IOPscience | ADS
• Sills, A., Pinsonneault, M. H., & Terndrup, D. M. 2000, ApJ, 534, 335 First citation in article | IOPscience | ADS
• Stryker, L. L. 1993, PASP, 105, 1081 First citation in article | IOPscience | ADS
• Tassoul, J.-L. 1978, Theory of Rotating Stars (Princeton: Princeton Univ. Press) First citation in article
• . 2000, Stellar Rotation (Cambridge: Cambridge Univ. Press) First citation in article | Crossref