Even though plenty of symbiotic stars (SySts) have been found in the Galactic field and nearby galaxies, not a single one has ever been confirmed in a Galactic globular cluster (GC). We investigate the lack of such systems in GCs for the first time by analysing 144 GC models evolved with the MOCCA code, which have different initial properties and are roughly representative of the Galactic GC population. We focus here on SySts formed through the wind-accretion channel, which can be consistently modelled in binary population synthesis codes. We found that the orbital periods of the majority of such SySts are sufficiently long (${\gtrsim10^3}$ d) so that, for very dense GC models, dynamical interactions play an important role in destroying their progenitors before the present day (${\sim11-12}$ Gyr). In less dense GC models, some SySts are still predicted to exist. However, these systems tend to be located far from the central parts (${\gtrsim70}$ per cent are far beyond the half-light radius) and are sufficiently rare (${\lesssim1}$ per GC per Myr), which makes their identification rather difficult in observational campaigns. We propose that future searches for SySts in GCs should be performed in the outskirts of nearby low-density GCs with sufficiently long half-mass relaxation times and relatively large Galactocentric distances. Finally, we obtained spectra of the candidate proposed in $\omega$ Cen (SOPS IV e-94) and showed that this object is most likely not a SySt.

Over the last few decades, exhaustive surveys of extra Galactic globular clusters (EGGCs) have become feasible. Only recently, kinematical information of globular clusters (GCs) were available through Gaia DR2 spectroscopy and also proper motions. On the other hand, simulations of GCs can provide detailed information about the dynamical evolution of the system. We present a preliminary study of EGGCs' properties for different dynamical evolutionary stages. We apply this study to 12 Gyr-old GCs simulated as part of the MOCCA Survey Database. Mimicking observational limits, we consider only a subssample of the models in the database, showing that it is possible to represent observed Galactic GCs. In order to distinguish between different dynamical states of EGGCs, at least three structural parameters are necessary. The best distinction is achieved by considering the central parameters, those being observational core radius, central surface brightness, ratio between central and half-mass velocity dispersion, or similarly considering the central color, the central V magnitude and the ratio between central and half-mass radius velocity dispersion, although such properties could be prohibitive with current technologies. A similar but less solid result is obtained considering the average properties at the half-light radius, perhaps accessible presently in the Local Group. Additionally, we mention that the color spread in EGGCs due to internal dynamical models, at fixed metallcity, could be just as important due to the spread in metallicity.

All ten LIGO/Virgo binary black hole (BH-BH) coalescences reported from the O1/O2 runs have near zero effective spins. There are only three potential explanations of this fact. If the BH spin magnitudes are large then (i) either both BH spin vectors must be nearly in the orbital plane or (ii) the spin angular momenta of the BHs must be oppositely directed and similar in magnitude. Or, (iii) the BH spin magnitudes are small. We test the third hypothesis within the framework of the classical isolated binary evolution scenario of the BH-BH merger formation. We test three models of angular momentum transport in massive stars: a mildly efficient transport by meridional currents (as employed in the Geneva code), an efficient transport by the Tayler-Spruit magnetic dynamo (as implemented in the MESA code), and a very-efficient transport (as proposed by Fuller et al.) to calculate natal BH spins. We allow for binary evolution to increase the BH spins through accretion and account for the potential spin-up of stars through tidal interactions. Additionally, we update the calculations of the stellar-origin BH masses, include revisions to the history of star formation and to the chemical evolution across cosmic time. We find that we can match simultaneously the observed BH-BH merger rate density, BH masses, and effective spins. Models with efficient angular momentum transport are favored. The updated stellar-mass weighted gas-phase metallicity evolution now used in our models appears to be a key in better reproducing the LIGO/Virgo merger rate estimate. Mass losses during the pair-instability pulsation supernova phase are likely overestimated if the merger GW170729 hosts a BH more massive than 50 Msun. We also estimate rate of BH-NS mergers from recent LIGO/Virgo observations. Our updated models of BH-BH, BH-NS and NS-NS mergers are now publicly available at www.syntheticuniverse.org.

The photometric properties that we could observe for Extra-Galactic Globular Clusters (EGGCs) are the integrated light of the system and for nearby EGGCs it also is possible to measure both half-light radii and the color spatial distribution, e.g. for areas smaller and larger than the half-light radius. No information about the internal dynamical state of the system could be directly obtained from observations. On the other hand, simulations of Globular Clusters (GCs) can provide detailed information about the dynamical evolution of the system. We present a preliminary study of EGGCs' photometric properties for different dynamical evolutionary stages. We apply this study to 12 Gyr old GCs simulated as part of the MOCCA Survey Database. We determine the magnitudes in different bands from their projected snapshots using the Flexible Stellar Population Synthesis (FSPS) code and we measure the half-light radii from the surface brightness.

We investigate the dissolution process for dynamically evolving star clusters embedded in an external tidal field by exploring the MOCCA Survey Database I, with focus on the presence and evolution of a stellar-mass black hole subsystem. We argue that the presence of a black hole subsystem can lead to the dissolution of tidally filling star clusters and this can be regarded as a third type of cluster dissolution mechanism (in addition to well-known mechanisms connected with strong mass loss due to stellar evolution and mass loss connected with the relaxation process). This third process is characterized by abrupt cluster dissolution connected with the loss of dynamical equilibrium. The abrupt dissolution is powered by strong energy generation from a stellar-mass black hole subsystem accompanied by tidal stripping. Additionally, we argue that such a mechanism should also work for even tidally under-filling clusters with top-heavy initial mass function. Observationally, star clusters which undergo dissolution powered by the third mechanism would look as a 'dark cluster' i.e. composed of stellar mass black holes surrounded by an expanding halo of luminous stars (Banerjee & Kroupa 2011), and they should be different from 'dark clusters' harbouring intermediate mass black holes as discussed by Askar et al. (2017a). An additional observational consequence of an operation of the third dissolution mechanism should be a larger than expected abundance of free floating black holes in the Galactic halo.

We investigate the dissolution process of star clusters embedded in an external tidal field and harboring a subsystem of stellar-mass black hole. For this purpose we analyzed the MOCCA models of real star clusters contained in the Mocca Survey Database I. We showed that the presence of a stellar-mass black hole subsystem in tidally filling star cluster can lead to abrupt cluster dissolution connected with the loss of cluster dynamical equilibrium. Such cluster dissolution can be regarded as a third type of cluster dissolution mechanism. We additionally argue that such a mechanism should also work for tidally under-filling clusters with a top-heavy initial mass function.

Sizeable number of stellar-mass black holes (BHs) in globular clusters (GCs) can strongly influence the dynamical evolution and observational properties of their host cluster. Using results from a large set of numerical simulations, we identify the key ingredients needed to sustain a sizeable population of BHs in GCs up to a Hubble time. We find that while BH natal kick prescriptions are essential in determining the initial retention fraction of BHs in GCs, the long-term survival of BHs is determined by the size, initial central density and half-mass relaxation time of the GC. Simulated GC models that contain many BHs are characterized by relatively low central surface brightness, large half-light and core radii values. We also discuss novel ways to compare simulated results with available observational data to identify GCs that are most likely to contain many BHs.

We have been investigating populations of cataclysmic variables (CVs) in a set of more than 300 globular cluster (GC) models evolved with the MOCCA code. One of the main questions we have intended to answer is whether most CVs in GCs are dynamically formed or not. Contrary to what has been argued for a long time, we found that dynamical destruction of primordial CV progenitors is much stronger in GCs than dynamical formation of CVs. In particular, we found that, on average, the detectable CV population is predominantly composed of CVs formed via a typical common envelope phase (> 70 per cent). However, core-collapsed models tend to have higher fractions of bright CVs than non-core-collapsed ones, which suggests then that the formation of CVs is indeed slightly favoured through strong dynamical interactions in core-collapsed GCs, due to the high stellar densities in their cores.

