STAR CLUSTER SIMULATIONS: THE STATE OF THE ART SVERRE J. AARSETH Institute of Astronomy, Madingley Road, Cambridge, England 9 9 9 Abstract. 1 Thispaperconcentratesonfourkeytoolsforperformingstarclustersimulations n developed during the last decade which are sufficient to handle all the relevant a dynamical aspects. First we discuss briefly the Hermite integration scheme which J is simple to use and highly efficient for advancing the single particles. The main 7 numericalchallengeisindealingwithweaklyandstronglyperturbedhardbinaries. A newtreatment of theclassical Kustaanheimo-Stiefel two-body regularization has 1 provedtobemoreaccurateforstudyingbinariesthanpreviousalgorithmsbasedon v divideddifferencesorHermiteintegration.ThisformulationemploysaTaylorseries 9 expansion combined with the Stumpff functions, still with one force evaluation per 6 step,whichgivesexactsolutionsforunperturbedmotionandisatleastcomparable 0 tothepolynomialmethodsforlargeperturbations.Stronginteractionsbetweenhard 1 0 binariesandsingle starsorotherbinariesarestudiedbychainregularization which 9 ensures a non-biased outcome for chaotic motions. A new semi-analytical stability 9 criterionforhierarchicalsystemshasbeenadoptedandthelong-termeffectsonthe / inner binary are now treated by averaging techniques for cases of interest. These h modificationsdescribeconsistentchangesoftheorbitalvariablesduetolargeKozai p cycles and tidal dissipation. The range of astrophysical processes which can now - o be considered by N-body simulations include tidal capture, circularization, mass r transferbyRoche-lobeoverflowaswell asphysicalcollisions, wherethemasses and t s radii of individual stars are modelled by syntheticstellar evolution. a : v Keywords: Numerical methods, KS-regularization, N-body problem i X r a 1. Introduction The study of self-gravitational N-body systems by direct integration poses many technical challenges which must be addressed. However, progressduringthelastdecadenowenablessuchproblemstobetackled with confidence. In this personal review of recent developments, we concentrate on four main numerical tools which appear to be sufficient forthetaskinhand.Thecorrespondingalgorithmsmaybesummarized under the following headings: • Hermite integration • Two-body regularization • Chain regularization • Hierarchical systems (cid:13)c 2008 Kluwer Academic Publishers. Printed in the Netherlands. 2 SverreJ.Aarseth These topics are discussed briefly in the subsequent sections, to- gether with an outline of current applications. Given adequate tools, a massive effort is still required in order to develop an efficient star cluster simulation code but these aspects are beyond the scope of the present contribution. 2. Hermite Scheme AlthoughtheHermiteintegrationschemewasdevelopedforthespecial- purposeHARPcomputer (Makino1991), itis alsoprovinghighly effec- tive for standard workstations as well as conventional supercomputers. Since coding is now somewhat simpler than for the traditional di- vided difference formulation (Ahmad and Cohen 1973, Aarseth 1985), it should be considered the method of choice for direct N-body simu- lations. It may also be remarked that Hermite integration is actually more accurate than divided differences for the same order. The main idea is again to employ a fourth-order force polynomial but now the two first terms are evaluated by explicit summation over all N particles, thereby enabling two corrector terms to be formed. At first sight it may seem rather expensive to extend the full summation to the force derivative since this also requires prediction of velocities. However, simplicity as well as increased accuracy combine to outweigh the drawback of extra operations, particularly if block-step predictions are introduced. We expand a Taylor series for the force F and its first derivative F(1) for each particle up to the third derivative about the reference time t as F = F +F(1)t+ 1F(2)t2+ 1F(3)t3, (1) 0 0 2 0 6 0 F(1) = F(1)+F(2)t+ 1F(3)t2. (2) 0 0 2 0 (2) SubstitutingF from(2)into(1)andsimplifyingthenyieldsthethird 0 derivative corrector 6 F(3) = (2(F −F)+(F(1)+F(1))t) . (3) 0 0 0 t3 Similarly, substituting (3) into (1) gives the second derivative corrector 2 F(2) = (−3(F −F)−(2F(1) +F(1))t) . (4) 0 0 0 t2 (1) Using F and F evaluated at the beginningof a time-step, the co- 0 0 ordinatesandvelocities arefirstpredictedtoorderF(1) forallparticles. Following determination of the new F and F(1) by summation over all StarClusterSimulations 3 the contributions, the two higher derivatives are obtained by (3) and (4). This gives rise to corrector terms for coordinates and velocities given by ∆r = 1 F(2)∆t4+ 1 F(3)∆t5, i 24 0 120 0 ∆v = 1F(2)∆t3+ 1 F(3)∆t4. (5) i 6 0 24 0 Given the high-order derivatives, individual time-steps can now be assigned in the usual way from some suitable convergence criterion. The overheads of predicting N coordinates and velocities at each time-step can be reduced considerably by adopting so-called hierarchi- cal time-steps (McMillan 1986), where the indicated values are trun- cated to be factor 2 commensurate. The apparent inefficiency of just a few particles sharing the same (small) step and yet requiring one full prediction is compensated by having a distribution of discrete levels (typically 16 for N ≃ 104) such that the number of predictions is sig- nificantly reduced with respect to the continuous case (say by factor of 100).Thisschemeisparticularlysuitableforthespecial-purposeHARP computers butlendsitself equally well to other architectures, including parallel supercomputers (Spurzem 1998). Somewhat surprisingly, the workstation code NBODY6 which is based on the Ahmad-Cohen (1973) neighbour scheme (Makino and Aarseth 1992, Aarseth 1994) is in fact slightlyfasterandmorestablethantheolderNBODY5codeforN = 1000 single particles and the same number of steps. 3. Two-Body Regularization The early 1970’s saw the introduction of the Kustaanheimo-Stiefel (1965) regularization for treating close encounters and hard binaries in N-body simulations (Bettis and Szebehely 1972, Aarseth 1972) and theelegant KSmethod has proved to bevery resilient. However, even a regularized two-body solution is subject to small but systematic errors whenstudiedoverlongtimes.Inordertoavoidthisundesirablefeature, the concept of energy stabilization has been tried for weak perturba- tions (Aarseth 1985). Although this procedureensures that the orbit is constrained to have the correct energy arising from the perturbation, the corresponding angular momentum is no longer conserved so well. The subsequent exploitation of adiabatic invariance (Mikkola and Aarseth 1996) by the so-called slow-down principle tends to alleviate this imperfection since now one KS orbit may represent a number of physical periods by augmenting the perturbation itself and neglecting short-period effects. As for the earlier claim that a time-symmetric KS 4 SverreJ.Aarseth methodwouldbesuperior(Funato etal.1996), itnowappearsthatthe requirement of variable time-steps for perturbed orbits cannot be ac- commodated (Kokubo et al. 1998). So far there is no evidence that the resulting eccentricities of cluster binaries studied by the stabilization schemecannotbetrusted,especiallybearinginmindthatthelong-term evolution of most binaries is predominantly subject to discrete changes ofarandomnature.Thecaseoflong-livedhierarchicalsystemsdeserves special consideration, however, but here additional effects should also be considered, as discussed in a subsequent section. AnalternativeKSregularizationschemehasbeenpresentedrecently (Mikkola and Aarseth 1998) which achieves a high accuracy without extracost.ThisnewapproachisbasedontheideaofatruncatedTaylor series,whereadditional correction termsrepresenttheneglected higher orders and which yields exact solutions in the unperturbed case. The new algorithm is again of Hermite type and will be outlined in the following. First, coordinates and velocities of the perturbers are predicted in the usual way (i.e. to first order), whereas the regularized coordinates and velocities (U,U′) are predicted to highest order. Here U(4),U(5) include the modified Stumpff (1962) functions ∞ (−z)k c˜ (z) = n! , (6) n (n+2k)! Xk=0 wheretheargumentisrelatedtothetime-stepbyz = −1h∆τ2 andhis 2 thespecificbindingenergy.Thesecoefficients onlydeviateslightlyfrom unity and a twelfth-order expansion (re-evaluated every step) appears sufficient. After transforming the physical coordinates and velocities to global values, the predictor cycle is completed by evaluating the perturbing acceleration F as well as its explicit derivative F˙. Because of the insufficient accuracy of the predicted deviation from unperturbed motion at the end of a step, the corrector cycle employs an iteration. Setting Ω = −1h, the basic equation of motion takes the 2 familiar form U(2) = −ΩU+ 1rLTF, (7) 2 whereL(U) is a 4×3 linear matrix andr = U·U is the separation. We express the new KS acceleration and its derivative (where F′ = rF˙) at the start of a step as (2) (2) U = −Ω U +f , (8) 0 0 0 0 U(3) = −Ω U′ +f(3), (9) 0 0 0 0 where f(2) = 1rQ, with Q = LTF, is the perturbed force contribution 0 2 evaluated after the previous predictor cycle. StarClusterSimulations 5 The two next Taylor series terms are constructed from the Hermite scheme.Usingthecurrentvalueofh(andΩ),predictedtofourthorder, we form the new perturbative functions at the end of the step f(2) = (Ω −Ω)U+ 1rQ, (10) 0 2 f(3) = (Ω −Ω)U′−Ω′U+ 1r′Q+ 1rQ′, (11) 0 2 2 (4) (5) from which the corrector derivatives f ,f are recovered by the 0 0 Hermite rule (Makino 1991). (4) (5) The expressions for U and U are readily formed in analogy 0 0 with Eqs. (8) and (9) which yield (4) (2) (4) U = −Ω U +f , (12) 0 0 0 0 (5) (3) (5) U = −Ω U +f . (13) 0 0 0 0 FromEqs.(8)-(13), theprovisionalsolution forU,U′ is thenobtained by the general expression (cf. Mikkola and Aarseth 1998), which con- tains the Stumpff functions. The treatment of the energy remains the same as for standard Hermite based on Ω′ = −U′·Q and the physical timeisobtainedfromintegratingt′ = U·Uwhichalsoinvolves Stumpff functions. Substituting for U(2), we write the second derivative as Ω(2) = Ω U·Q−f(2) ·Q−U′·Q′. (14) 0 Thetwo corrector terms constructed from Ω′ and Ω(2) are added to the predicted value without any Stumpff functions to yield an improved solution for Ω at the start of the next iteration or at the end point. Subsequent iterations repeat the procedure above, starting from Eq.(10)withoutre-evaluatingthephysicalperturbationanditsderiva- tive. Thus the new values of Eqs. (10) and (11) are based on the improved solution for U,U′ and r,r′, as well as the new Ω. In the present treatment, one iteration yields a significant improvement for modest perturbations and experience so far indicates that this may also be sufficient for strong interactions because of the shortening of the stepsize ∆τ (cf. Aarseth 1994). The corrector cycle ends by specifying new derivatives for use in the next prediction, as well as saving the perturbative derivatives (10) and (11) required for the Hermite scheme. This is completed by re- initializing Eqs. (8) and (9) at the end point, substituting f(2),f(3) as well as the iterated values of Ω′ and Ω(2). It is advantageous to employ the corrected values of r and r′ for this purpose; the re-evaluation of f(2) and f(3) is fast and also benefits the final quantities U(2) and U(3) to be used in the next prediction. A more accurate expression of the 6 SverreJ.Aarseth fourth KS derivative at the end of the interval is obtained by including the next order by U(4) = U(4) +U(5)∆τ, (15) 0 0 and similarly for the third derivative of the energy, Ω(3) = Ω(3) +Ω(4)∆τ. (16) 0 Theaboveschemehasbeenimplementedinthestateoftheartcodes NBODY4 and NBODY6 and has proved itself in large-scale simulations. Accuracy tests obtained by a toy code shows that high accuracy can be obtained with 30 steps per orbit for relatively weak perturbations, which is about half that required by the old stabilization scheme. A significant part of this gain is due to the modifications by the Stumpff coefficients, although the basic Taylor series (or Encke-type) formu- lation is also considerably more accurate than the standard method. The number of operations for a typical step is not much larger in the new method, including the overhead for the Stumpfffunctions and one iteration in the corrector. Hence the computational effort is less for typical calculations, although this depends on the actual number of perturbers. Finally, we remark that the Stumpff method also includes the slow-down scheme in exactly the same way as before. 4. Chain Regularization The concept of chain regularization is simple, yet the mathematical formulation is quite technical and this has acted as an impediment to wider usage. However, it enables new types of problems to be studied and is therefore worth the extra effort. The basic idea is a generaliza- tion of three-body regularization (Aarseth and Zare 1974) which treats two perturbed KS solutions with respect to a common reference body, where each two-body solution is described by regular equations. Thus an extension to four participating bodies merely introduces one more perturbed KS solution, although the formalism is somewhat different (Mikkola and Aarseth 1990). Once the step from three to four particles has been mastered, the general case becomes feasible (Mikkola and Aarseth 1993). The essential feature of chain regularization is that dominant in- teractions along the chain itself are treated as perturbed KS solutions and all the other attractions are included as perturbations. Hence it becomes imperative to select the chain vectors in such a manner as to minimizethe perturbations.Since weare dealing with dynamical inter- actions, the chain vectors need to be redrawn in response to changing StarClusterSimulations 7 configurations.Fortuitously,alltherelevantdecision-makingconstitute a minor overhead here since the integration is carried out by the high- order Bulirsch-Stoer (1966) scheme and a certain elasticity is tolerated as regards switching to more favourable chain vectors. Theequations of motion arederived fromaregularized Hamiltonian of the form Γ⋆ = g(H −E), (17) where H is expressed in terms of the coordinates and momenta and E is the internal system energy. Here the function g is given by the corresponding time transformation dt = gdτ (18) andchoosingtheinverseLagrangianenergy(L = T+Φ)ensuresregular solutions for any chain separation R . k The treatment begins by selecting a compact subsystem of three or fourparticles;i.e.so-calledB+S orB+B type.Externalperturbersare chosen in analogy with the KS implementation and the internal inte- gration includes any perturbation effect which also tends to change the total energy according to its separate equation of motion. At the same time, the c.m. motions are advanced by the standard Hermite scheme with due attention to the slightly modified form of the corresponding acceleration which requires a differential correction. The analogy with KS does not hold in one important respect since the chain membership may change before termination occurs. Thus an initial subsystem of four members may lose one member due to ejection, or an approaching perturber - a single particle or binary - may beadded.Alternatively, themembershipmay alsochange through physical collision. All the relevant corrections and re-initializations are performedin situ.Hencetheuseofchain variables isalsohighlybenefi- cial for the evaluation of nearly singular quantities. Chain termination usually occurs when a binary becomes well separated from one or two other members in which case the binary is accepted for KS treatment, whereastheremainingmembershipisinitializedbytheHermitescheme or even as a second KS system. The actual decision-making also takes into account the cluster environment and is therefore quite involved. Clustersimulationsofprimordialbinariesfrequentlyinvolveinterac- tions of two binaries where the size of one is much less than the other. In such cases even the powerful chain method becomes prohibitive because the shortest period is a small fraction of the local crossing time. Fortunately the principle of slow-down applied to weakly per- turbed KS binaries can also be employed here (Mikkola and Aarseth 1996). Thispermitsaconsistentstudyofbinarieswitharbitrarilyshort 8 SverreJ.Aarseth periodswhichwouldotherwisehaveto betreated asinertsystems.The implementation itself differs from the KS case since here we adjust the slow-down factor continuously according to the maximum apocentre perturbationexerted bytheotherchain members,ratherthanchoosing an appropriate discrete level (factor of 2) at each apocentre passage. Sincethestronginteractionsstudiedbythechainmethodareusually of short duration, the simulation code only allows one such case to be considered at a time for technical reasons. However, there is provision for studying one triple as well as one quadruple system by unperturbed three-body (Aarseth and Zare 1974) and chain (Mikkola and Aarseth 1990) regularization. Given a few hundred critical events in a typical clustersimulation,thelatterproceduresareusuallynotneededbutthis may change with the addition of more primordial binaries. 5. Hierarchical Systems The Solar neighbourhood contains many examples of multiple systems where the inner component of a binary is itself a binary, and levels of higher multiplicity also exist. Likewise, hard binaries in star clusters may acquire an outer component with sufficiently small eccentricity to be stable over many orbits. Hierarchical triples may be formed by the classical three-body capture mechanism in which the binary itself acts mainly as a point-mass. However, in clusters with significant binary populations such systems are more likely to form in strong interactions between two binaries since this involves two-body encounters. The sec- ond formation process was already identified in scattering experiments with colliding binaries which yielded a high percentage of positive out- comes (Mikkola 1983). Thus one way for such triples to become stable requires the impact parameter to exceed some critical value and yet be sufficiently small for the weakest binary to be disrupted, but other processes are also favoured, including exchange. Given a newly formed hierarchical triple, the question of long-term stability naturally arises. Depending on the period ratio, the direct calculation of a perturbed inner binary can be quite time-consuming even with KS regularization. However, since the corresponding semi- majoraxismaynotbesubjecttoanyseculareffectsitbecomespossible to adopt the centre-of-mass approximation and thereby only neglect cyclical changes oftheeccentricity. Various empiricalcriteriahavebeen obtained by fitting theresults of systematic three-bodycalculations for a restricted set of parameters (Harrington 1977, Eggleton and Kiseleva 1995). Based on these results, the so-called merger procedure has been employed for some time (Aarseth 1985). Thus provided the stability StarClusterSimulations 9 condition is satisfied, the inner binaryis replaced by its combined mass to facilitate KS treatment of the outer orbit. A more rigorous approach based on correspondence with the chaos boundary in the binary-tides problem (Mardling 1995) has yielded a semi-analyticalstabilitycriterionwhichholdsforquitelargemassratios and arbitrary outer eccentricities (Mardling and Aarseth 1998). Here the critical outer pericentre distance is given in terms of the inner semi-major axis, a , by in 2/5 (1+e ) Rout = C (1+q ) out a , (19) p (cid:20) out (1−e )1/2(cid:21) in out where q = m /(m + m ) is the outer mass ratio, e is the cor- out 3 1 2 out responding eccentricity and C ≃ 2.8 is determined empirically. This criterion is only valid for coplanar prograde motion and still ignores a weak dependence on the inner eccentricity. However, the general case of inclined orbits exhibit increased stability so that Eq. (19) represents an upper limit. Further tests suggests an inclination correction factor f = 1−0.3i/180 (with i in degrees) which has been adopted in practi- cal simulations; this is also in qualitative agreement with the original stability condition for retrograde orbits (Harrington 1972). Themerger treatment is only allowed while the pericentre condition is satisfied, after which the inner binary is re-initialized. A further refinement is included when the outer component itself is a binary. In the case of a B + B configuration, the smallest binary plays the role of the outer body in a triple. Since the correspond- ing chaos boundary is not very sensitive to a second extended object (Mardling 1991), we adopt an additional correction factor f = f + 1 0.1min(a /a ,a /a ), with a representing the second semi-major in 2 2 in 2 axis.Wealsomentionherethatevendoublehierarchiesmaybeformed, where a system of type B+S or B+B itself acquires an outer bound component. Such configurations do occur occasionally and procedures have therefore been developed for their special treatment. The criterion (19) above is concerned with long-term stability and hence the absence of escape. However, it is also of interest to consider thepossibilityofexchangebetweentheoutercomponentandonemem- beroftheinnerbinary.Accordingtoclassicaldevelopments (Zare1977, Szebehely and Zare 1977), the critical value for exchange in a coplanar prograde triple is given by G2f2(ρ)g(ρ) (c2E) = − , (20) crit 2(m +m +m ) 1 2 3 where c is the angular momentum and the functions f(ρ),g(ρ) are expressed in terms of the masses, with ρ determined by iteration from 10 SverreJ.Aarseth afifth-orderalgebraicequationforthecollinear equilibriumpoints.Nu- mericaltests showthatthechaos boundarygiven byEq.(19)lies above theexchange boundarywhen the masses are comparable and the latter only begins to overlap above q ≃ 5. Application of the exchange out criterion is therefore less useful in practical calculations. We also note that once an exchange occurs the final evolution will inevitably lead to escape. The long-term evolution of a hierarchical triple is characterized by cyclicoscillationsoftheinnereccentricity wheretheamplitudedepends on the inclination. The so-called Kozai effect (Kozai 1962) has received much attention recently in connection with external planetary systems butthereisalsoanearlyexamplefromN-bodysimulations(vanAlbada 1968) which points to the relevance for star clusters. Various analytical tools have been employed in order to model this process in some detail, includingtidaldissipationforhigheccentricities (MardlingandAarseth 1999). Among the useful quantities which can be calculated theoreti- cally (Heggie 1995) are the time-scale for a complete oscillation, T , K as well as its maximum value, e . max Since the time-scale for the Kozai cycle is usually much greater thantheKeplerperiod,themerger schemeforhierarchical triples lends itself particularly well to a semi-analytical treatment. At present only systems with e > 0.8 are considered since smaller amplitudes are max less likely to result in tidal activity. We have used a double averaging procedure (Eggleton 1997, Mardling and Eggleton 1998) to calculate the evolution of such systems in terms of the inner Runge-Lenz vector andangularmomentumvector. Thussomeexamplesshowthatinclina- tionsnear90◦ may inducetidalcircularization even ifoblateness effects are included. Clearly further developments of this experimental ap- proachisneededinordertoimprovethemodellingofthesecomplicated processes. 6. Astrophysical Applications The realistic simulation of star cluster dynamics requires a variety of astrophysical processes to be considered. In particular, the implemen- tation of consistent stellar evolution enables the study of mass loss and finite-size effects. This is an ongoing project which has been outlined elsewhere (Aarseth 1996) and now employs an improved description of Roche mass transfer and physical collisions (Tout et al. 1997). Partic- ular emphasis has been devoted to the modelling of chaotic motions and tidal circularization which form a link between an initial binary distribution and the Roche stage (Mardling and Aarseth 1999). In