Monte Carlo simulations of star clusters – II. Tidally limited, multimass systems with stellar evolution

A revision of Stodółkiewicz's Monte Carlo code is used to simulate evolution of large star clusters. The new method treats each superstar as a single star and follows the evolution and motion of all individual stellar objects. A survey of the evolution of N -body systems influenced by the tidal field of a parent galaxy and by stellar evolution is presented. The process of energy generation is realized by means of appropriately modified versions of Spitzer's and Mikkola's formulae for the interaction cross-section between binaries and field stars and binaries themselves. The results presented are in good agreement with theoretical expectations and the results of other methods (Fokker–Planck, Monte Carlo and N -body). The initial rapid mass loss, resulting from stellar evolution of the most massive stars, causes expansion of the whole cluster and eventually leads to the disruption of less bound systems ( W ₀=3). Models with larger W ₀ survive this phase of evolution and then undergo core collapse and subsequent post-collapse expansion, like isolated models. The expansion phase is eventually reversed when tidal limitation becomes important. The results presented are the first major step in the direction of simulating evolution of real globular clusters by means of the Monte Carlo method.     

Monte Carlo simulations of star clusters — I. First Results

A revision of Stodoíkiewicz's Monte Carlo code is used to simulate evolution of star clusters. The new method treats each superstar as a single star and follows the evolution and motion of all individual stellar objects. The first calculations for isolated, equal-mass N -body systems with three-body energy generation according to Spitzer's formulae show good agreement with direct N -body calculations for N  = 2000, 4096 and 10 000 particles. The density, velocity, mass distributions, energy generation, number of binaries, etc., follow the N -body results. Only the number of escapers is slightly too high compared with N -body results, and there is no level-off anisotropy for advanced post-collapse evolution of Monte Carlo models as is seen in N -body simulations for N  ≤ 2000. For simulations with N  > 10 000 gravothermal oscillations are clearly visible. The calculations of N   2000, 4096, 10 000, 32 000 and 100 000 models take about 2, 6, 20, 130 and 2500 h, respectively. The Monte Carlo code is at least 10⁵ times faster than the N -body one for N  = 32 768 with special-purpose hardware. Thus it becomes possible to run several different models to improve statistical quality of the data and run individual models with N as large as 100 000. The Monte Carlo scheme can be regarded as a method which lies in the middle between direct N -body and Fokker–Planck models and combines most advantages of both methods.     

Orbital evolution of long-period comets I. Trivariate Monte Carlo simulation

The mean orbital evolution of long-period comets for 16 representative initial orbits to short-period comets is calculated by a Monte Carlo method. First, trivariate perturbation distributions of barycentric Kepler energy, total angular momentum, and its z component in single encounters of comets with Jupiter are obtained numerically. Their characteristics are examined in detail and the distributions are found to be simple, symmetric, and easy to handle. Second, utilizing these distributions, we have done trivariate Monte Carlo simulations of the orbital evolution of long-period comets, with special emphasis on high-inclination orbits. About half of the 16 initial orbits are traced up to 5000 returns. For each of these orbits, the mean values of semimajor axis, perihelion distance, and inclination; their standard deviations, survival, and capture rates; as well as time scales of orbital evolution are calculated as functions of return number. Survival rates of the initial orbits with high inclination (～90°) and small perihelion distance (～1–2 AU) have been found to be only two or three times smaller than those of the main-source orbits of short-period comets established quantitatively by Everhart. The time scales of orbitsl evolution of the former, however, are nearly 10 times longer than the latter. There is a general trend that, for smaller perihelion distance, the survival efficiency becomes higher. The results of this paper should be considered a basis for a succeeding paper (Paper II) in which the physical lifetime of comets will be determined, and a comparison with the orbital data will be done.     

Long-term evolution of Oort Cloud comets: capture of comets