In this paper we explore the interplay between intermediate-mass black holes (IMBHs) and their nursing globular clusters (GCs), taking advantage of over 2000 Monte Carlo GC models. We find that the average density of IMBHs sphere of influence can be uniquely connected to the host GCs luminosity and half-light radius via a fundamental plane. We propose a statistical approach to systematically identify potential Galactic GCs harbouring either an IMBH or a massive subsystem comprised of stellar BHs. Our models show that the IMBH is often bound to a stellar companion or a stellar BH, which can lead to tidal disruption events or to low-frequency gravitational waves. We show that GCs orbiting close to the Galactic Centre have a larger probability to witness IMBH formation during their early evolution. These low-orbit GCs can deliver several IMBHs into the galaxy innermost regions, with potential impact on both electromagnetic and GW emission. We discuss potential connections between IMBHs and SMBHs inhabiting galactic nuclei, exploring the possibility that in some cases they share similar formation pathways.

We briefly describe and discuss the set-up of the project MOCCA Survey Database I. The database contains more than 2000 Monte Carlo models of evolution of real star cluster performed with the MOCCA code. Then, we very briefly discuss results of analysis of the database regarding the following projects: formation of intermediate mass black holes, abrupt cluster dissolution harboring black hole subsystems, retention fraction of black hole - black hole mergers, and tidal disruption events with intermediate mass black holes.

Machine learning is a powerful technique, becoming increasingly popular in astrophysics. In this paper, we apply machine learning to more than a thousand globular cluster (GC) models simulated as part of the 'MOCCA-Survey Database I' project in order to correlate present-day observable properties with the presence of a subsystem of stellar mass black holes (BHs). The machine learning model is then applied to available observed parameters for Galactic GCs to identify which of them that are most likely to be hosting a sizeable number of BHs and reveal insights into what properties lead to the formation of BH subsystems. With our machine learning model, we were able to shortlist 21 Galactic GCs that are most likely to contain a BH subsystem. We show that the clusters shortlisted by the machine learning classifier include those in which BH candidates have been observed (M22, M10 and NGC 3201) and that our results line up well with independent simulations and previous studies that manually compared simulated GC models with observed properties of Galactic GCs. These results can be useful for observers searching for elusive stellar mass BH candidates in GCs and further our understanding of the role BHs play in GC evolution. In addition, we have released an online tool that allows one to get predictions from our model after they input observable properties.

About a century ago, Jeans (1919) discovered that if binary stars reach a state approximating energy equipartition, for example through many dynamical encounters that exchange energy, their eccentricity distribution can be described by : dN/de = 2e. This is referred to as the thermal eccentricity distribution, and has been widely used for initial conditions in theoretical investigations of binary stars. However, observations suggest that the eccentricity distributions of most observed binaries, and particularly those with masses < 5 Msun, are flatter than thermal and follow more closely to a uniform distribution. Nonetheless, it is often argued that dynamical interactions in a star cluster would quickly thermalize the binaries, which could justify imposing a thermal eccentricity distribution at birth for all binaries. In this paper we investigate the validity of this assumption. We develop our own rapid semi-analytic model for binary evolution in star clusters, and also compare with detailed N-body and Monte Carlo star cluster models. We show that, for nearly all binaries, dynamical encounters fail to convert an initially uniform eccentricity distribution to thermal within a star cluster's lifetime. Thus, if a thermal eccentricity distribution is observed, it is likely imprinted upon formation rather than through subsequent long-term dynamical processing. Theoretical investigations that initialize all binaries with a thermal distribution will make incorrect predictions for the evolution of the binary population. Such models may overpredict the merger rate for binaries with modest orbital separations by a factor of about two.

We investigate here populations of cataclysmic variables (CVs) in a set of 288 globular cluster (GC) models evolved with the MOCCA code. This is by far the largest sample of GC models ever analysed with respect to CVs. Contrary to what has been argued for a long time, we found that dynamical destruction of primordial CV progenitors is much stronger in GCs than dynamical formation of CVs, and that dynamically formed CVs and CVs formed under no/weak influence of dynamics have similar white dwarf mass distributions. In addition, we found that, on average, the detectable CV population is predominantly composed of CVs formed via typical common envelope phase (CEP) ($\gtrsim70$ per cent), that only $\approx2-4$ per cent of all CVs in a GC is likely to be detectable, and that core-collapsed models tend to have higher fractions of bright CVs than non-core-collapsed ones. We also consistently show, for the first time, that the properties of bright and faint CVs can be understood by means of the pre-CV and CV formation rates, their properties at their formation times and cluster half-mass relaxation times. Finally, we show that models following the initial binary population proposed by Kroupa and set with low CEP efficiency better reproduce the observed amount of CVs and CV candidates in NGC 6397, NGC 6752 and 47 Tuc. To progress with comparisons, the essential next step is to properly characterize the candidates as CVs (e.g. by obtaining orbital periods and mass ratios).

Anisotropy of gravitational wave (GW) emission results in a net momentum gained by the black hole (BH) merger product, leading to a recoil velocity up to $\sim10^3\text{ km s}^{-1}$, which may kick it out of a globular cluster (GC). We estimate GW kick retention fractions of merger products assuming different models for BH spin magnitude and orientation (MS0 - random, MS1 - spin as a function of mass and metalicity, MS2 - constant value of $0.5$). We check how they depend on BH-BH merger time and properties of the cluster. We analyze the implications of GW kick retention fractions on intermediate massive BH (IMBH) formation by repeated mergers in a GC. We also calculate final spin of the merger product, and investigate how it correlates with effective spin of the binary. We used data from MOCCA (MOnte Carlo Cluster simulAtor) GC simulations to get a realistic sample of BH-BH mergers, assigned each BH spin value according to a studied model, and calculated recoil velocity and final spin based on most recent theoretical formulas. We discovered that for physically motivated models, GW kick retention fractions are about $30\%$ and display small dependence on assumptions about spin, but are much more prone to cluster properties. In particular, we discovered a strong dependence of GW kick retention fractions on cluster density. We also show that GW kick retention fractions are high in final life stages of the cluster, but low at the beginning. Finally, we derive formulas connecting final spin with effective spin for primordial binaries, and with maximal effective spin for dynamical binaries.

The dense environment of globular clusters (GCs) can facilitate the formation of binary black holes (BBHs), some of which can merge with gravitational waves (GW) within the age of the Universe. We have performed a survey of Monte-Carlo simulations following the dynamical evolution of GCs with different masses, sizes and binary fractions and explored the impact of the host GC properties on the formation of BBH mergers. We find that the number of BBH mergers from GCs is determined by the GC's initial mass, size and primordial binary fraction. We identify two groups of BBH mergers: a primordial group whose formation does not depend on cluster's dynamics and a dynamical group whose formation is driven by the cluster's dynamical evolution. We show how the BBH origin affects the BBH mergers' main properties such as the chirp mass and merging time distributions. We provide analytic expressions for the dependence of the number of BBH mergers from individual GCs on the main cluster's structural properties and the time evolution of the merger rates of these BBHs. These expressions provide an essential ingredient for a general framework allowing to estimate the merger rate density. Using the relations found in our study, we find a local merger rate density of 0.18$-$1.8 ${\rm Gpc}^{-3}{\rm yr}^{-1}$ for primordial BBH mergers and 0.6$-$18 ${\rm Gpc}^{-3}{\rm yr}^{-1}$ for dynamical BBH mergers, depending on the GC mass and size distributions, initial binary fraction and the number density of GCs in the Universe.

The predicted and observed space density of cataclysmic variables (CVs) have been for a long time discrepant by at least an order of magnitude. The standard model of CV evolution predicts that the vast majority of CVs should be period bouncers, whose space density has been recently measured to be $\rho \lesssim 2 \times 10^{-5}$ pc$^{-3}$. We performed population synthesis of CVs using an updated version of the Binary Stellar Evolution (BSE) code for single and binary star evolution. We find that the recently suggested empirical prescription of consequential angular momentum loss (CAML) brings into agreement predicted and observed space densities of CVs and period bouncers. To progress with our understanding of CV evolution it is crucial to understand the physical mechanism behind empirical CAML. Our changes to the BSE code are also provided in details, which will allow the community to accurately model mass transfer in interacting binaries in which degenerate objects accrete from low-mass main-sequence donor stars.

