Predictions for Triple Stars with and without a Pulsar in Star Clusters
aa r X i v : . [ a s t r o - ph ] A p r Mon. Not. R. Astron. Soc. , 1–11 (2005) Printed 30 October 2018 (MN L A TEX style file v2.2)
Predictions for Triple Stars with and without a Pulsar inStar Clusters
Michele Trenti ⋆ , Scott Ransom ⋆ , Piet Hut ⋆ and Douglas C. Heggie ⋆ Space Telescope Science Institute, Baltimore, MD 21210, U.S.A. National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA, 22903, U.S.A. Institute for Advanced Study, Princeton, NJ 08540, U.S.A. School of Mathematics and Maxwell Insitute of Mathematical Sciences, University of Edinburgh, King’s Buildings,Edinburgh EH9 3JZ, Scotland, U.K.
Accepted ; Received ; in original form
ABSTRACT
Though about 80 pulsar binaries have been detected in globular clusters so far, nopulsar has been found in a triple system in which all three objects are of comparablemass. Here we present predictions for the abundance of such triple systems, and for themost likely characteristics of these systems. Our predictions are based on an extensiveset of more than 500 direct simulations of star clusters with primordial binaries, anda number of additional runs containing primordial triples. Our simulations employ anumber N tot of equal mass stars from N tot = 512 to N tot = 19661 and a primordialbinary fraction from 0 − N = 19661 that include a mass spectrum with a turn-off mass at 0.8 M ⊙ ,appropriate to describe the old stellar populations of galactic globular clusters. Basedon our simulations, we expect that typical triple abundances in the core of a densecluster are two orders of magnitude lower than the binary abundances, which in itselfalready suggests that we don’t have to wait too long for the first comparable-masstriple with a pulsar to be detected. Key words: globular clusters – pulsars, general – methods: N-body simulations –stellar dynamics.
Numerical simulations of star clusters routinely pro-duce dynamically formed triple stars, especially in thepresence of primordial binaries (McMillan et al. 1990;Heggie & Aarseth 1992; McMillan & Hut 1994). Given thatmost globular clusters have a significant population of pri-mordial binaries, we would expect some pulsar binaries inglobular clusters to be actually part of a triple system. Insuch a case, most likely all stars have comparable masses,given the fact that encounters between binaries, and also be-tween binaries and single stars, are likely to expel the lighteststar, and hence subsequent encounters tend to concentratethe most massive stars in binaries and triples. Such a systemwe refer to as a pulsar triple , i.e. a triple system in whichall components have comparable mass and one is a pulsar.(This phrase is preferable to the alternative triple pulsar ,which would strongly suggest a triple system in which all ⋆ E-mail addresses: [email protected] (MT); [email protected](SR); [email protected] (PH); [email protected] (DCH) three components are pulsars.). The pulsar that is detectedthrough its radio emission can then either be one of the in-ner two stars, or it can be the third star in orbit around theother two. The possible set of configurations is quite large:each of the other two stars can be a neutron star, a massivewhite dwarf, possibly a main sequence star, and perhapseven a black hole. So far, no pulsar triple with compara-ble masses has been discovered, among the more than 120pulsars found in globular clusters.Ideally one would like to explore the formation of pulsartriples using detailed N-body simulations which include di-rect integration of the orbits and stellar evolution (thereforeusing one particle per star). Unfortunately this approach isnot computationally feasible today nor it will be in the nearfuture. In fact, even if we consider a pulsar rich globularcluster such as Terzan 5, the population synthesis models ofIvanova et al. (2007), which include a simplified treatmentof the dynamics of the system, predict a total of ≈ c (cid:13) Trenti et al. see Aarseth 2001; see also Sec. 6). This immediately high-lights that, in order to gather enough statistics, one wouldneed to run several tens of simulations of star clusters suchas Terzan 5, which has a present mass of 3 . · M ⊙ and anumber of stars N > · . This is outside the capabilitiesof current hardware, including the GRAPE6: the largest di-rect N-body simulations carried out to date for at least a re-laxation time are those by Baumgardt & Makino (2003) andhave N=131072 (see Hut & Trenti 2007). Thus one is forcedto introduce some simplications to address the formation oftriples stars with a pulsar. As a first step we use equal-masses calculations as discussed in Sec. 3. These providea reasonable approximation for the dynamics of old globu-lar clusters, where the turn-off mass is below 1 M ⊙ and thetypical mass ratio between two stars in the cluster is around2:1. Such an approach appears to be preferable to the useof simulations that start with very massive stars and up toseveral thousands particles, as is instead appropriate for thestudy of young open clusters (e.g. see de la Fuente Marcos1996; Hurley et al. 2005). In addition, we validate our equal-masses results against simulations with a mass spectrumthat has a turn-off mass at 0 . M ⊙ , consistent with the ageof stellar populations in galactic globular clusters.The aim of this paper is two-fold. First, we characterizethe typical fraction of triples that are formed dynamicallyin a star cluster with a sizeable population of primordialbinaries, comparing their number density with that mea-sured at late times in simulations that start with primordialtriples. Second, we take advantage of the fact that triple sys-tems with one or more pulsars have most likely a dynamicaland not primordial origin. We therefore apply our results tomake an estimate of the likely time we have to wait untilthe first pulsar triple will be found as well as to predict itsmost likely characteristics.The paper is organized as follows. In the next twosections we present an overview of observations and sim-ulations, respectively. In Sec. 4 we present the initialconditions for the simulations presented in this paper.Sec. 5 briefly summarizes the global evolution of a starcluster with primordial binaries, extensively discussed inHeggie, Trenti & Hut (2006); Trenti, Heggie & Hut (2007);Trenti et al. (2007). In Sec. 6 we develop a simple modelto describe the formation and destruction of triples. Weuse the model to interpret the observed triples fraction inour runs, attempting to characterize its N -dependence. InSec. 7 we discuss in detail the properties of these dynam-ically formed triples. Sec. 8 describes the results of someadditional runs where we have started with the presence ofprimordial triples. Sec. 9 presents the prospects for detectionof a triple system with a pulsar and Sec 10 sums up. There are currently ∼
130 pulsars known in globular clus-ters, the majority of which are in binary systems (for a re-cent review, see Camilo & Rasio 2005) . Since almost all ofthese systems are millisecond pulsars (MSPs) with extraor-dinary timing precisions (pulse arrival times are typically See P. Freire’s cluster pulsar catalog at ∼ pfreire/GCpsr.html measured to 10 − µ s), if they were part of a triple sys-tem, the other two components would be easily detected,even if one of the components was low-mass or in a long-period orbit (Thorsett et al. 1999).Since binary MSPs typically have orbital periods of sev-eral hours to several days, a pulsar triple system with theMSP in the internal orbit would likely be initially identi-fied as a normal binary pulsar system. Timing observationswould expose the third body in the external orbit after atime equal to a few percent of the external orbital period,which is likely to be less than a year.Alternatively, a pulsar triple system could have theMSP in a relatively long-period (months to years) externalorbit around a compact inner binary likely comprised of neu-tron stars, massive white dwarfs,or a combination of both,and possibly main-sequence stars. The MSP in such a sys-tem would be initially identified as a long-period binary orpossibly even an isolated pulsar. Timing observations wouldrelatively quickly reveal the MSP’s external orbit, but un-less strong orbital perturbations due to the inner binary arepresent, the internal orbit might never be detected and theMSP “companion” will be assumed to be a single massivestar or compact object instead of a binary.It is possible that stellar interactions could produce atriple system where two of the members are radio MSPs. TheMSP in the external orbit would be relatively easy to detectdue to the likely small orbital accelerations present over theshort time intervals (i.e. hours) used for pulsar searches. Foran MSP in the compact inner orbit, though, strong orbitalaccelerations would make the initial detection of the MSPvery difficult, likely requiring specialized algorithms and ex-tremely large amounts of computing (Ransom et al. 2003).Finally one could ask what is the a-priori chance to finda triple system with two of the three bodies being pulsar.We expect that this probability is significantly smaller thanthat of finding a triple system with just one pulsar. In factamong the ∼
100 binaries containing a pulsar, none is a dou-ble pulsar. This suggests that the probability for a triplewith two pulsars is at least two orders of magnitude smallerthan the probability of finding a single pulsar triple. Theexpected number of triples with two pulsars ( N T pp ) can beestimated as: N T pp = η N p N s N B ps N B N T , (1)where N p is the number of pulsars, N s is the number of singlestars, N B is the number of binaries, N B ps is the number ofbinaries with one pulsar, N T is the total number of triples ofthe system and η is the fraction of pulsars not bound to anycompanion in globular clusters ( η ≈ . N p ≈ N s ≈ , N B ps ≈ N B ≈ , N T ≈ we obtain N T pp ≈ − .Therefore it is unlikely that any double-pulsar triple will bedetected in globular clusters. The dense cores of globular clusters represent natural lab-oratories for studying exotic stellar populations, frequentlyliving in multiple systems, such as blue stragglers, low massX-ray binaries, cataclysmic variables, millisecond pulsars,stellar and, possibly, intermediate mass black holes. Most of c (cid:13) , 1–11 riple Stars with and without a Pulsar in Star Clusters these objects form through the combined effects of stellarevolution and stellar dynamics, especially due to interactionand evolution of binary stars, that may constitute up to 50%of the core mass of a globular cluster (e.g., see Albrow et al.2001; Bellazzini et al. 2002; Pulone et al. 2003). Observa-tional evidence for triple stars is also starting to accumulate(e.g., see Raboud & Mermilliod 1998; Sterzik & Tokovinin2002).From the theoretical point of view, a detailed model-ing of dense stellar systems is extremely challenging due tothe huge dynamical range involved. Even neglecting stellarevolution and hydrodynamics and focusing on gravitationalinteractions only, the local orbital timescale of a binary staris several orders of magnitude smaller than the global re-laxation timescale (hard binaries have an orbital period ofa few hours, while the half mass relaxation time may be upto a few billion years).Numerical simulations with primordial binaries havethus been performed either using approximate algorithmssuch as Fokker Planck or Monte Carlo methods (Gao et al.1991; Giersz & Spurzem 2000; Fregeau et al. 2003), thathave to rely on physical inputs on the interaction cross sec-tions (but see Fregeau & Rasio 2007 and Giersz & Spurzem2003 for Fokker Planck codes that treats binary interactionsvia direct N-Body integration), or using direct N-Body codesbut employing, until recently, only a limited number of parti-cles N ≈ (McMillan et al. 1990; Heggie & Aarseth 1992;McMillan & Hut 1994; de la Fuente Marcos 1996). Studiesof hierarchical systems in open clusters have been carried outby de la Fuente Marcos et al. (1997). A few direct numericalsimulations with N ≈ aimed at studying the formationand evolution evolution of hierarchical systems have beenrecently performed by Aarseth & Mardling (2000); Aarseth(2001) and by Aarseth (2000, 2004).To enhance the statistical base and remove some ofthe limitations of previous direct N-body investigations,we have started a project (Heggie, Trenti & Hut 2006;Trenti, Heggie & Hut 2007, hereafter HTH,THH respec-tively; see also Trenti et al. 2007) to study the evolutionof star clusters with primordial binaries by means of sev-eral hundred direct numerical simulations employing up to N tot = 19661 particles. Our aim is to provide a milestone forcomparison and validations of approximate methods and toinvestigate in detail the complex dynamical interactions ofmultiple stellar systems, that cannot so easily be character-ized in terms of approximate treatments like Fokker Planckor Monte Carlo methods.In this paper we take advantage of the extensive set ofnumerical simulations (more than 500) that we have per-formed (see HTH,THH) and we estimate the frequency andthe orbital properties of dynamically formed triple stars.We discuss their dependence on the dynamical state of thestar cluster (pre and post-collapse), on the initial binaryratio and on the number of particles N . We complementthese simulations (i) with previously unpublished runs thatinclude a mass spectrum appropriate to represent the colli-sional evolution of old stellar populations in globular clus-ters (that is with a turn-off mass below 1 M ⊙ ) and (ii) with afew equal mass runs with a significant population of primor-dial triples. Although most of our simulations are limitedto equal mass stars without stellar evolution, we briefly dis-cuss in the concluding section the applicability of our results to the prediction of the properties of globular cluster triplesystems of comparable mass that include one, or more, mil-lisecond pulsars. The general properties of the numerical simulations consid-ered in this paper to study the properties of triple systemshave been presented in detail in HTH and THH. To summa-rize, we use Aarseth’s NBODY-6 code (Aarseth 2003) con-sidering stars of equal mass and no stellar evolution. Theinitial distribution is either a Plummer model (HTH; in thiscase the system is considered isolated) or King models withconcentration parameter W = 3 , ,
11 (THH; here the ef-fects of a self-consistent galactic tidal field are also taken intoaccount). The number of objects N used ranges from 256 to16384; here N = N s + N b , where N b and N s are the initialnumber of binaries and single stars, respectively. The totalnumber of stars is N tot = N s + 2 N b . In our standard runs(see HTH and THH), the primordial binary population hasa fractional abundance ( f = N b / ( N s + N b )) from 0 to 50%,with an interal binding energy distribution flat in log scalein the range ≈ kT –680 kT , where (3 / kT is the meankinetic energy per particle of the system (the binaries beingreplaced by their barycenter). In addition we have performeda few simulations with N = 8192 and f = 20% that startfrom a King W = 7 model (with a self consistent tidal field)and adopt an extended binding energy range (from 5 kT to10 kT ) for primordial binaries. The aim of these runs isto specifically investigate the formation of triples originatedfrom neutron star binaries with a degenerate companion.In addition to these equal mass runs we also considertwo tidally limited runs ( W = 7 and W = 11 King models)with a mass spectrum. For these simulations, which have N = 16384 and f = 10%, the individual particle mass isdrawn from an initial mass function appropriate to studythe late evolutionary stages of star clusters, when stars are ≈
10 Gyr old. To build the initial mass function for therun we first consider a Miller & Scalo (1979) mass functionin the interval [0 . , M ⊙ , and then we proceed to takeinto account the effect of stellar evolution for massive starsby reducing the mass of particles above the turn-off (as-sumed at m toff = 0 . M ⊙ ) accordingly to the prescriptionof Hurley et al. (2000).All our results are presented using standard units(Heggie & Mathieu 1986) in which G = M = − E T ≡ G is the gravitational constant, M the total mass,and E T the total energy of the system of bound objects. Inother words, E T does not include the internal binding energyof the binaries, only the kinetic energy of their center-of-mass motion and the potential energy contribution whereeach binary is considered to be a point mass. The cor-responding unit of time is t d = GM / / ( − E T ) / ≡ t rh = 0 . Nr / h ln (0 . N ) . We recall that as our simulations consider gravitational in- c (cid:13) , 1–11 Trenti et al. teractions only, our results, expressed in terms of the dimen-sionless units described above, can be applied to any physi-cal choice for the total mass and scale radius of the system.The dependence of the triples’ properties on the number ofparticles employed is discussed in Sec. 6 and provides a wayto extrapolate the results of our simulations to a number ofparticles realistic for standard globular clusters.For reference in physical units, a globular cluster (de-scribed by a Plummer model) with N = 3 · stars, a totalmass of M = 3 · M ⊙ and a half-mass radius of 4 pc has ahalf-mass relaxation time t rh ≈ . · yr ; in this cluster abinary, formed by equal mass stars each of mass 1 M ⊙ , withbinding energy of 1 kT has a semi-major axis of ≈ AU and an orbital period p ∼
20 yrs.
The evolution of the system in our runs is driven by thebalance of two competing phenomena: the tendency of thecore to contract (undergoing a gravothermal collapse) andthe generation of energy due to binary-binary and binary-single interactions, that eventually halts the core contractionand fuels the half mass radius expansion. We can broadlyidentify two phases, extensively described in HTH and THH:an initial core radius adjustment and a quasi-steady binaryburning phase.The first “adjustment” transient is characterized by theevolution of the core toward a quasi-equilibrium configura-tion that will be maintained during the subsequent binaryburning phase. In this phase we may have, depending onthe initial conditions, a contraction (e.g., starting from aPlummer model, or from King models with W .
7) or an expansion , if the core is initially too small (Fregeau et al.2003 and the discussion of the W = 11 runs in THH).For equal mass stars this transient lasts up to ≈ t rh (0)and the expected core radius value at the end of this initialadjustment is well modeled in terms of general theoreticalconsiderations (Vesperini & Chernoff 1994). Essentially theefficiency of the production of energy due to binary burningdepends on the core density: if the core density is too low,the energy production is inefficient and the core shrinks; onthe other hand a very high core density leads to an excessiveenergy generation, with a consequent core expansion.After the “adjustment” transient, the details of the ini-tial conditions are largely erased from the system and a selfsimilar evolution sets in: the cluster expands by keeping theratio of the core to half mass radius almost constant. Thefuel for the expansion is provided by the hardening and de-struction of the primordial binary population, that can sus-tain this phase for about 100 t rh (0) in isolated models (thusfor a time much longer than the age of the universe for a typ-ical globular cluster) if the initial binary fraction is above10%.As a result of four body interactions, basically happen-ing in the core of the system only, where the interactionprobability is highest, a (small) fraction of stable triple starsis formed. In the next Section we discuss the properties ofthese multiple systems. Our runs start without primordial triples, but in a fewrelaxation times stable triples (identified as stable inNBODY6 using a semi-analytic approach based on the bi-nary tides problem - see Aarseth & Mardling 2000 andMardling & Aarseth 2001) are formed via dynamical inter-actions in binary-binary encounters following ejection of onestar. We recall that stable triples cannot be formed dueto three body encounters (see, for example, Heggie & Hut2003).After a first rapid rise of the number of triples thathappens in the first few relaxation times, a quasi stationaryregime sets in and the number of triples is approximatelyproportional to the number of binaries. This regime lastsfor about 30 relaxation times, a time longer than the ageof the universe for a typical globular cluster. Eventually, inisolated runs, as the evolution proceeds and the reservoir ofprimordial binaries is being depleted, the number of triplesalso drops. Fig 1 illustrates the evolution of the number oftriples in an isolated run (N=16384, 10% primordial bina-ries, starting from a Plummer model) up to ≈ t rh (0). Thedrop in N t /N for t > t rh (0) is clear, while the fraction oftriples to binaries declines more slowly. Fig. 2 shows insteadthe number of triples for a run that includes a galactic tidalfield, starting from a King W = 7 model with 20% pri-mordial binaries and N=16384. The tidal field destroys thecluster in about 30 t rh (0), thus even in the last stages of theevolution the binary fraction remains high (as single starspreferentially escape from the cluster) and so does the triplefraction. The asymptotic fraction of triple to single stars forour runs with equal-mass particles is given in Tab. 1. Thishas been obtained by averaging the triple to single stars ra-tio after the core contraption ( t > t cc ) and, in the case oftidal field runs, stopping when the number of stars in thesystem is reduced below 10% of its initial value.The evolution of the system is very similar even whena mass spectrum is introduced in the simulations (see Fig. 3and Tab. 2). Within the statistical uncertainties associatedto the small number of triples present at any time in thesystem ( N t . M ⊙ .The evolution of the number of triples can be under-stood in terms of the balance between the formation ( r f )and the destruction ( r d ) rate for triples. In first approxima-tion we can write locally: r f = α f · ρ b , (2) r d = ρ t ( α dB · ρ b + α dS · ρ s ) , (3)where ρ s , ρ b , ρ t are the density of singles, binaries and triplesrespectively, while the α i are numerical coefficients that de-pends on the cross section for the process considered.With the aid of Eqs. (2-3) we can explain the behaviorof Fig. 2. During the initial core contraction phase the for-mation rate r f is increasing as binaries accumulate in thecore, while the destruction rate r d is smaller due to the rel-ative absence of triples, whose numbers therefore increase. c (cid:13) , 1–11 riple Stars with and without a Pulsar in Star Clusters Figure 1.
Evolution of the number of triples with respect to thenumber of particles N (upper panel), to the number of binaries(middle panel) and number of binaries to the number of particles(lower panel) for a simulation with N = 16384 starting froma Plummer model (isolated, no tidal field) and 10% primordialbinaries. Due to the absence of a tidal field the mass loss timescaleis significantly longer than in the run presented in Fig. 2 and thedrop in the number of triples at later stages of the evolution isapparent. See HTH for more details on the run. Figure 2.
Evolution of the number of triples like in Fig. 1 but fora simulation with N = 16384 (with a tidal field) starting from aKing W = 7 model and 20% primordial binaries. The time is inunits of the initial half mass relaxation time t rh (0). After about30 t rh (0) the cluster is completely dissolved by the tidal field. SeeTHH for more details on the run. Figure 3.
Evolution of the number of triples like in Fig. 1 but fortwo simulations with N = 16384 (with a tidal field) starting froma King W = 7 model (red solid line) and W = 11 model (bluedashed line) with 10% primordial binaries. These simulations haveparticle masses drawn from a Miller & Scalo (1979) mass functionevolved to have a turn-off mass at 0 . M ⊙ . Around core collapse the formation rate reaches its peak andin the mean time an equilibrium regime sets in where thedensity of triples can be obtained by equating the formationand destruction rate ( r f = r d ): ρ t = α f · ρ b α dB · ρ b + α dS · ρ s . (4)We can further assume, based on geometrical arguments,that the cross section for a binary-triple encounter is greaterthan that of a triple-single encounter. In addition we recallthat after core collapse in simulations starting with f . ρ t = α f · ρ b α dB · ρ b = α f α dB ρ b . (5)This theoretical prediction is nicely confirmed in Fig. 4,where the abundance of triples after core collapse for runsstarting from N = 4096 is shown as a function of the frac-tion of primordial binaries: here the structural properties ofthe system are nearly identical for f &
10% (see Figs. 17 &18 in HTH) and there is a linear dependence of N t on f .Finally at later stages of the evolution of the system(i.e. well after the end of core-contraction phase), ρ b in thecore drops as binaries are being depleted, and the numberof triples also decreases.The fractional abundance of triples as a function of thenumber of particles used, for runs starting with f = 10%,either from a Plummer or from a King W = 7 model, isdepicted in Figs. 5-8. Especially for low N runs, before thecore collapse the triple abundance is lower than after it. c (cid:13) , 1–11 Trenti et al.
Figure 4.
Fractional abundance of triples averaged after the corecontraction for isolated runs starting from a Plummer model anddifferent primordial binary ratios. Each dot represents a differentrun, while the squares, with a 1 σ error bar, are the average valuesof the fractional abundances of triples associated with a givenprimordial binary ratio f . The number of triples appears to belinearly proportional to the primordial binary fraction. This is because the time for core collapse (in units of therelaxation time) increases with N (see HTH, Fig. 17), sothat when N is low enough the system has not yet reacheda balance between the formation and destruction rate oftriples. In the runs starting from a King model we observea marginally higher ratio of triples. This is due to the factthat the runs starting from a King model have a tidal field(see THH), so that the system steadily loses stars in itsoutskirts and the ratio of triples over total number of stars isenhanced in the post-collapse phase with respect to isolatedruns starting from a Plummer model.The N dependence of the number of triples at fixed ini-tial binary ratio is more difficult to characterize than thedependence on the binary ratio at fixed N . In fact, while inthe latter case the structural properties of the cluster arefixed, by varying N the core to half mass ratio is chang-ing (see Vesperini & Chernoff 1994, and HTH, THH), sothat the analysis in terms of Eqs. (2-3) is complicated bythe necessity to integrate over the density profile and bythe N-dependence of the α coefficients that describe the effi-ciency of the different formation/destruction channels. Em-pirically we can also note the intrinsic large scatter in themeasured triples abundances (see the dots representing in-dividual runs in Figs. 5-8), so that it is hard to draw afirm conclusion. Tentatively, the N dependence is logarith-mic at most, especially at lower N (in fact the core radiusdecreases approximately as 1 /log (0 . · N )). However for thelast three points (N=4096,8192,16384) the data are consis-tent with a constant value. In fact the decrease in the coreradius is compensated by a higher central density and byan increased number of binaries over singles in the core (see Figure 5.
Fractional abundance of triples averaged after the corecontraction for isolated runs starting from a Plummer model with10% primordial binaries and different numbers of particles. Sym-bols are as in Fig. 4
Figure 6.
Like Fig. 5 but for simulations with a tidal field startingfrom a King W = 7 model. Fig. 21 in HTH), so that the fractional abundance of triplesshould stay approximately constant.The evolution proceeds in a qualitatively similar wayin our runs initialized with an extended binary binding en-ergy range (up to 10 kT ). The only difference is that, allother conditions being equal, triples formation is enhancedby a factor between 1.5 and 2 when harder primordial bi-naries are present. In fact, binaries that have binding en- c (cid:13) , 1–11 riple Stars with and without a Pulsar in Star Clusters Figure 7.
Fractional abundance of triples averaged before thecore contraction for isolated runs starting from a Plummer modelwith 10% primordial binary and different number of particles.
Figure 8.
Like Fig. 7 but for simulations with a tidal field startingfrom a King W = 7 model. ergies above ≈ kT are dynamically inert and behaveessentially like singles during gravitational interactions withother stars. Therefore the encounter of one of these ultra-hard binaries with a regular (hard) binary can lead to anexchange encounter where one component of the regular bi-nary is replaced by the ultra-hard pair, much like during asingle-regular binary interaction. A more efficient produc-tion channel for stable triples is therefore available in thesesimulations. Therefore neutron star binaries with a degen- Table 1.
Asymptotic triple fraction in our equal mass runs. Thefirst column reports the number of particles N , the second theinitial density profile (Isolated Plummer model - P l or tidallylimited king model - W ), the binary fraction f is in the thirdcolumn. The fourth and fifth entries are the average asymptotictriple fraction and its standard deviation, computed over N runs realizations (last entry).N Model f h N t /N i σ N t /N N runs
512 Pl 0 0.80E-04 0.11E-03 48512 Pl 10 0.12E-02 0.68E-03 48512 Pl 20 0.22E-02 0.81E-03 47512 Pl 50 0.45E-02 0.14E-02 251024 Pl 0 0.35E-04 0.31E-04 491024 Pl 10 0.82E-03 0.39E-03 461024 Pl 20 0.19E-02 0.62E-03 491024 Pl 50 0.47E-02 0.13E-02 472048 Pl 0 0.13E-04 0.16E-04 462048 Pl 10 0.46E-03 0.16E-03 352048 Pl 20 0.17E-02 0.44E-03 392048 Pl 50 0.41E-02 0.93E-03 374096 Pl 0 0.50E-05 0.37E-05 74096 Pl 10 0.70E-03 0.27E-03 184096 Pl 20 0.18E-02 0.54E-03 114096 Pl 50 0.40E-02 0.62E-03 68192 Pl 0 0.18E-05 N/A 18192 Pl 10 0.63E-03 0.31E-04 28192 Pl 20 0.17E-02 N/A 116384 Pl 10 0.71E-03 N/A 1512 W = 7 10 0.16E-02 0.11E-02 161024 W = 7 10 0.14E-02 0.50E-03 152048 W = 7 10 0.87E-03 0.22E-03 154096 W = 7 10 0.65E-03 0.27E-03 68192 W = 7 10 0.76E-03 0.26E-03 416384 W = 7 10 0.81E-03 0.28E-04 216384 W = 7 20 0.15E-02 N/A 116384 W = 11 10 0.54E-03 N/A 116384 W = 11 20 0.22E-02 N/A 1 Table 2.
Asymptotic triple abundance for tidally limited runswith a Miller & Scalo mass spectrum. Columns as in Tab. 1.N Model f h N t /N i σ N t /N N runs W = 3 10 0.98E-03 N/A 116384 W = 7 10 0.68E-03 N/A 116384 W = 11 10 0.61E-03 N/A 1 erate companion and with short orbital periods have an en-hanced probability of ending up being in a triple systemcompared to main sequence binaries.To summarize, in a typical globular cluster ( N = 3 · )we expect an order of at least 100 triples if we conservativelyassume a binary fraction of 10% (with a standard bindingenergy range, up to contact main sequence binaries) and N t /N b ≈ · − (this accounts for the effect of a mass spec-trum and tidal field in setting N T /N b ; see Fig. 3). The num-ber approximately doubles by adopting an extended bindingenergy range that includes binding energies reached by neu-tron star binaries with a degenerate companion. c (cid:13) , 1–11 Trenti et al.
Figure 9.
Internal ( e int ) and external ( e ext ) eccentricity distri-bution for triples measured in a series of simulations with a tidalfield starting from King W = 3 , ,
11 models with 20% primor-dial binaries.
While for the majority of the standard simulations in oursample we recorded only the number of stable triples, wehave a small number of high resolution (N=16384) simula-tions with complete information on the orbital propertiesof triples, saved every 10 t d . In addition, all the extendedbinding energy range simulations have recorded informationof the orbital properties of multiple systems.For a sample of 3 simulations with N=16384 and f =20%, Fig. 9 shows the distribution of internal and exter-nal eccentricities. The inner eccentricity has a distributionpeaked at high eccentricities which derives (i) from the in-put thermal distribution of the eccentricity for primordialbinaries and (ii) from Kozai (1962) perturbations inducedby the outer body, which lead to growth of the inner binaryeccentricity (Aarseth 2001; Aarseth & Mardling 2000). Notethat our simulations are limited to Newtonian dynamics, sowe are missing relativistic corrections, which may alter theorbits of the inner binary (Aarseth 2007). The outer compo-nent of a stable triple is instead biased toward more circularorbits. This means that there is a selective disruption oftriple systems with eccentric outer orbits.For the same simulations the period distribution is re-produced in Fig. 10. Here we can note that the internalcomponent has an orbital period between two and three or-ders of magnitude shorter than the outer component. Theinternal orbit corresponds on average to a binary with bind-ing energy of a few hundreds of kT while the external orbithas a binding energy of a few kT . The distribution of orbitalperiods has a lower limit at log ( p/p ) = − .
25, where p isthe orbital period of a 1kT binary. This lower limit corre- Figure 10.
Internal (dotted) and external (solid) period distri-bution - in units of the period for a 1kT binary ( p ) - for triplesmeasured in a series of simulations with a tidal field startingfrom King W = 3 , ,
11 models with 20% primordial binariesand N = 16384. sponds to a binding energy of ≈ kT , i.e. a quasi-contactbinary with 1 M ⊙ main sequence components.The dynamical interactions in the star cluster lead toa fast destruction of triples with external binding energybelow 1 kT ; the destruction of the typical triples found in oursimulations (with external binding energy of 5 kT ) happenson a timescale of about 10 half mass relaxation times.The eccentricity distribution in our extended kT rangeruns is similar to that plotted in Fig. 9, showing only amarginally less marked circularization of the outer orbit.The distribution of the orbital periods is instead different(see Fig. 11). Stable triples are composed preferentially ofan inner ultra hard binary ( E b & kT ) and an externalcomponent that, on average, has a period only marginallyshorter than in the standard runs.The limited statistics in our simulations do not allow usto empirically constrain the N-dependence of the dynamicalproperties of triples. On the basis of theoretical modeling wedo not expect significant variations with N for single massstar clusters, so that the orbital properties discussed hereshould be representative for simulations of globular clusterswith a realistic number of particles. However one impor-tant caveat must be made as there is no guarantee thatthese properties are representative for the triple populationformed in a realistic model that includes a mass spectrumand stellar evolution. If a star cluster originates with a (small) population of pri-mordial triples, how does the number of surviving primordialtriples after several billions of years compare with the num- c (cid:13) , 1–11 riple Stars with and without a Pulsar in Star Clusters Figure 11.
Internal (dotted) and external (solid) period distribu-tion - like in Fig. 10 - for triples measured in a series of simulationswith a tidal field starting from King W = 7 models with 20%primordial binaries, N = 8192 and an extended initial bindingenergy range ([5 kT : 10 kT ]), representative for neutron starbinaries with a degenerate companion. ber of dynamically formed triples? To answer this questionwe have performed a few simulations with N = 4096 start-ing from an isolated Plummer model and with a number ofprimordial triples ( f trip = 5% and f trip = 10%, where f trip is the initial number fraction of triples). These runs are, tothe best of our knowledge, the first attempt to investigatethe dynamical evolution of a star cluster with a number ofparticles and of primordial triples greater than the few hun-dreds particles used by van den Berk et al. (2006).We started these runs by initializing the inner binarycomponents with the same procedure adopted in Heggie,Trenti & Hut (2006), i.e. in an energy range [5 : 680] kT . Anexternal component has then been added using the triplesinitialization subroutine within NBODY6 (Aarseth 2003),with an energy in the range [0 . kT when the externalenergy is defined as 2 m/a ext with m being the single particlemass and a ext being the external semi-axis. This means thatthe softest external component has a semi-axis 100 timeslarger than the softest inner binary. The evolution of thenumber of triples is shown in Fig. 12: after a slow start(due to the initial large core, and therefore low density oftriples in the core) primordial triples are steadily burned inthe first 15 t rh (0) with a rate approximately proportional tothe number of remaining triples (note that if multiplied bya factor 2 the curve for f trip = 5% reproduces closely thebehavior of the curve for f trip = 10%). This agrees well withthe expectation based on Eq. 3.At later times the destruction rate slows down slightlyin Fig. 12 due to the expansion of the half mass radius, whichincreases the instantaneous relaxation time with respect tothe initial one used in the figure. At about 25 t rh (0), a timelarger than the typical age of a galactic globular cluster, the Figure 12.
Evolution of the fraction of triples f trip in two runswith N = 4096 starting from an isolated Plummer model and f trip = 5% , number of primordial triples is higher than the number ofdynamically formed triples from 10 −
20% primordial bina-ries runs, if the initial fraction of triples is higher than aboutone percent. If we add also primordial binaries in runs withprimordial triples, the destruction rate of triples is slightlyenhanced, as expected from Eq. 3 (we have verified this witha N=1024 run with f trip = 3% and 7% primordial binaries). As described in Sec. 4, for typical cluster parameters and1 M ⊙ stars, a 1kT binary has a semi-major axis of ∼
10 AU,and an orbital period of p ∼
20 yrs. Such an orbit is on thehigh-end of the external orbit distribution as shown in Fig-ure 9. For an MSP in such a 1kT orbit, the instantaneousorbitally induced period derivative would be ˙ P orbit ∼ − (compared to a typical intrinsic ˙ P PSR ∼ − , and 1-yr tim-ing accuracies of ∆ ˙ P ∼ − ), and would be identified withina couple months of timing observations . Shorter orbital pe-riods would cause significantly larger period derivatives since˙ P orbit ∝ p − / , where p is the orbital period. There is no ev-idence for such orbitally induced accelerations in any of the ∼
60 globular cluster binary pulsars that have been or arecurrently being timed. This time is longer than what would be estimated by simplyconsidering the six orders of magnitude difference in the accu-racy of the 1-yr timing vs. the required accuracy for detection ofa triple system. This is because there are many degeneracies (forexample a positional offset) that can mask as a period deriva-tive and that require a few months of observations before beingexcluded.c (cid:13) , 1–11 Trenti et al.
Over the past several years, intensive globular clusterpulsar searches have uncovered numerous systems that arealmost certainly due to exchange interactions. These sys-tems include at least 10 recycled pulsars in highly eccen-tric ( e > .
25) orbits with periods between ∼ −
30 days(e.g. Freire et al. 2004; Ransom et al. 2005). Such systemsare very similar to the external orbits found in the triplesdescribed in Figure 9. However, precise timing observa-tions rule out the existence of internal orbits composed ofat least one main-sequence star to a high degree of con-fidence. Surveys have also uncovered another class of ex-change products, the 3 pulsar − “main-sequence” binariesNGC6397 A, Terzan 5 P, and Terzan 5 ad (D’Amico et al.2001; Ransom et al. 2005; Hessels et al. 2006). These sys-tems have circular orbits of duration 0.3 − − R ⊙ ) and mass-function constraints indicatelikely main-sequence dwarf companions.The simulations we present suggest there should be ap-proximately 100 times more binaries than triples in a typicalglobular cluster (see Figure 2). Given that we know of ∼
10 DISCUSSION
In this paper we have presented a systematic investigation ofthe frequency of dynamically formed triples for an extensiveset of direct simulations starting with a significant fractionof primordial binaries.We have shown that on a timescale of a few relaxationtimes the abundance of triples reaches an equilibrium valuethat depends linearly on the primordial density of binariesand that is about two orders of magnitude smaller. Sim-ulations including a tidal field present a relatively largerfraction of triples ( N t /N ) with respect to isolated runs withsimilar initial conditions. The presence of a mass spectrumtends instead to reduce the triple fraction.Stable triple stars in our standard runs, representativefor main sequence stars, are primarily formed due to fourbody encounters that lead to the escape of a single star.The formation mechanism influences the properties of thesesystems, which are primarily made of a hard inner binary,and which in a typical globular cluster would have an aver-age orbital period of a few days, accompanied by an externalstar with a period on average between 100 and 1000 timeslonger. The inner component is also, on average, on a slightlymore eccentric orbit than the outer member of the triple. Infact eccentric outer orbits are preferentially perturbed bygravitational encounters or break up the inner binary.Neutron stars binaries with a degenerate companion canreach binding energy well above the 700 kT limit of two main sequence stars. To quantify the rate of formation of stabletriples involving two degenerate stars we have carried out aseries of simulations with binding energy distributed up to10 kT . In these runs we observe an enhanced triple abun-dance, as an ultra-hard binary made of two degenerate starsessentially behaves like a single body during gravitationalencounters and can efficiently interact with a standard bi-nary to form a triple via an exchange-like encounter. In atypical globular cluster, the inner orbital period may in thiscase be of the order of a few hours, while the outer compo-nent of the system is expected to be of the order of a fewyears, much like in the case of triples with main sequencestars.Triple stars with at least one pulsar component are ex-pected to have a dynamical origin, as it is highly unlikelythat a pre-existing triple system can avoid being disruptedby the Supernova explosion that forms the pulsar. Thereforewe expect that the properties of pulsar triples are in agree-ment with those derived by runs starting with primordialbinaries only, where all the triples have been formed troughstellar encounters (though a single pulsar could become amember of a primordial triple via an exchange process). Todate, ∼
100 pulsar binaries are known in the globular clustersystem (Camilo & Rasio 2005). Based on our idealized sim-ulations, we expect that the discovery of a pulsar in a triplesystem of comparable masses is likely to happen soon. How-ever, we note that the majority of the known pulsar binarieshave companions of mass ∼
10% of the pulsar mass ratherthan ∼ M ⊙ . How this difference affects our conclusions iscurrently unknown and will demand more detailed studieswith a spectrum of stellar masses beyond that used in thisstudy.Higher order hierarchical systems, such as quadrupletsand quintuplets are also occasionally observed during ourruns, and their average number density is about two ordersof magnitude smaller than that of triples. This makes amore detailed characterization of their properties extremelychallenging, at least for simulations with only a limitednumber of particles like ours. A survey aimed at quantifyingthe observed fraction of multiple stars in globular clusterswould allow us to understand if the observed frequencyof triples and higher order systems is consistent withdynamical production from primordial binaries. Our inves-tigation suggests that a triple to binary ratio up to ≈
2% isconsistent with a purely dynamical origin. If the observednumber of triples is higher, then primordial triples have tobe introduced in the theoretical modeling, as investigatedfor small N open clusters ( N = 182) by van den Berk et al.(2006). Acknowledgements
We thank the referee for a careful reading of themanuscript and John Fregeau for useful suggestions. MTwas supported in part by NASA award HST-AR.10982.
REFERENCES
Aarseth S., 2000, in The Formation of Hierarchical Systemsin Star Clusters Schielicke, R. E. ed., AGM, 16, T06AAarseth, S. & Mardling R. A. 2000, astro-ph/0011514 c (cid:13) , 1–11 riple Stars with and without a Pulsar in Star Clusters Aarseth S., 2001, in Dynamics of Star Clusters and theMilky Way, ASP Conference Series, 228, 111Aarseth S., 2003, Gravitational N-body Simulations. Cam-bridge University PressAarseth S. 2004, in “The environment and evolution of dou-ble and multiple stars, ed. Allen & Scarfe, IAU Colloquim191.Aarseth, S. 2007, MNRAS, 378, 285Albrow M. D., Gilliland R. L., Brown T. M., EdmondsP. D., Guhathakurta P., Sarajedini A., 2001, Ap.J., 559,1060Baumgardt, H. and Makino, J. 2003, MNRAS, 340, 227Bellazzini M., Fusi Pecci F., Messineo M., Monaco L., RoodR. T., 2002, A.J., 123, 1509Camilo, F. and Rasio, F.A. 2005, in ASP Conf. Ser. 328:Binary Radio Pulsars, ed. F.A. Rasio & I.H. Stairs, 147D’Amico, N. and Possenti, A. and Manchester, R. N. andSarkissian, J. and Lyne, A. G. and Camilo, F. 2001, Ap.J.,561, L89de la Fuente Marcos, R. 1996, A&A, 314, 453de La Fuente Marcos, R. and Aarseth, S. J. and Kiseleva,L. G. and Eggleton, P. P. 1997, Astrophysics and SpaceScience Library, 223, 165van den Berk, J., Portegies Zwart, S. F. and McMillan,S. L. W. 2006, MNRAS, submitted, astro-ph/0607456Freire, P. (2005), in Binary Radio Pulsars, ASP ConferenceSeries, F. A. Rasio and I. H. Stairs eds. Vol. 328, 405Freire, P. C. and Gupta, Y. and Ransom, S. M. andIshwara-Chandra, C. H. 2004, Ap.J., 606, L53Fregeau J. M., G¨urkan M. A., Joshi K. J., Rasio F. A.,2003, Ap.J., 593, 772Fregeau, J. M. and Rasio F. A. 2007, Ap.J., 658, 1047Fregeau J. M., G¨urkan M. A., Rasio F. A., 2005, to ap-pear in Few-Body Problem: Theory and Computer Sim-ulations, ed C. Flynn, Annales Universitatis Turkuensis;also astro-ph/0512032Gao B., Goodman J., Cohn H., Murphy B., 1991, Ap.J.,370, 567Giersz M., Spurzem R., 2000, MNRAS, 317, 581Giersz M., Spurzem R., 2003, MNRAS, 343, 781Heggie D. C., Aarseth S. J., 1992, MNRAS, 257, 513Heggie D. C., Mathieu R. D., 1986, in LNP 267: The Use ofSupercomputers in Stellar Dynamics Standardised Unitsand Time Scales. p 233Heggie, D.C & Hut, P. (2003) “The Gravitational Million-Body Problem: A Multidisciplinary Approach to StarCluster Dynamics” 2003, Cambridge University PressHeggie, D.C., Trenti, M. & Hut, P. (2006), MNRAS, 368,677Hessels, J. W. T. and Ransom, S. M. and Stairs, I. H. andFreire, P. C. C. and Kaspi, V. M. and Camilo, F. 2006,Science, 311, 1901Hurley, J. R. and Pols, O. R. and Tout, C. A. 2000, MN-RAS, 315, 543Hurley, J. R. and Pols, O. R. and Aarseth, S. J.and Tout,C. A. 2005, MNRAS, 363, 293Hut, P. & Trenti, M. 2007, Scholarpedia Encyclopedia ofAstrophysics Kozaim Y. 1962, A.J., 67, 591 Ivanova, N. and Heinke, C. and Rasio, F. A. and Belczyn-ski, K. and Fregeau, J. 2007, MNRAS, submitted, astro-ph/0706.4096Lorimer, D. R., 2005, Living Rev. Relativity, 8, 7 (Novem-ber 9, 2005 version)Mardling, R. A. & Aarseth, S. 2001, MNRAS, 321, 3McMillan S., Hut P., 1994, Ap.J., 427, 793McMillan S., Hut P., Makino J., 1990, Ap.J., 362, 522Mikkola, S. 1983, MNRAS, 203, 1107Mikkola, S. 1984, MNRAS, 207, 115Miller, G. E. and Scalo J. E. 1979, Ap.J.S., 41, 513Pulone L., De Marchi G., Covino S., Paresce F., 2003,A&A, 399, 121Raboud, R. & Mermilliod, J. C. 1998, A&A, 329, 101Ransom, S. M. and Hessels, J. W. T. and Stairs, I. H. andFreire, P. C. C. and Camilo, F. and Kaspi, V. M. andKaplan, D. L., 2005, Science, 307, 892Ransom, S. M., Cordes, J. M., & Eikenberry, S. S. 2003,Ap.J., 589, 911Spitzer, L. 1987, Dynamical Evolution of Globular Clus-ters. Princeton, NJ, Princeton University Press, 1987Sterzik, M. F. & Tokovinin, A. A. 2002, A&A, 384, 1030Thorsett, S. E. and Arzoumanian, Z. and Camilo, F. andLyne, A. G. 1999, Ap.J., 523, 763Trenti, M., Heggie, D.C. & Hut, P. (2007), MNRAS, 374,344Trenti, M., Ardi, E., Mineshige, S. & Hut, P. (2007), MN-RAS 374, 857Vesperini, E., & Chernoff, D. F. 1994, Ap.J., 431, 231This paper has been typeset from a TEX/ L A TEX file preparedby the author. c (cid:13)000