We test different possibilities for the origin of short-period comets captured from the Oort Cloud. We use an efficient Monte Carlo simulation method that takes into account non-gravitational forces, Galactic perturbations, observational selection effects, physical evolution and tidal splittings of comets. We confirm previous results and conclude that the Jupiter family comets cannot originate in the spherically distributed Oort Cloud, since there is no physically possible model of how these comets can be captured from the Oort Cloud flux and produce the observed inclination and Tisserand constant distributions. The extended model of the Oort Cloud predicted by the planetesimal theory consisting of a non-randomly distributed inner core and a classical Oort Cloud also cannot explain the observed distributions of Jupiter family comets. The number of comets captured from the outer region of the Solar system are too high compared with the observations if the inclination distribution of Jupiter family comets is matched with the observed distribution. It is very likely that the Halley-type comets are captured mainly from the classical Oort Cloud, since the distributions in inclination and Tisserand value can be fitted to the observed distributions with very high confidence. Also the expected number of comets is in agreement with the observations when physical evolution of the comets is included. However, the solution is not unique, and other more complicated models can also explain the observed properties of Halley-type comets. The existence of Jupiter family comets can be explained only if they are captured from the extended disc of comets with semimajor axes of the comets   a <5000 au  . The original flattened distribution of comets is conserved as the cometary orbits evolve from the outer Solar system era to the observed region.     

Possible Contribution of Mature γ-ray Pulsars to Cosmic-ray Positrons

We restudy the possible contribution of mature gamma-ray pulsars to cosmic ray positrons based on the new version of outer gap model. In this model, the inclination angle and average properties of the outer gap are taken into account, and more mature pulsars can have the outer gap and emit high energy photons. Half of the primary particles in the outer gaps will flow back toward the star surface and emit synchrotron photons, which can produce electron/positron pairs by the cascade of pair production. Some of these pairs will escape from the light cylinder and be accelerated to relativistic energies in the pulsar wind driven by low-frequency electromagnetic waves. Using a Monte Carlo method, we obtain a sample of mature gamma-ray pulsars and then calculate the production of the positrons from these pulsars. The observed excess of cosmic positrons can be well explained by this model.     

Multidimensional simulations of radiative transfer in Type Ia supernovae

A three-dimensional Monte Carlo code for modelling radiation transport in Type Ia supernovae is described. In addition to tracking Monte Carlo quanta to follow the emission, scattering and deposition of radiative energy, a scheme involving volume-based Monte Carlo estimators is used to allow properties of the emergent radiation field to be extracted for specific viewing angles in a multidimensional structure. This eliminates the need to compute spectra or light curves by angular binning of emergent quanta. The code is applied to two test problems to illustrate consequences of multidimensional structure on the modelling of light curves. First, elliptical models are used to quantify how large-scale asphericity can introduce angular dependence to light curves. Secondly, a model which incorporates complex structural inhomogeneity, as predicted by modern explosion models, is used to investigate how such structure may affect light-curve properties.     

Planetary impact probabilities for long-period comets

Planetary impact probabilities for long-period (near-parabolic) comets are determined by averaging Öpik's equations over inclination and perihelion distance for each planet. These averaged values compare well with the results of more elaborate Monte Carlo calculations. The impact probabilities are proportional to the square of the normalized capture radius of each planet, which in turn is a function of the planet's radius and mass, so that the major planets have the highest impact probabilities. Encounter velocities have an average value of $\text{3}\text{1}\text{2}$ times the planetary orbital velocity but the most probable encounter velocities are slightly higher than this for the terrestrial planets and slightly lower for the major planets. Comparison of the impact probabilities with the cratering record, corrected for gravity and velocity effects, indicates that long-period comets may account for 3 to 9% of the observed large crattes (diameter > 10 km) on the terrestrial planets. The inclination and perihelion properties of the impact probabilities obtained from numerical averaging provide a simple method for determining the impact probabilities for nonuniform distributions. The perihelion distribution of long period comets from J. A. Fernandez ((1981) Astron. Astrophys.96, 26–35) results in a crater production rate quite similar throughout the solar system, unlike that of a uniform perihelion distribution.     