There have been increasing theoretical speculations and observational indications that certain globular clusters (GCs) could contain a sizeable population of stellar mass black holes (BHs). In this paper, we shortlist at least 29 Galactic GCs that could be hosting a subsystem of BHs (BHS). In a companion paper, we analysed results from a wide array of GC models (simulated with the MOCCA code for cluster simulations) that retained few tens to several hundreds of BHs at 12 Gyr and showed that the properties of the BHS in those GCs correlate with the GC's observable properties. Building on those results, we use available observational properties of 140 Galactic GCs to identify 29 GCs that could potentially be harbouring up to a few hundreds of BHs. Utilizing observational properties and theoretical scaling relations, we estimate the density, size and mass of the BHS in these GCs. We also calculate the total number of BHs and the fraction of BHs contained in a binary system for our shortlisted Galactic GCs. Additionally, we mention other Galactic GCs that could also contain significant number of single BHs or BHs in binary systems.

In this paper, we discuss how globular clusters (GCs) structural and observational properties can be used to infer the presence of a black hole system (BHS) inhabiting their inner regions. We propose a novel way to identify the BHS size, defined as the GC radius containing a mass contributed equally from stars and stellar BHs. Using this definition, similar to the well-known concept of "influence radius", we found a "fundamental plane" connecting the BHS typical density with the GC central surface density profile, total luminosity and observational half-mass radius. Our approach allows us to define a unique way to connect the observational GCs parameters with their dark content. Comparing our results with observed Milky Way GCs, we found that many of them likely host, at the present time, as many as several hundreds of BHs. These BHS are characterized by a relatively low typical density, $\rho_\bhs \sim 10-10^5\Ms$ pc$^{-3}$ and composed of relatively massive BHs, with average masses in the range $m_\bhs = 14-22\Ms$. We also show that a similar approach can be used to find Milky Way GCs potentially hosting an intermediate-mass black hole.

The first neutron star-neutron star (NS-NS) merger was discovered on August 17, 2017 through gravitational waves (GW170817) and followed with electromagnetic observations. This merger was detected in an old elliptical galaxy with no recent star formation. We perform a suite of numerical calculations to understand the formation mechanism of this merger. We probe three leading formation mechanisms of double compact objects: classical isolated binary star evolution, dynamical evolution in globular clusters and nuclear cluster formation to test whether they are likely to produce NS-NS mergers in old host galaxies. Our simulations with optimistic assumptions show current NS-NS merger rates at the level of 10^-2 yr^-1 from binary stars, 5 x 10^-5 yr^-1 from globular clusters and 10^-5 yr^-1 from nuclear clusters for all local elliptical galaxies (within 100 Mpc^3). These models are thus in tension with the detection of GW170817 with an observed rate 1.5 yr^-1 (per 100 Mpc^3; LIGO/Virgo estimate). Our results imply that either (i) the detection of GW170817 by LIGO/Virgo at their current sensitivity in an elliptical galaxy is a statistical coincidence; or that (ii) physics in at least one of our three models is incomplete in the context of the evolution of stars that can form NS-NS mergers; or that (iii) another very efficient (unknown) formation channel with a long delay time between star formation and merger is at play.