Single star scattering on an open cluster

We examine the scattering of single stars from an open star cluster. The probability of the capture of a star by a star cluster is dependent on the velocity and mass of the star, and the stars that are not captured experience a velocity change. For low-velocity stars there is an exponential decrease of the capture probability with the initial velocity, and the velocity change decreases almost linearly. For high-velocity stars there is a v ⁻⁶ dependence for the capture probability, and a v ⁻¹ dependence for the velocity change. Analytical estimations, Monte Carlo and full N -body simulations are all in good agreement.     

Monte Carlo simulations of star clusters – IV. Calibration of the Monte Carlo code and comparison with observations for the open cluster M67

We outline the steps needed in order to incorporate the evolution of single and binary stars into a particular Monte Carlo code for the dynamical evolution of a star cluster. We calibrate the results against N -body simulations, and present models for the evolution of the old open cluster M67 (which has been studied thoroughly in the literature with N -body techniques). The calibration is done by choosing appropriate free code parameters. We describe, in particular, the evolution of the binary, white dwarf and blue straggler populations, though not all channels for blue straggler formation are represented yet in our simulations. Calibrated Monte Carlo runs show good agreement with results of N -body simulations not only for global cluster parameters, but also for, for example, binary fraction, luminosity function and surface brightness. Comparison of Monte Carlo simulations with observational data for M67 shows that it is possible to get reasonably good agreement between them. Unfortunately, because of the large statistical fluctuations of the numerical data and uncertainties in the observational data the inferred conclusions about the cluster initial conditions are not firm.     

Monte Carlo simulations of star clusters – III. A million-body star cluster

A revision of Stodółkiewicz's Monte Carlo code is used to simulate the evolution of million-body star clusters. The new method treats each superstar as a single star and follows the evolution and motion of all individual stellar objects. The evolution of N -body systems influenced by the tidal field of a parent galaxy and by stellar evolution is presented. All models consist of 1 000 000 stars. The process of energy generation is realized by means of appropriately modified versions of Spitzer's and Mikkola's formulae for the interaction cross-section between binaries and field stars and binaries themselves. The results presented are in good agreement with theoretical expectations and the results of other methods. During the evolution, the initial mass function (IMF) changes significantly. The local mass function around the half-mass radius closely resembles the actual global mass function. At the late stages of evolution, the mass of the evolved stars inside the core can be as high as 97 per cent of the total mass in this region. For the whole system, the evolved stars can compose up to 75 per cent of the total mass. The evolution of cluster anisotropy strongly depends on initial cluster concentration, IMF and the strength of the tidal field. The results presented are the first step in the direction of simulating the evolution of real globular clusters by means of the Monte Carlo method.     

The Hill stability of a binary or planetary system during encounters with a third inclined body

The dynamical interaction of a binary or planetary system and a third body moving on a parabolic orbit inclined to the system is discussed in terms of Hill stability for the full three-body problem. The situation arises in binary star disruption and exchange, in extrasolar planetary system disruption, exchange and capture. It is found that increasing the inclination of the third body decreases the Hill regions of stability. This makes exchange or disruption of the component masses more likely as does increasing the eccentricity of the binary.
The stability criteria are applied to determine possible disruption and capture distances for currently known extrasolar planetary systems.     

Geometric characterization of the Arjuna orbital domain