We study the gravitational wave (GW) frequency and chirp mass distribution of binary black hole (BBH) mergers assembled through three-body interactions in globular clusters (GCs), when GW emission at the 2.5 post-Newtonian (PN) level is included in the $N$-body equation-of-motion (EOM). From performing $\sim 2.5\times10^{6}$ PN binary-single interactions based on GC data from the `MOCCA-Survey Database I' project, and by the use of analytical methods, we find that $5-10\%$ of all the three-body assembled GC BBH mergers have a GW frequency at formation that is $\gtrsim 10^{-1}$ Hz, implying they enter the LIGO band without having drifted through the LISA band first. If PN terms are not included in the EOM, one finds instead that all BBH mergers drifts through both LISA and LIGO. As the fraction of BBH mergers that only show up in LIGO is expected to be $\sim 0\%$ for standard field binary BBH mergers, future joint measurements with LISA and LIGO can be used to gain insight into the formation of BBH mergers.

We estimate the population of eccentric gravitational wave (GW) binary black hole (BBH) mergers forming during binary-single interactions in globular clusters (GCs), using ~ 800 GC models that were evolved using the MOCCA code for star cluster simulations as part of the MOCCA-Survey Database I project. By re-simulating binary-single interactions (only involving 3 BHs) extracted from this set of GC models using an N-body code that includes GW emission at the 2.5 post-Newtonian level, we find that ~ 10% of all the BBHs assembled in our GC models that merge at present time form during chaotic binary-single interactions, and that about half of this sample have an eccentricity > 0.1 at 10 Hz. We explicitly show that this derived rate of eccentric mergers is ~ 100 times higher than one would find with a purely Newtonian N-body code. Furthermore, we demonstrate that the eccentric fraction can be accurately estimated using a simple analytical formalism when the interacting BHs are of similar mass; a result that serves as the first successful analytical description of eccentric GW mergers forming during three-body interactions in realistic GCs.

We introduce and present results from the COCOA (Cluster simulatiOn Comparison with ObservAtions) code that has been developed to create idealized mock photometric observations using results from numerical simulations of star cluster evolution. COCOA is able to present the output of realistic numerical simulations of star clusters carried out using Monte Carlo or \textit{N}-body codes in a way that is useful for direct comparison with photometric observations. In this paper, we describe the COCOA code and demonstrate its different applications by utilizing globular cluster (GC) models simulated with the MOCCA (MOnte Carlo Cluster simulAtor) code. COCOA is used to synthetically observe these different GC models with optical telescopes, perform PSF photometry and subsequently produce observed colour magnitude diagrams. We also use COCOA to compare the results from synthetic observations of a cluster model that has the same age and metallicity as the Galactic GC NGC 2808 with observations of the same cluster carried out with a 2.2 meter optical telescope. We find that COCOA can effectively simulate realistic observations and recover photometric data. COCOA has numerous scientific applications that maybe be helpful for both theoreticians and observers that work on star clusters. Plans for further improving and developing the code are also discussed in this paper.

In order to allow a better understanding of the origin of Galactic field populations, dynamical equivalence of stellar-dynamical systems has been postulated by Kroupa and Belloni et al. to allow mapping of solutions of the initial conditions of embedded clusters such that they yield, after a period of dynamical processing, the Galactic field population. Dynamically equivalent systems are defined to initially and finally have the same distribution functions of periods, mass ratios and eccentricities of binary stars. Here we search for dynamically equivalent clusters using the {\sc mocca} code. The simulations confirm that dynamically equivalent solutions indeed exist. The result is that the solution space is next to identical to the radius--mass relation of Marks \& Kroupa, $\left( r_h/{\rm pc} \right)= 0.1^{+0.07}_{-0.04}\, \left( M_{\rm ecl}/{\rm M}_\odot \right)^{0.13\pm0.04}$. This relation is in good agreement with the observed density of molecular cloud clumps. According to the solutions, the time-scale to reach dynamical equivalence is about 0.5~Myr which is, interestingly, consistent with the lifetime of ultra-compact HII regions and the time-scale needed for gas expulsion to be active in observed very young clusters as based on their dynamical modelling.

Numerical simulations were carried out to study the origin of multiple stellar populations in the intermediate-age clusters NGC 411 and NGC 1806 in the Magellanic Clouds. We performed NBODY6++ simulations based on two different formation scenarios, an ad hoc formation model where second-generation (SG) stars are formed inside a cluster of first-generation (FG) stars using the gas accumulated from the external intergalactic medium and a minor merger model of unequal mass ($M_{SG}$/$M_{FG}$ ~5-10%) clusters with an age difference of a few hundred million years. We compared our results such as the radial profile of the SG-to-FG number ratio with observations on the assumption that the SG stars in the observations are composed of cluster members, and confirmed that both the ad hoc formation and merger scenarios reproduce the observed radial trend of the SG-to-FG number ratio which shows less centrally concentrated SG than FG stars. It is difficult to constrain the formation scenario for the multiple populations by only using the spatial distribution of the SG stars. SG stars originating from the merger scenario show a significant velocity anisotropy and rotational features compared to those from the ad hoc formation scenario. Thus, observations aimed at kinematic properties like velocity anisotropy or rotational velocities for SG stars should be obtained to better understand the formation of the multiple populations in these clusters. This is, however, beyond current instrumentation capabilities.

Colour-magnitude diagrams (CMDs) are powerful tools that might be used to infer stellar properties in globular clusters (GCs), for example, the binary fraction and their mass ratio ($q$) distribution. In the past few years, observations have revealed that q distributions of GC main-sequence (MS) binaries are generally flat, and a distribution characterized by a strong increase towards q $\approx$ 1 is not typical in GCs. In numerical simulations of GC evolution with the initial binary population (IBP) described by Kroupa, synthetic CMD colour distributions exhibit a peak associated with binaries that have q $\approx$ 1. While the Kroupa IBP reproduces binary properties in star-forming regions, clusters and the Galactic field, the peak in the q distribution towards q $\approx$ 1 observed for GC simulations is not consistent with distributions derived from observations. The objective of this paper is to refine and further improve the physical formulation of pre-main-sequence eigenevolution proposed by Kroupa in order to achieve CMD colour distributions of simulated GC models similar to those observed in real GCs, and to get a similarly good agreement with binary properties for late-type binaries in the Galactic field. We present in this paper a modified Kroupa IBP, in which early-type stars follow observational distributions, and late-type stars are generated according to slightly modified pre-main-sequence eigenevolution prescriptions. Our modifications not only lead to a qualitatively good agreement with respect to long-term observations of late-type binaries in the Galactic field, but also resolve the above-mentioned problem related to binary distributions in GC models.

The mass evaporation rate of globular clusters evolving in a strong Galactic tidal field is derived through the analysis of large, multi-mass $N$-body simulations. For comparison, we also study the same evaporation rates using MOCCA Monte Carlo models for globular cluster evolution. Our results show that the mass evaporation rate is a dynamical value; that is, far from a constant single number found in earlier analytical work and commonly used in the literature. Moreover, the evaporation rate derived with these simulations is higher than values previously published. These models also show that the value of the mass evaporation rate depends on the strength of the tidal field. We give an analytical estimate of the mass evaporation rate as a function of time and galactocentric distance xi(R_{GC},t). Upon extrapolating this formula to smaller R_{GC} values, our results provide tentative evidence for a very high xi value at small R_{GC}. Our results suggest that the corresponding mass loss in the inner Galactic potential could be high and it should be accounted for when star clusters pass within it. This has direct relevance to nuclear cluster formation/growth via the infall of globular clusters through dynamical friction. As an illustrative example, we estimate how the evaporation rate increases for a ~10^5 solar masses globular cluster that decays through dynamical friction into the galactic centre. We discuss the findings of this work in relation to the formation of nuclear star clusters by inspiraling globular clusters.

In this third of a series of papers related to cataclysmic variables (CVs) and related objects, we analyse the population of CVs in a set of 12 globular cluster models evolved with the MOCCA Monte Carlo code, for two initial binary populations (IBPs), two choices of common-envelope phase (CEP) parameters, and three different models for the evolution of CVs and the treatment of angular momentum loss. When more realistic models and parameters are considered, we find that present-day cluster CV duty cycles are extremely-low ($\lesssim 0.1$ per cent) which makes their detection during outbursts rather difficult. Additionally, the IBP plays a significant role in shaping the CV population properties, and models that follow the Kroupa IBP are less affected by enhanced angular momentum loss. We also predict from our simulations that CVs formed dynamically in the past few Gyr (massive CVs) correspond to bright CVs (as expected), and that faint CVs formed several Gyr ago (dynamically or not) represent the overwhelming majority. Regarding the CV formation rate, we rule out the notion that it is similar irrespective of the cluster properties. Finally, we discuss the differences in the present-day CV properties related to the IBPs, the initial cluster conditions, the CEP parameters, formation channels, the CV evolution models, and the angular momentum loss treatments.

The detection of intermediate mass black holes (IMBHs) in Galactic globular clusters (GCs) has so far been controversial. In order to characterize the effectiveness of integrated-light spectroscopy through integral field units, we analyze realistic mock data generated from state-of-the-art Monte Carlo simulations of GCs with a central IMBH, considering different setups and conditions varying IMBH mass, cluster distance, and accuracy in determination of the center. The mock observations are modeled with isotropic Jeans models to assess the success rate in identifying the IMBH presence, which we find to be primarily dependent on IMBH mass. However, even for a IMBH of considerable mass (3% of the total GC mass), the analysis does not yield conclusive results in 1 out of 5 cases, because of shot noise due to bright stars close to the IMBH line-of-sight. This stochastic variability in the modeling outcome grows with decreasing BH mass, with approximately 3 failures out of 4 for IMBHs with 0.1% of total GC mass. Finally, we find that our analysis is generally unable to exclude at 68% confidence an IMBH with mass of $10^3~M_\odot$ in snapshots without a central BH. Interestingly, our results are not sensitive to GC distance within 5-20 kpc, nor to mis-identification of the GC center by less than 2'' (<20% of the core radius). These findings highlight the value of ground-based integral field spectroscopy for large GC surveys, where systematic failures can be accounted for, but stress the importance of discrete kinematic measurements that are less affected by stochasticity induced by bright stars.

This is the second in a series of papers associated with cataclysmic variables (CVs) and related objects, formed in a suite of simulations for globular cluster evolution performed with the MOCCA Monte Carlo code. We study the properties of our simulated CV populations throughout the entire cluster evolution. We find that dynamics extends the range of binary CV progenitor properties, causing CV formation from binary progenitors that would otherwise not become CVs. The CV formation rate in our simulations can be separated into two regimes: an initial burst ($\lesssim$ 1 Gyr) connected with the formation of the most massive white dwarfs, followed by a nearly constant formation rate. This result holds for all models regardless of the adopted initial conditions, even when most CVs form dynamically. Given the cluster age-dependence of CV properties, we argue that direct comparisons to observed Galactic field CVs could be misleading, since cluster CVs can be up to 4 times older than their field counterparts. Our results also illustrate that, due mainly to unstable mass transfer, some CVs that form in our simulations are destroyed before the present-day. Finally, some field CVs might have originated from globular clusters, as found in our simulations, although the fraction of such escapers should be small relative to the entire Galactic field CV population.

The formation and evolution of X-ray sources in globular clusters is likely to be affected by the cluster internal dynamics and the stellar interactions in the cluster dense environment.Several observational studies have revealed a correlation between the number of X-ray sources and the stellar encounter rate and provided evidence of the role of dynamics in the formation of X-ray binaries. We have performed a survey of Monte-Carlo simulations aimed at exploring the connection between the dynamics and formation of cataclysmic variables (CVs) and the origin of the observed correlation between the number of these objects, $N_{\rm cv}$, and the stellar encounter rate, $\Gamma$.The results of our simulations show a correlation between $N_{\rm cv}$ and $\Gamma$ as found in observational data, illustrate the essential role played by dynamics, and shed light on the dynamical history behind this correlation. CVs in our simulations are more centrally concentrated than single stars with masses close to those of turn-off stars, although this trend is stronger for CVs formed from primordial binaries undergoing exchange encounters, which include a population of more massive CVs absent in the group of CVs formed from binaries not suffering any component exchange.

We describe the dynamical evolution of a unique type of dark star cluster model in which the majority of the cluster mass at Hubble time is dominated by an intermediate-mass black hole (IMBH). We analyzed results from about 2000 star cluster models (Survey Database I) simulated using the Monte-Carlo code MOCCA and identified these dark star cluster models. Taking one of these models, we apply the method of simulating realistic "mock observations" by utilizing the COCOA and SISCO codes to obtain the photometric and kinematic observational properties of the dark star cluster model at 12 Gyr. We find that the perplexing Galactic globular cluster NGC 6535 closely matches the observational photometric and kinematic properties of the dark star cluster model presented in this paper. Based on our analysis and currently observed properties of NGC 6535, we suggest that this globular cluster could potentially harbour an IMBH. If it exists, the presence of this IMBH can be detected robustly with proposed kinematic observations of NGC 6535.

In this first of a series of papers, we utilize results for around two thousand star cluster models simulated using the MOCCA code for star cluster evolution (Survey Database I) to determine the astrophysical properties and local merger rate densities for coalescing binary black holes (BBHs) originating from globular clusters (GCs). We extracted information for all coalescing BBHs that escape the cluster models and subsequently merge within a Hubble time along with BBHs that are retained in our GC models and merge inside the cluster via gravitational wave (GW) emission. By obtaining results from a substantial number of realistic star cluster models that cover different initial parameters, we have an extremely large statistical sample of BBHs with stellar mass and massive stellar BH ($\lesssim 100M_{\odot}$) components that merge within a Hubble time. Using this data, we estimate local merger rate densities for these BBHs originating from GCs to be at least 5.4 ${\rm Gpc}^{-3}\,{\rm yr}^{-1}$

In this paper, which is the first in a series of papers associated with cataclysmic variables and related objects, we introduce the CATUABA code, a numerical machinery written for analysis of the MOCCA simulations, and show some first results by investigating the present-day population of cataclysmic variables in globular clusters. Emphasis was given on their properties and the observational selection effects when observing and detecting them. In this work we analysed in this work six models, including three with Kroupa distributions of the initial binaries. We found that for models with Kroupa initial distributions, considering the standard value of the efficiency of the common envelope phase adopted in BSE, no single cataclysmic variable was formed only via binary stellar evolution, i. e., in order to form them, strong dynamical interactions have to take place. We show and explain why this is inconsistent with observational and theoretical results. Our results indicate that the population of cataclysmic variables in globular clusters is, mainly, in the last stage of their evolution and observational selection effects can drastically change the expected number of observed cataclysmic variables. We show that the probability of observing them during the outbursts is extremely small and conclude that the best way of looking for cataclysmic variables in globular clusters is by searching for variabilities during quiescence, instead of during outbursts. For that, one would need a very deep observation which could reach magnitudes $\gtrsim$ 27 mag. Finally, we argue that cataclysmic variables in globular clusters are not necessarily magnetic.

We introduce and present preliminary results from COCOA (Cluster simulatiOn Comparison with ObservAtions) code for a star cluster after 12 Gyrs of evolution simulated using the MOCCA code. The COCOA code is being developed to quickly compare results of numerical simulations of star clusters with observational data. We use COCOA to obtain parameters of the projected cluster model. For comparison, a FITS file of the projected cluster was provided to observers so that they could use their observational methods and techniques to obtain cluster parameters. The results show that the similarity of cluster parameters obtained through numerical simulations and observations depends significantly on the quality of observational data and photometric accuracy.

We discuss a new scenario for the formation of intermediate mass black holes in dense star clusters. In this scenario, intermediate mass black holes are formed as a result of dynamical interactions of hard binaries containing a stellar mass black hole, with other stars and binaries. We discuss the necessary conditions to initiate the process of intermediate mass black hole formation and the influence of an intermediate mass black hole on the host global globular cluster properties. We discuss two scenarios for intermediate mass black hole formation. The SLOW and FAST scenarios. They occur later or earlier in the cluster evolution and require smaller or extremely large central densities, respectively. In our simulations, the formation of intermediate mass black holes is highly stochastic. In general, higher formation probabilities follow from larger cluster concentrations (i.e. central densities). We further discuss possible observational signatures of the presence of intermediate mass black holes in globular clusters that follow from our simulations. These include the spatial and kinematic structure of the host cluster, possible radio, X-ray and gravitational wave emissions due to dynamical collisions or mass-transfer and the creation of hypervelocity main sequence escapers during strong dynamical interactions between binaries and an intermediate mass black hole. All simulations discussed in this paper were performed with the MOCCA Monte Carlo code. MOCCA accurately follows most of the important physical processes that occur during the dynamical evolution of star clusters but, as with other dynamical codes, it approximates the dissipative processes connected with stellar collisions and binary mergers.

We very briefly discuss proposed in the literature possible scenarios for intermediate mass black holes (IMBH) formation in globular clusters. We also discuss the results of the MOCCA simulations of about 2000 models (BigSurvey) regarding the distribution of events connected with electromagnetic and gravitational radiations, namely: mass transfer on IMBH, collisions and mergers with IMBH and mergers with IMBH due to gravitational radiation. The rates of these events are very small, so their observation is very improbable.

In this first investigation of the MOCCA database with respect to cataclysmic variables, we found that for models with Kroupa initial distributions, considering the standard value of the efficiency of the common-envelope phase adopted in BSE, no single cataclysmic variable was formed only via binary stellar evolution, i. e., in order to form them, strong dynamical interactions have to take place. Our results also indicate that the population of cataclysmic variables in globular clusters are, mainly, in the last stage of their evolution and observational selection effects can change drastically the expected number and properties of observed cataclysmic variables.

The paper presents an analysis of properties of populations of blue stragglers (BSs) in evolving globular clusters, based on numerical simulations done with the MOCCA code for various initial globular clusters conditions. We find that various populations of BSs strongly depend on the initial semi-major axes distributions. With a significant number of compact binaries, the number of evolutionary BSs can be also significant. In turn, for semi-major axes distributions preferring binaries with wider orbits, dynamical BSs are the dominant ones. Their formation scenario is very distinct: for wide binaries the number of dynamical interactions is significantly larger. Most interactions are weak and increase only slightly the eccentricities. However, due to a large number of such interactions, the eccentricities of a number of binaries finally get so large that the stars collide. We study how larger initial clusters' concentrations influence the BSs. Besides the expected increase of the number of dynamically created BSs (for denser GCs the probabilities of strong dynamical interactions and collisions are higher), we find that the number of the evolutionary BSs is not affected even by very high initial concentrations. This has a very important implication on observations - it supports the theory that the evolutionary BSs are the result of the unperturbed evolution of the primordial binaries. In addition, the paper presents the evolution of the ratio between the number of BSs in binaries and as single stars R_B/S. For a vast diversity of models, the ratio R_B/S approaches the value ~0.4. Additionally, we identified two subgroups which differ in the initial semi-major axes distributions. The first group starts with a high ratio R_B/S, it decreases with time and settles around 0.4. The second group starts with lower values of the ratio R_B/S ... . (etc., abstract continues)

The paper presents an analysis of formation mechanism and properties of spatial distributions of blue stragglers in evolving globular clusters, based on numerical simulations done with the MOCCA code. First, there are presented N-body and MOCCA simulations which try to reproduce the simulations presented by Ferraro (2012). Then, the agreement between N-body and the MOCCA code is shown. Finally, we discuss the formation process of the bimodal distribution. We report that so-called bimodal spatial distribution of blue stragglers is a very transient feature. It is formed for one snapshot in time and it can easily vanish in the next one. Moreover, we show that the radius of avoidance proposed by Ferraro (2012) goes out of sync with the apparent minimum of the bimodal distribution after about two half-mass relaxation times. This finding creates a real challenge for the dynamical clock, which uses this radius to determine the dynamical age of globular clusters. Additionally, the paper discusses a few important problems concerning the apparent visibilities of the bimodal distributions which have to be taken into account while studying the spatial distributions of blue stragglers.

Introducing the DRAGON simulation project, we present direct $N$-body simulations of four massive globular clusters (GCs) with $10^6$ stars and 5$\%$ primordial binaries at a high level of accuracy and realism. The GC evolution is computed with NBODY6++GPU and follows the dynamical and stellar evolution of individual stars and binaries, kicks of neutron stars and black holes, and the effect of a tidal field. We investigate the evolution of the luminous (stellar) and dark (faint stars and stellar remnants) GC components and create mock observations of the simulations (i.e. photometry, color-magnitude diagrams, surface brightness and velocity dispersion profiles). By connecting internal processes to observable features we highlight the formation of a long-lived 'dark' nuclear subsystem made of black holes (BHs), which results in a two-component structure. The inner core is dominated by the BH subsystem and experiences a core collapse phase within the first Gyr. It can be detected in the stellar (luminous) line-of-sight velocity dispersion profiles. The outer extended core - commonly observed in the (luminous) surface brightness profiles - shows no collapse features and is continuously expanding. We demonstrate how a King (1966) model fit to observed clusters might help identify the presence of post core-collapse BH subsystems. For global observables like core and half-mass radii the direct simulations agree well with Monte-Carlo models. Variations in the initial mass function can result in significantly different GC properties (e.g. density distributions) driven by varying amounts of early mass loss and the number of forming BHs.

We will discuss the evolution of star clusters with an large initial binary fraction, up to 95%. The initial binary population is chosen to follow the invariant orbital-parameter distributions suggested by Kroupa (1995). The Monte Carlo MOCCA simulations of star cluster evolution are compared to the observations of Milone et al. (2012) for photometric binaries. It is demonstrated that the observed dependence on cluster mass of both the binary fraction and the ratio of the binary fractions inside and outside of the half mass radius are well recovered by the MOCCA simulations. This is due to a rapid decrease in the initial binary fraction due to the strong density-dependent destruction of wide binaries described by Marks, Kroupa & Oh (2011). We also discuss a new scenario for the formation of intermediate mass black holes in dense star clusters. In this scenario, intermediate mass black holes are formed as a result of dynamical interactions of hard binaries containing a stellar mass black hole, with other stars and binaries. We will discuss the necessary conditions to initiate the process of intermediate mass black hole formation and the dependence of its mass accretion rate on the global cluster properties.

In this paper, we discuss the origin of the observed correlation between cluster concentration c and present-day mass function (PDMF) slope {\alpha} reported by De Marchi, Paresce & Pulone. This relation can either be reproduced from universal initial conditions combined with some dynamical mechanism(s) that alter(s) the cluster structure and mass function over time, or it must arise early on in the cluster lifetime, such as during the gas-embedded phase of cluster formation. Using a combination of Monte Carlo and N-body models for globular cluster evolution performed with the MOCCA and NBODY6 codes, respectively, we explore a number of dynamical mechanisms that could affect the observed relation. For the range of initial conditions considered here, our results are consistent with an universal initial binary fraction \approx 10% (which does not, however, preclude 100%) and an universal initial stellar mass function resembling the standard Kroupa distribution. Most of the dispersion observed in the c-{\alpha} relation can be attributed to two-body relaxation and Galactic tides. However, dynamical processes alone could not have reproduced the dispersion in concentration, and we require at least some correlation between the initial concentration and the total cluster mass. We argue that the origin of this trend could be connected to the gas-embedded phase of cluster evolution.

(abridged) In this paper, we constrain the properties of primordial binary populations in Galactic globular clusters using the MOCCA Monte Carlo code for cluster evolution. Our results are compared to the observations of Milone et al. (2012) using the photometric binary populations as proxies for the true underlying distributions, in order to test the hypothesis that the data are consistent with an universal initial binary fraction near unity and the binary orbital parameter distributions of Kroupa (1995). With the exception of a few possible outliers, we find that the data are to first-order consistent with the universality hypothesis. Specifically, the present-day binary fractions inside the half-mass radius r$_{\rm h}$ can be reproduced assuming either high initial binary fractions near unity with a dominant soft binary component as in the Kroupa distribution combined with high initial densities (10$^4$-10$^6$ M$_{\odot}$ pc$^{-3}$), or low initial binary fractions ($\sim$ 5-10%) with a dominant hard binary component combined with moderate initial densities near their present-day values (10$^2$-10$^3$ M$_{\odot}$ pc$^{-3}$). This apparent degeneracy can be broken using the binary fractions outside r$_{\rm h}$- only high initial binary fractions with a significant soft component combined with high initial densities can contribute to reproducing the observed anti-correlation between the binary fractions outside r$_{\rm h}$ and the total cluster mass. We further illustrate using the simulated present-day binary orbital parameter distributions and the technique introduced in Leigh et al. (2012) that the relative fractions of hard and soft binaries can be used to further constrain the initial cluster density and mass-density relation. Our results favour an initial mass-density relation of the form r$_{\rm h} \propto$ M$_{\rm clus}^{\alpha}$ with $\alpha <$ 1/3.

Populations of multiple stars inside clustered regions are known to change through dynamical interactions. The efficiency of binary disruption is thought to be determined by stellar density. King and collaborators recently investigated the multiplicity properties in young star forming regions and in the Galactic field. They concluded that stellar density-dependent modification of a universal initial binary population (the standard or null hypothesis model) cannot explain the observations. We re-visit their results, analyzing the data within the framework of different model assumptions, namely non-universality without dynamical modification and universality with dynamics. We illustrate that the standard model does account for all known populations if regions were significantly denser in the past. Some of the effects of using present-day cluster properties as proxies for their past values are emphasized and that the degeneracy between age and density of a star forming region can not be omitted when interpreting multiplicity data. A new analysis of the Corona Australis region is performed within the standard model. It is found that this region is likely as unevolved as Taurus and an initial density of $\approx190M_{\odot}\;pc^{-3}$ is required to produce the presently observed binary population, which is close to its present-day density.

We present the results of a photometric search for variable stars in the core of the Galactic globular cluster M4. The input data are a large and unprecedented set of deep Hubble Space Telescope WFC3 images (large program GO-12911; 120 orbits allocated), primarily aimed at probing binaries with massive companions by detecting their astrometric wobbles. Though these data were not optimised to carry out a time-resolved photometric survey, their exquisite precision, spatial resolution and dynamic range enabled us to firmly detect 38 variable stars, of which 20 were previously unpublished. They include 19 cluster-member eclipsing binaries (confirming the large binary fraction of M4), RR Lyrae, and objects with known X-ray counterparts. We improved and revised the parameters of some among published variables.

Using a Monte Carlo code, we construct a dynamic evolutionary model of the Galactic globular cluster M22 (NGC6656). The initial conditions are chosen so that, after about 12Gyr of stellar and dynamical evolution, the model is an approximate fit to the surface brightness and velocity dispersion profiles of the cluster, to its mass function, and to the current binary fraction. Depending on the distribution of black hole natal kicks, we predict that the present-day population of stellar-mass black holes ranges from about 40 (no kicks) down to essentially zero (kicks distributed like those of neutron stars). Provided that natal kicks do not eject all new black holes, it is suggested that clusters with a present-day half-mass relaxation time above about 1Gyr are the ones that may still retain an appreciable population of black holes.

The M4 Core Project with HST is designed to exploit the Hubble Space Telescope to investigate the central regions of M4, the Globular Cluster closest to the Sun. In this paper we combine optical and near-infrared photometry to study multiple stellar populations in M4. We detected two sequences of M-dwarfs containing ~38% (MS_I) and ~62% (MS_II) of MS stars below the main-sequence (MS) knee. We compare our observations with those of NGC2808, which is the only other GCs where multiple MSs of very low-mass stars have been studied to date. We calculate synthetic spectra for M-dwarfs, assuming the chemical composition mixture inferred from spectroscopic studies of stellar populations along the red giant branch, and different Helium abundances, and we compare predicted and observed colors. Observations are consistent with two populations, one with primordial abundance and another with enhanced nitrogen and depleted oxygen.

We describe a major upgrade of a Monte Carlo code which has previously been used for many studies of dense star clusters. We outline the steps needed in order to calibrate the results of the new Monte Carlo code against $N$-body simulations for large $N$ systems, up to $N=200000$. The new version of the Monte Carlo code (called MOCCA), in addition to the features of the old version, incorporates the direct Fewbody integrator (Fregeau et al. 2004) for three- and four-body interactions, and a new treatment of the escape process based on Fukushige & Heggie (2000). Now stars which fulfil the escape criterion are not removed immediately, but can stay in the system for a certain time which depends on the excess of the energy of a star above the escape energy. They are called potential escapers. With the addition of the Fewbody integrator the code can follow all interaction channels which are important for the rate of creation of various types of objects observed in star clusters, and ensures that the energy generation by binaries is treated in a manner similar to the $N$-body model. There are at most three new parameters which have to be adjusted against $N$-body simulations for large $N$: two (or one, depending on the chosen approach) connected with the escape process, and one responsible for the determination of the interaction probabilities. The values adopted for the free parameters have at most a weak dependence on $N$. They allow MOCCA to reproduce $N$-body results with reasonable precision, not only for the rate of cluster evolution and the cluster mass distribution, but also for the detailed distributions of mass and binding energy of binaries. Additionally, the code can follow the rate of formation of blue stragglers and black hole - black hole binaries.

We present an overview of the ongoing Hubble Space Telescope large program GO-12911. The program is focused on the core of M4, the nearest Galactic globular cluster, and the observations are designed to constrain the number of binaries with massive companions (black holes, neutron stars, or white dwarfs) by measuring the ``wobble'' of the luminous (main-sequence) companion around the center of mass of the pair, with an astrometric precision of ~50 micro-arcseconds. The high spatial resolution and stable medium-band PSFs of WFC3/UVIS will make these measurements possible. In this work we describe: (i) the motivation behind this study, (ii) our observing strategy, (iii) the many other investigations enabled by this unique data set, and which of those our team is conducting, and (iv) a preliminary reduction of the first-epoch data-set collected on October 10, 2012.

We introduce an improved code for simulations of star clusters, called MOCCA. It combines the Monte Carlo method for star cluster evolution and the Fewbody code to perform scattering experiments. The Fewbody was added in order to track more precisely dynamical interactions between objects which can lead to creations of various exotic objects observed in the star clusters, like Blue Stragglers Stars (BSS). The MOCCA code is currently one of the most advanced codes for simulating real size star clusters. It follows the star cluster evolution closely to N-body codes but is much faster. We show that the MOCCA code is able to follow the evolution of BSS with details. It is a suitable tool to perform full scale evolution of real star clusters and detail comparison with observations of exotic star cluster objects like BSS. This paper is the first one of the series of papers about properties of BSS in star clusters. This type of stars is particularly interesting today, because by studying them one can get important constrains on a link between the stellar and dynamical evolution of star clusters. We discuss here first results concerning BSS for an arbitrary chosen test model. We investigate properties of BSS which characterize different channels of formation like masses, semi-major axes, eccentricities, and orbital periods. We show how BSS from different channels change their types, and discuss initial and final positions of BSS, their bimodal distribution in the star cluster, lifetimes and more.

We study the compact binary population in star clusters, focusing on binaries containing black holes, using a self-consistent Monte Carlo treatment of dynamics and full stellar evolution. We find that the black holes experience strong mass segregation and become centrally concentrated. In the core the black holes interact strongly with each other and black hole-black hole binaries are formed very efficiently. The strong interactions, however, also destroy or eject the black hole-black hole binaries. We find no black hole-black hole mergers within our simulations but produce many hard escapers that will merge in the galactic field within a Hubble time. We also find several highly eccentric black hole-black hole binaries that are potential LISA sources, suggesting that star clusters are interesting targets for space-based detectors. We conclude that star clusters must be taken into account when predicting compact binary population statistics.

We use a self-consistent Monte Carlo treatment of stellar dynamics to investigate black hole binaries that are dynamically ejected from globular clusters to determine if they will be gravitational wave sources. We find that many of the ejected binaries have initially short periods and will merge within a Hubble time due to gravitational wave radiation. Thus they are potential sources for ground-based gravitational wave detectors. We estimate the yearly detection rate for current and advanced ground-based detectors and find a modest enhancement over the rate predicted for binaries produced by pure stellar evolution in galactic fields. We also find that many of the ejected binaries will pass through the longer wavelength Laser Interferometer Space Antenna (LISA) band and may be individually resolvable. We find a low probability that the Galaxy will contain a binary in the LISA band during its three-year mission. Some such binaries may, however, be detectable at Mpc distances implying that there may be resolvable stellar-mass LISA sources beyond our Galaxy. We conclude that globular clusters have a significant effect on the detection rate of ground-based detectors and may produce interesting LISA sources in local group galaxies.

We describe Monte Carlo models for the dynamical evolution of the massive globular cluster 47 Tuc (NGC 104). The code includes treatments of two-body relaxation, most kinds of three- and four-body interactions involving primordial binaries and those formed dynamically, the Galactic tide, and the internal evolution of both single and binary stars. We arrive at a set of initial parameters for the cluster which, after 12Gyr of evolution, gives a model with a fairly satisfactory match to surface brightness and density profiles, the velocity dispersion profile, the luminosity function in two fields, and the acceleration of pulsars. Our models appear to require a relatively steep initial mass function for stars above about turnoff, with an index of about 2.8 (where the Salpeter mass function has an index of 2.35), and a relatively flat initial mass function (index about 0.4) for the lower main sequence. According to the model, the current mass is estimated at 0.9 million solar masses, of which about 34% consists of remnants. We find that primordial binaries are gradually taking over from mass loss by stellar evolution as the main dynamical driver of the core. Despite the high concentration of the cluster, core collapse will take at least another 20Gyr.

M4 and NGC 6397 are two very similar galactic globular clusters, which differ mainly in their surface brightness profile. M4 has a classic King-like profile, whereas NGC 6397 has a more concentrated profile, which is often interpreted as that of a post-core collapse cluster. Heggie & Giersz (2008), however, found that M4 is also a post-core collapse cluster, and Giersz & Heggie (2009) concluded that the main reason for the difference between the two surface brightness profiles is fluctuations. This conclusion was reached on the basis of Monte Carlo models, however, and in the present Letter we verify that similar fluctuations occur in N-body models. The models were initialised by generating initial conditions from the Monte Carlo model of NGC 6397 at the simulated age of 12Gyr, and one was followed for 1Gyr. The new models help to clarify the nature of the fluctuations, which have the nature of semi-regular oscillations with a time scale of order 10^8yr. They are influenced by the dynamical role which is played by primordial binaries in the evolution of the core.

We describe Monte Carlo models for the dynamical evolution of the nearby globular cluster NGC 6397. The code includes treatments of two-body relaxation, most kinds of three- and four-body interactions involving primordial binaries and those formed dynamically, the Galactic tide, and the internal evolution of both single and binary stars. We arrive at a set of initial parameters for the cluster which, after 12Gyr of evolution, gives a model with a fairly satisfactory match to the surface brightness profile, the velocity dispersion profile, and the luminosity function in two fields. We describe in particular those aspects of the evolution which distinguish this cluster from M4, which has a roughly similar mass and Galactocentric distance, but a qualitatively different surface brightness profile. Within the limitations of our modelling, we conclude that the most plausible explanation for the difference is fluctuations: both clusters are post-collapse objects, but sometimes have resolvable cores and sometimes not.

At least 10-15% of nearby sun-like stars have known Jupiter-mass planets. In contrast, very few planets are found in mature open and globular clusters such as the Hyades and 47 Tuc. We explore here the possibility that this dichotomy is due to the post-formation disruption of planetary systems associated with the stellar encounters in long-lived clusters. One supporting piece of evidence for this scenario is the discovery of freely floating low-mass objects in star forming regions. We use two independent numerical approaches, a hybrid Monte Carlo and a direct $N$-body method, to simulate the impact of the encounters. We show that the results of numerical simulations are in reasonable agreement with analytical determinations in the adiabatic and impulsive limits. They indicate that distant stellar encounters generally do not significantly modify the compact and nearly circular orbits. However, moderately close stellar encounters, which are likely to occur in dense clusters, can excite planets' orbital eccentricity and induce dynamical instability in systems which are closely packed with multiple planets. We compute effective cross sections for the dissolution of planetary systems and show that, for all initial eccentricities, dissolution occurs on time scales which are longer than the dispersion of small stellar associations, but shorter than the age of typical open and globular clusters. Although it is much more difficult to disrupt short-period planets, close encounters can excite modest eccentricity among them, such that subsequent tidal dissipation leads to orbital decay, tidal inflation, and even disruption of the close-in planets.

We describe Monte Carlo models for the dynamical evolution of the nearby globular cluster M4. The code includes treatments of two-body relaxation, three- and four-body interactions involving primordial binaries and those formed dynamically, the Galactic tide, and the internal evolution of both single and binary stars. We arrive at a set of initial parameters for the cluster which, after 12Gyr of evolution, gives a model with a satisfactory match to the surface brightness profile, the velocity dispersion profile, and the luminosity function in two fields. We describe in particular the evolution of the core, and find that M4 (which has a classic King profile) is actually a post-collapse cluster, its core radius being sustained by binary burning. We also consider the distribution of its binaries, including those which would be observed as photometric binaries and as radial-velocity binaries. We also consider the populations of white dwarfs, neutron stars, black holes and blue stragglers, though not all channels for blue straggler formation are represented yet in our simulations.

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 e.g. binary fraction, luminosity function and surface brightness. Comparison of Monte Carlo simulations with observational data for M67 shows that 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.

Astronomers have constructed models of globular clusters for over 100 years. These models mainly fall into two categories: (i) static models, such as King's model and its variants, and (ii) evolutionary models. Most attention has been given to static models, which are used to estimate mass-to-light ratios and mass segregation, and to combine data from proper motions and radial velocities. Evolutionary models have been developed for a few objects using the gaseous model, the Fokker-Planck model, Monte Carlo models and N-body models. These models have had a significant role in the search for massive black holes in globular clusters, for example. In this presentation the problems associated with these various techniques will be summarised, and then we shall describe new work with Giersz's Monte Carlo code, which has been enhanced recently to include the stellar evolution of single and binary stars. We describe in particular recent attempts to model the nearby globular cluster M4, including predictions on the spatial distribution of binary stars and their semi-major axis distribution, to illustrate the effects of about 12 Gyr of dynamical evolution. We also discuss work on an approximate way of predicting the "initial" conditions for such modelling.

We outline the steps needed in to calibrate the Monte Carlo code in order to perform large scale simulations of real globular clusters. We calibrate the results against $N$-body simulations for $N = 2500$, 10000 and for the old open cluster M67. The calibration is done by choosing appropriate free code parameters.

A revision of Stodolkiewicz'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. 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. During the evolution, the initial mass function (IMF) changes significantly. The local mass function (LMF) around the half--mass radius closely resembles the actual global mass function (GMF). At the late stages of evolution the mass of the evolved stars inside the core can be as high as 97% of the total mass in this region. For the whole system, the evolved stars can compose up to 67% 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.

We present new models of the evolution and dissolution of star clusters evolving under the combined influence of internal relaxation and external tidal fields, using the anisotropic gaseous model based on the Fokker-Planck approximation, and a new escaper loss cone model. This model borrows ideas from loss cones of stellar distributions near massive black holes, and describes physical processes related to escaping stars by a simple model based on two timescales and a diffusion process. We compare our results with those of direct $N$-body models and of direct numerical solutions of the orbit-averaged Fokker-Planck equation. For this comparative study we limit ourselves to idealized single point mass star clusters, in order to present a detailed study of the physical processes determining the rate of mass loss, core collapse and other features of the system's evolution. With the positive results of our study the path is now open in the future to use the computationally efficient gaseous models for future studies with more realism (mass spectrum, stellar evolution).

This is a summary paper of MODEST-2, a workshop held at the Astronomical Institute ``Anton Pannekoek'' in Amsterdam, 16-17 December 2002. MODEST is a loose collaboration of people interested in MOdelling DEnse STellar systems, particularly those interested in modelling these systems using all the available physics (stellar dynamics, stellar evolution, hydrodynamics and the interplay between the three) by defining interfaces between different codes. In this paper, we summarize 1) the main advances in this endeavour since MODEST-1; 2) the main science goals which can be and should be addressed by these types of simulations; and 3) the most pressing theoretical and modelling advances that we identified.

Spherically symmetric equal mass star clusters containing a large amount of primordial binaries are studied using a hybrid method, consisting of a gas dynamical model for single stars and a Monte Carlo treatment for relaxation of binaries and the setup of close resonant and fly-by encounters of single stars with binaries and binaries with each other (three- and four-body encounters). What differs from our previous work is that each encounter is being integrated using a highly accurate direct few-body integrator which uses regularized variables. Hence we can study the systematic evolution of individual binary orbital parameters (eccentricity, semi-major axis) and differential and total cross sections for hardening, dissolution or merging of binaries (minimum distance) from a sampling of several ten thousands of scattering events as they occur in real cluster evolution including mass segregation of binaries, gravothermal collapse and reexpansion, binary burning phase and ultimately gravothermal oscillations. For the first time we are able to present empirical cross sections for eccentricity variation of binaries in close three- and four-body encounters. It is found that a large fraction of three-body and four-body encounters results in merging. Previous cross sections obtained by Spitzer and Gao for strong encounters can be reproduced, while for weak encounters non-standard processes like formation of hierarchical triples occur.

We construct a simple dynamical model of the massive globular cluster Omega Cen. The model includes simple treatments of dynamical evolution in a galactic tide, and evolution of single stars. Binary stars and rotation are neglected. The model approximately fits observational data on the surface brightness profile, the profile of radial velocity dispersion, and the main sequence mass function at two radii.

A revision of Stod\{'o}{\l}kiewicz's Monte Carlo code is used to simulate evolution of large star clusters. 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 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, due to stellar evolution of the most massive stars, causes expansion of the whole cluster and eventually leads to the disruption of less bound systems ($W_0=3$). Models with larger $W_0$ 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.

A revision of Stod\'o{\l}kiewicz's Monte Carlo code is used to simulate evolution of large star clusters. The survey on the evolution of multi-mass N-body systems influenced by the tidal field of a parent galaxy and by stellar evolution is discussed. For the first time, the simulation on the "star-by-star" bases of evolution of 1,000,000 body star cluster is presented. \

The new approach outlined in Paper I (Spurzem \& Giersz 1996) 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 percent of the initial number of stars in the cluster). 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. In a self-consistent model it is the first time that such a realistically large number of binaries is evolving in a cluster with an even ten times larger number of single stars. Due to the Monte Carlo treatment of the binaries we can at every moment analyze their external and internal parameters in the cluster as in an N-body simulation.

A revision of Stodolkiewicz'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 10000 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 to 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 > 10000 gravothermal oscillations are clearly visible. The calculations of N=2000, 4096, 10000, 32000 and 100000 models take about 2, 6, 20, 130 and 2500 hours, respectively. The Monte-Carlo code is at least 10^5 times faster than the N-body one for N=32768 with special-purpose hardware (Makino 1996ab). 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 100000. 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.

We review recent progress in computational methods for studying the dynamical evolution of large star clusters. After describing some highlights of a new comparison of different models, including N-body models, we mention in detail recent developments in three codes which allow models with anisotropic distributions of velocity to be studied: (i) a Fokker-Planck code using finite differences, (ii) a Monte Carlo code and (iii) a gaseous model code.

We describe results from large numbers of $N$-body simulations containing from $250$ to $1000$ stars each. The distribution of stellar masses is a power law, and the systems are isolated. While the collapse of the core exhibits the expected segregation of different masses, we find that the post-collapse evolution is, at a first approximation, homologous. This is quite surprising because there is no reason for supposing that mass segregation should not continue to have a substantial effect on the evolution of the cluster. In fact the spatial distribution of the mean stellar mass is nearly static throughout the post-collapse regime, except for the overall expansion of the systems, and this helps to explain why the post-collapse evolution is nearly self-similar. Self-similarity is also exhibited by the distribution of anisotropy and the profile of departures from equipartition, which show little change during the post-collapse phase. The departures from energy equipartition and isotropy are small in the core and increase with radius. During post-collapse evolution massive stars (mainly) are removed from the system by binary activity. This effect dominates the preferential escape of low-mass stars due to standard two-body relaxation processes.

This paper presents and analyses statistical results from a large number of $N$-body simulations of isolated systems with equal masses, in which $250\le N\le 2000$. It concentrates on the phase starting around the end of core collapse, when binaries start to play a crucial role. Interactions of hard binaries are consistent with the Spitzer (1987) cross section. The evolution of the half-mass radius after core collapse nearly follows classical theory, The evolution of the inner parts of the system around the time of core bounce is consistent with these simplified models provided that the continuous production of energy, as usually assumed, is replaced by a model of {\sl stochastic} energy production. Similarly, post-collapse evolution of the core requires a modest recalibration of the coefficient of energy generation, especially for small $N$. The distributions of velocity dispersion and anisotropy become remarkably homologous soon after core bounce. The bound mass of the systems very nearly follows a power law with time. A small number of escapers, presumed to be those associated with binary activity, dominate the energy which is carried off. The ``internal" energy of escaping binaries is consistent with theoretical expectations.

We compare the results for the dynamical evolution of star clusters derived from anisotropic gaseous models with the data from N-body simulations of isolated and one-component systems, each having modest number of stars. The statistical quality of N-body data was improved by averaging results from many N-body runs, each with the same initial parameters but with different sequences of random numbers used to initialize positions and velocities of the particles. We study the development of anisotropy, the spatial evolution and energy generation by three-body binaries and its N-dependence. We estimate the following free parameters of anisotropic gaseous models: the time scale for collisional anisotropy decay and the coefficient in the formulae for energy generation by three-body binaries. To achieve a fair agreement between N-body and gaseous models for the core in pre- as well as in post-collapse only the energy generation by binaries had to be varied by N. We find that anisotropy has considerable influence on the spatial structure of the cluster particularly for the intermediate and outer regions.

We study the dynamical evolution of idealised stellar systems by averaging results from many $N$-body simulations, each having modest numbers of stars. For isolated systems with stars of uniform mass, we discuss aspects of evolution up to the point of core collapse: relaxation and its $N$-dependence, the evolution of the density profile, the development of the velocity dispersion and anisotropy, and the rate of stellar escape.