Arjuna‐type orbits are characterized by being Earth‐like, having both low‐eccentricity and low‐inclination. Objects following these trajectories experience repeated trappings in the 1:1 commensurability with the Earth and can become temporary Trojans, horseshoe librators, quasi‐satellites, and even transient natural satellites. Here, we review what we know about this peculiar dynamical group and use a Monte Carlo simulation to characterize geometrically the Arjuna orbital domain, studying its visibility both from the ground and with the European Space Agency Gaia spacecraft. The visibility analysis from the ground together with the discovery circumstances of known objects are used as proxies to estimate the current size of this population. The impact cross‐section of the Earth for minor bodies in this resonant group is also investigated. We find that, for ground‐based observations, the solar elongation at perigee of nearly half of these objects is less than 90°. They are best observed by space‐borne telescopes, but Gaia is not going to improve significantly the current discovery rate for members of this class. Our results suggest that the size of this population may have been underestimated by current models. On the other hand, their intrinsically low encounter velocities with the Earth induce a 10–1000‐fold increase in the impact cross‐section with respect to what is typical for objects in the Apollo or Aten asteroid populations. We estimate that their probability of capture as transient natural satellites of our planet is about 8 %. (© 2015 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

On the periodicity hypothesis of the ages of large impact craters

An analysis is made of the periodicity hypothesis of the ages of large craters, based on the compilation by Grieve with the addition of recently identified craters. A method earlier proposed by Broadbent is used to derive a period, and the significance of the derived period is tested by a Monte Carlo experiment. In accordance with the result of Stothers, the ages of large craters  ( D >30 km)  are shown to exhibit a period close to 37.5 Myr. Monte Carlo experiments show, however, that the derived period is far from being statistically significant. A subset of crater data earlier adopted by Napier for the purpose of similar investigation is also tested, and it is shown that they also exhibit a similar period at an almost identical level of confidence. A brief discussion is made of the relation between the derived period and that associated with faunal mass extinctions.     

Early scattering by Jupiter and its collision effects in the terrestrial zone

When Jupiter was on the order of three to ten Earth masses in size, there undoubtedly was a considerably larger mass of condensed matter in its zone, since Jupiter would have perturbed most of it to other parts of the solar system. Monte Carlo studies indicate a significant portion would have crossed the Earth's orbit. If the Earth and Moon had not yet fully formed, the probability of Earth-zone planetesimals being hit by this Jupiter-scattered material was high. Further Monte Carlo models of these collisions and their products indicate a significant portion of matter was heated to melting, even if less than 5% of the relative kinetic energy went into heat. The models include capture probabilities by an embryo Earth and a protolunar swarm. Because heat energy is correlated with comminution energy, and because the capture probability of the swarm is mass-dependent while the embryo's is not, the protolunar material suffered much higher heating on the average than did the proto-Earth material.     

Reflected iron line from a source above a Kerr black hole accretion disc

In this paper we present a fully relativistic approach to modelling both the continuum emission and the reflected fluorescent iron line from a primary X-ray source near a Kerr black hole. The X-ray source is located above an accretion disc orbiting around the black hole. The source is assumed to be a static point source located on an arbitrary position above the disc, on or off the axis of rotation. We carry out Monte Carlo simulations in order to estimate the iron line spectrum as well as its equivalent width. Because of the gravitational lensing effect, an enhancement of the iron line is expected when the primary source is located close to the central black hole. We find that for a source located on the axis of rotation the enhancement is relatively modest. An observer at inclination 30° would measure an equivalent width of ∼300 eV in the extreme case of a maximally rotating black hole and a source located at height 1.5 gravitational radii from the centre. This corresponds to an equivalent width enhancement factor of about 2 compared with the classical value where no lensing effect comes into play. However, when allowing the source to be located off the axis of rotation, much stronger enhancement can be obtained. In the extreme case of a maximally rotating black hole and a source located just above the approaching side of the disc, an observer at inclination 30° could measure an equivalent width as high as ∼1.5 keV (i.e., ∼10 times the classical value). We also find that observers located at high inclination angles observe a stronger line than observers at low inclination angles.     

How rapidly do neutron stars spin at birth? Constraints from archival X-ray observations of extragalactic supernovae

Traditionally, studies aimed at inferring the distribution of birth periods of neutron stars are based on radio surveys. Here we propose an independent method to constrain the pulsar spin periods at birth based on their X-ray luminosities. In particular, the observed luminosity distribution of supernovae (SNe) poses a constraint on the initial rotational energy of the embedded pulsars, via the     correlation found for radio pulsars, and under the assumption that this relation continues to hold beyond the observed range. We have extracted X-ray luminosities (or limits) for a large sample of historical SNe observed with Chandra , XMM and Swift , which have been firmly classified as core-collapse SNe. We have then compared these observational limits with the results of Monte Carlo simulations of the pulsar X-ray luminosity distribution for a range of values of the birth parameters. We find that a pulsar population dominated by millisecond periods at birth is ruled out by the data.     

The evolution of binary fractions in globular clusters

We study the evolution of binary stars in globular clusters using a new Monte Carlo approach combining a population synthesis code ( startrack ) and a simple treatment of dynamical interactions in the dense cluster core using a new tool for computing three- and four-body interactions ( fewbody ). We find that the combination of stellar evolution and dynamical interactions (binary–single and binary–binary) leads to a rapid depletion of the binary population in the cluster core. The maximum binary fraction today in the core of a typical dense cluster such as 47 Tuc, assuming an initial binary fraction of 100 per cent, is only ∼ 5–10 per cent. We show that this is in good agreement with recent Hubble Space Telescope observations of close binaries in the core of 47 Tuc, provided that a realistic distribution of binary periods is used to interpret the results. Our findings also have important consequences for the dynamical modelling of globular clusters, suggesting that 'realistic models' should incorporate much larger initial binary fractions than has usually been the case in the past.     

Towards multiple-star population synthesis

The multiplicities of stars, and some other properties, were collected recently by Eggleton & Tokovinin, for the set of 4559 stars with Hipparcos magnitude brighter than 6.0 (4558 excluding the Sun). In this paper I give a numerical recipe for constructing, by a Monte Carlo technique, a theoretical ensemble of multiple stars that resembles the observed sample. Only multiplicities up to eight are allowed; the observed set contains only multiplicities up to seven. In addition, recipes are suggested for dealing with the selection effects and observational uncertainties that attend the determination of multiplicity. These recipes imply, for example, that to achieve the observed average multiplicity of 1.53, it would be necessary to suppose that the real population has an average multiplicity slightly over 2.0.
This numerical model may be useful for (i) comparison with the results of star and star cluster formation theory, (ii) population synthesis that does not ignore multiplicity above 2 and (iii) initial conditions for dynamical cluster simulations.     

Solar flare neutron production and the angular dependence of the capture gamma-ray emission

The depth dependence of the production of neutrons and capture gamma-ray line emission are calculated by Monte Carlo simulation of the nuclear processes taking place when flare-accelerated ions interact with the solar atmosphere. The calculations also give the heliocentric-angular dependence of the 2.223 MeV neutron capture line emission as a function of accelerated-ion energy spectrum and angular distribution. These results are compared with observations to determine the energy spectrum shape and total ion number for various flares.     

A stochastic Monte Carlo approach to model real star cluster evolution – II. Self-consistent models and primordial binaries

The new approach outlined in Paper I to follow the individual formation and evolution of binaries in an evolving, equal point-mass star cluster is extended for the self-consistent treatment of relaxation and close three- and four-body encounters for many binaries (typically a few per cent of the initial number of stars in the cluster mass). The distribution of single stars is treated as a conducting gas sphere with a standard anisotropic gaseous model. A Monte Carlo technique is used to model the motion of binaries, their formation and subsequent hardening by close encounters, and their relaxation (dynamical friction) with single stars and other binaries. The results are a further approach towards a realistic model of globular clusters with primordial binaries without using special hardware. We present, as our main result, the self-consistent evolution of a cluster consisting of 300 000 equal point-mass stars, plus 30 000 equal-mass binaries over several hundred half-mass relaxation times, well into the phase where most of the binaries have been dissolved and evacuated from the core. The cluster evolution is about three times slower than found by Gao et al. Other features are rather comparable. At every moment we are able to show the individual distribution of binaries in the cluster.